PrecisionStrategy.mag

// This file is part of ExactpAdics
// Copyright (C) 2018 Christopher Doris
//
// ExactpAdics is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// ExactpAdics is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with ExactpAdics.  If not, see <http://www.gnu.org/licenses/>.


///# Precision strategies
///
/// All intrinsics in `ExactpAdics` either return the correct value or raise an error if this is not possible. In particular, some p-adic computations are not possible unless the inputs are known to sufficient precision, and if this occurs then we raise a *precision error*. Precision strategies are how we tell these computation how much we are willing to increase certain precisions before giving up.
/// 
///toc
///
///## Motivation
/// 
/// Consider we have a value `x` which is known to be 0 modulo 2:
/// 
/// ```
/// > K := ExactpAdicField(2);
/// > x := K!(2^19+1) - K!1;
/// > print x;
/// O(2)
/// ```
/// 
/// Suppose we want its valuation. All we know at the moment is that it is at least 1. We could try increasing the precision of `x` to 10 and find:
/// 
/// ```
/// > IncreaseAbsolutePrecision(x, 10);
/// > print x;
/// O(2^10)
/// ```
/// 
/// So its valuation is at least 10. Now try increasing its precision to 20:
/// 
/// ```
/// > IncreaseAbsolutePrecision(x, 20);
/// > print x;
/// 2^19 + O(2^20)
/// ```
/// 
/// And now we know the valuation is 19. But it was possible that the valuation was 100 or 1000 or Inifinity, so how often should we try increasing the precision, and by how much, before giving up? This decision is left to the user, and this information is encapsulated in a *precision strategy* which is essentially an ascending sequence of numbers. Most functions which use such strategies do so by increasing their inputs to each number in the sequence in turn until the computation can succeed.
/// 
/// ```
/// > K := ExactpAdicField(2);
/// > x := K!(2^19+1) - K!1;
/// > print x;
/// O(2)
/// > Valuation(x : Strategy := [1,5,10]);
/// Runtime error: weakly zero
/// > print x;
/// O(2^10)
/// > Valuation(x : Strategy := [1,5,10,20]);
/// 19
/// > print x;
/// 2^19 + O(2^20)
/// > Valuation(x : Stragegy := [1,5,10]);
/// 19
/// > print x;
/// 2^19 + O(2^10)
/// ```
/// 
/// This example shows that asking for the valuation of something with a strategy which doesn't go high enough will result in a precision error. And once the valuation is known, it remains known, even if we ask for it with a strategy which appears too small later.
/// 
///## Conventions
/// 
/// In this package we obey the convention that any intrinsic which accepts precision strategies will:
/// - Accept the strategies via named parameters including the word `Strategy`; e.g. `Strategy`, `MainStrategy`, `fStrategy`.
/// - Have a `Strategy` parameter which acts as a default value for any other parameters taking strategies.
/// - Use the named global strategy `"default"` as the default strategy.
/// 
/// For example, this is a valid intrinsic:
/// 
/// ```
/// intrinsic Foo(x : Strategy:="default", OtherStrategy:=Strategy) -> .
/// ```
/// 
/// Functions which accept strategies should be careful to propagate these out to any function calls which accept strategies.
/// 
///## Format of a strategy
/// 
/// A strategy is a strictly increasing sequence of integers, and is represented by one of the following:
/// - An integer `RngIntElt`, which is interpreted as a strategy of length 1.
/// - A sequence `SeqEnum` or list `List` of strategies, which is interpreted as the concatenation of those strategies.
/// - A string `MonStgElt`, which references a [globaly defined strategy](#global-strategies).
/// - A function `UserProgram`, which is called repeatedly on the last value. It either returns false, signifying the end of the strategy, or returns true and a new strategy. For example the strategy `func<n | n lt 100, 2*n>` keeps doubling the precision up to a maximum value of 100.
/// - An error `Err`, which causes that error to be raised.
/// - A tuple `Tup` which is a special instruction of one of the following forms:
///   - `<"limit", n>` which applies a limit to the rest of the strategy; that is, the first value in the strategy greater or equal to `n` is replaced by `n`, and then the strategy terminates.
///   - `<"randomize">` which randomizes the strategy as follows: if `m` was the last value in the strategy and `n` is the next value, then replace `n` by `Random(m+1,n)`.
///   - `<"exp", m>` which is shorthand for the strategy `func<n | true, Max(1, Ceiling(n*m))>` which multiplies the previous value by `m` and rounds up. `m` need not be an integer.
///   - `<"double">` which is shorthand for `<"exp", 2>`.


import "ExactpAdics.mag": IS_INSTANCE;

STATE := recformat<stack, prec, limit, randomize>;

declare verbose ExactpAdics_ExecutePrecisionStrategy, 1;

///## Execution

intrinsic ExactpAdics_StartPrecisionStrategy(strategy : InitialPrecision:=0) -> .
  {An object representing the initial state of the strategy.}
  return rec<STATE | stack:=[*strategy*], prec:=InitialPrecision, limit:=Infinity(), randomize:=false>;
end intrinsic;

intrinsic ExactpAdics_StepPrecisionStrategy(~ok, ~pr, ~state)
  {If the strategy can be stepped further, sets ok to true, sets pr to the next precision and sets state to the next state. Otherwise sets ok to false, pr is unchanged and state may change but can still not be stepped further.}

  // limit breached
  if state`prec ge state`limit then
    ok := false;
    return;
  end if;

  // keep popping stuff off the stack until done
  while #state`stack gt 0 do
    
    // pop a strategy off the stack
    strat := state`stack[#state`stack];
    state`stack := state`stack[1..#state`stack-1];

    // its type determines how to interpret it
    case Type(strat):

    // INTEGER: is a precision to return
    when RngIntElt:
      if strat le state`prec then
        error "precision did not increase";
      end if;
      state`prec := Min(state`randomize select Random(state`prec+1, strat) else strat, state`limit);
      ok := true;
      pr := state`prec;
      return;

    // LIST: add to the stack
    when List:
      state`stack cat:= Reverse(strat);
      continue;

    // SEQUENCE: add to the stack
    when SeqEnum:
      state`stack cat:= [*x : x in Reverse(strat)*];
      continue;

    // STRING: use the corresponding global named strategy
    when MonStgElt:
      Append(~state`stack, ExactpAdics_GetGlobalPrecisionStrategy(strat));
      continue;

    // FUNCTION: call it; if it returns false then do nothing, else it returns true and another strategy; push the new strategy and the function back on the stack
    when UserProgram:
      ok, s := strat(state`prec);
      if ok then
        state`stack cat:= [* strat, s *];
      end if;
      continue;

    // TUPLE: special instructions
    when Tup:
      if #strat ge 1 and Type(strat[1]) eq MonStgElt then
        case strat[1]:
        when "limit":
          if (#strat eq 2)
          and IS_INSTANCE(strat[2], {RngIntElt,Infty,FldRatElt,FldReElt}) then
            state`limit := strat[2];
            continue;
          end if;
        when "randomize":
          if #strat eq 1 then
            state`randomize := true;
            continue;
          elif #strat eq 2 and IS_INSTANCE(strat[2], BoolElt) then
            state`randomize := strat[2];
            continue;
          end if;
        when "exp":
          if (#strat ge 2)
          and (#strat le 3)
          and IS_INSTANCE(strat[2], {RngIntElt,FldRatElt,FldReElt})
          and strat[2] gt 1
          and ((#strat lt 3) or IS_INSTANCE(strat[3], RngIntElt)) then
            m := strat[2];
            base := (#strat lt 3) select 0 else strat[3];
            Append(~state`stack, func<pr | true, base + Max(1, Ceiling((pr-base)*m))>);
            continue;
          end if;
        when "double":
          Append(~state`stack, <"exp", 2>);
          continue;
        end case;
      end if;

    // ERROR: raise it (only useful for debugging)
    when Err:
      error strat;
      
    end case;
    error "bad strategy:", strat;

  end while;

  // stack empty
  ok := false;

end intrinsic;

intrinsic ExactpAdics_ExecutePrecisionStrategy(cb :: UserProgram, strategy : InitialPrecision:=0) -> BoolElt, RngIntElt, .
  {Executes the given precision strategy: calls cb(pr) for each pr in the strategy; if it ever returns (true,value) then this returns (true,pr,value), otherwise it returns (false,pr,_) where pr is the last precision used in the strategy.}
  st := ExactpAdics_StartPrecisionStrategy(strategy : InitialPrecision:=InitialPrecision);
  while true do
    ExactpAdics_StepPrecisionStrategy(~ok, ~pr, ~st);
    if not ok then
      return false, pr, _;
    end if;
    done, result := cb(pr);
    if done then
      if assigned result then
        return true, pr, result;
      else
        return true, pr, _;
      end if;
    end if;
  end while;
end intrinsic;

intrinsic ExactpAdics_ExecutePrecisionStrategy_Lazy(strategy, getDeps :: UserProgram, getValue :: UserProgram) -> ExactpAdics_Gettr
  {The (general) getter which, for each pr in the strategy calls getDeps(pr, ~getter) to get a dependent getter, which evaluates to gval, then calls getValue(gval, ~val) to get the value.}
  st := ExactpAdics_StartPrecisionStrategy(strategy);
  ExactpAdics_StepPrecisionStrategy(~ok, ~pr, ~st);
  if not ok then
    return ExactpAdics_ConstGetter([*false,pr,false*]);
  end if;
  return ExactpAdics_GeneralGetter([*st,pr*],
    procedure (~stpr, ~getter)
      getDeps(stpr[2], ~getter);
    end procedure,
    procedure (gval, ~stpr, ~val)
      getValue(gval, ~val2);
      if not assigned val2 then
        ExactpAdics_StepPrecisionStrategy(~ok, ~stpr[2], ~stpr[1]);
        if not ok then
          val := [*false, stpr[2], false*];
        end if;
      else
        val := [*true, stpr[2], val2*];
      end if;
    end procedure);
end intrinsic;

///## Global strategies
/// 
/// For user convenience, there is a table of globally defined strategies which may be referred to by name. There are 3 strategies provided by default, which may all be overridden:
/// - `"defaultLimit"`: Initially set to `<"limit", 100>`, this strategy doesn't have any values, but may be mixed in to other strategies to limit them.
/// - `"unlimitedDefault"`: Initially set to `[* 1, <"randomize">, <"double"> *]`, that is it starts at 1, and for each value `n`, the next value is `Random(n+1,2*n)`.
/// - `"default"`: Initially set to `[* "defaultLimit", "unlimitedDefault" *]`, that is it is the `unlimitedDefault` strategy with the limit imposed by `defaultLimit`, hence with the other defaults it starts at 1, increments `n` to `Random(n+1,2*n)` and stops after this gets greater than 100.
/// 
/// The string `"default"` is used as the default `Strategy` parameter by most instrinsics in this package, so for most users it suffices to simply set this global strategy to their preferred one and then use the intrinsics without any explicit `Strategy` parameters.

STRATS := NewStore();
StoreSet(STRATS, "defaultLimit", <"limit",100>);
StoreSet(STRATS, "unlimitedDefault", [* 1, <"randomize">, <"double"> *]);
StoreSet(STRATS, "default", ["defaultLimit", "unlimitedDefault"]);

intrinsic ExactpAdics_IsGlobalPrecisionStrategyDefined(name :: MonStgElt) -> BoolElt, .
  {Returns true if there is a global strategy with the given name, and if so, returns the strategy itself.}
  return StoreIsDefined(STRATS, name);
end intrinsic;

intrinsic ExactpAdics_SetGlobalPrecisionStrategy(name :: MonStgElt, strat)
  {Defines the global strategy with the given name to the given strategy.}
  StoreSet(STRATS, name, strat);
end intrinsic;

intrinsic ExactpAdics_GetGlobalPrecisionStrategy(name :: MonStgElt) -> .
  {Returns the global strategy with the given name.}
  ok, s := ExactpAdics_IsGlobalPrecisionStrategyDefined(name);
  error if not ok, "global exact p-adic precision strategy not defined:", name;
  return s;
end intrinsic;