MODULE PrimeConstants EXPORTS Main;

(*
  Computes various series based on prime numbers.
  
  The prime numbers up to "MaxPrime" are enumerated in reverse order.
  The terms of the series are computed (and added) as intervals, using 
  directed IEEE FP rounding.
  
  The tail of the series is estimated independently, also as interval.
  The tail estimates are on the folleowing observation.
  Let 

      r be any integer >= 2, 

      L be the product p[1]p[2]..p[r], and

      M any integer >= p[r],
      
      K be the number of integers in {M+1..M+L-1} which are 
        relatively prime with L (i.e. indivisible by
        p[1],.. p[r]).
      
      a[1],.. a[k] be those integers, in increasing order
        
      p[N] the largest prime <= M
        
  Then, for any i>=0, any j in [1..K], we have

     p[N + j + K*i]  >=  M + a[j] + L*i
     
  That's because the primes ">= M" are a subset of the set S of the
  ``"r"-subprimes'', the integers undivisible by "p[1],.. p[r]".
  Hence prime "p[N+j+K*i]" is no less than the "(j+K*i)"th element of
  S greater than "M", which we denote by "t".  The set "S" has period
  "L", with "K" elements per period; it follows that "t = M+a[j]+L*i".

  Hence, if "f(i,p)" is nonincreasing in p, 

     sum(f(i,prime[i]), i=N+1..oo) =
     
     sum(sum(f(N+j+K*i, p[N+j+K*i]), i=0..oo), j=1..K) <=
     
     sum(sum(f(N+j+K*i, M+a[j]+L*i), i=0..oo), j=1..K) <=
     
     sum(int(f(N+j+K*x, M+a[j]+L*x), x=-1..oo), j=1..K)

*)

IMPORT Wr, Thread, Fmt, Process, ParseParams;
IMPORT FloatMode, LongFloat;
FROM Stdio IMPORT stdout, stderr;

TYPE
  BOOL = BITS 8 FOR BOOLEAN; BOOLS = ARRAY OF BOOL;
  NAT  = CARDINAL;           NATS  = ARRAY OF NAT;
  NUM  = LONGREAL;           NUMS  = ARRAY OF NUM;
 
TYPE
  Interval = RECORD lo, hi: NUM END;

  Options = RECORD
      minPrime: NAT;
      maxPrime: NAT;
      chunkSize: NAT;
      debug: BOOL;
    END;

PROCEDURE DoIt() =
  BEGIN
    WITH
      o = GetOptions(),
      maxLoPrime = 2 * ISqrt(o.maxPrime),
      loPrime = GenPrimes(maxLoPrime)^
    DO
      WITH pm = loPrime[LAST(loPrime)] DO 
        <* ASSERT o.maxPrime DIV pm <= pm *>
      END;
      ComputeConstant(PSqrTerm, PSqrTail, PSqrForm, o, loPrime);
      (* ComputeConstant(PCubTerm, PCubTail, PCubForm, o, loPrime); *)
      (* ComputeConstant(PQtrTerm, PQtrTail, PQtrForm, o, loPrime); *)
      (* ComputeConstant(PQuiTerm, PQuiTail, PQuiForm, o, loPrime); *)
      (* ComputeConstant(PSexTerm, PSexTail, PSexForm, o, loPrime); *)
      (* ComputeConstant(PPowTerm, PPowTail, PPowForm, o, loPrime); *)
      (* ComputeConstant(TestTerm, TestTail, TestForm, o, loPrime); *)
    END;
  END DoIt;
  
TYPE 
  TermProc = PROCEDURE (i, p: NAT): Interval;
    (*
      Must estimate the term "term(i,p)",
      assuming "p = prime[i]". *)
      
  TailProc = PROCEDURE (numP, maxP, K, L: NAT): Interval;
    (*
      Must estimate "sum(term(numP + K*i, prime[numP + K*i]), i=0..oo)"
      assuming that "prime[numP + K*i] >= maxP + L*i".
    *)
  
PROCEDURE ComputeConstant(
    term: TermProc;         (* The term estimator *)
    tail: TailProc;         (* The tail estimator *)
    form: TEXT;             (* The series form *)
    READONLY o: Options;    (* The options *)
    READONLY loPrime: NATS; (* The small prime number table *)
  ) =
  VAR nPrimes: CARDINAL;
      s: Interval;
  BEGIN
    Message("term = " & form);
    WITH
      sieveChunk = NEW(REF BOOLS, o.chunkSize)^,
      nChunks = o.maxPrime DIV o.chunkSize,
      c = NEW(REF ARRAY OF Interval, nChunks + 1)^
    DO
      nPrimes := 0;
      FOR ic := 0 TO nChunks-1 DO
        WITH
          minP = 1 + ic * o.chunkSize, 
          maxP = MIN(minP + o.chunkSize - 1, o.maxPrime)
        DO
          Message("Enumerating primes in range [" & FI(minP) & ".." & FI(maxP) & "]");
          ComputeSieveChunk(loPrime, minP, maxP, sieveChunk);
          Message("Computing partial sum...");
          c[ic] := ComputeChunkSum(
            term, nPrimes, o.minPrime, o.maxPrime,
            minP, maxP, sieveChunk, o.debug
          );
          Message(FI(nPrimes) & " primes found so far");
          PrintRange("c", c[ic].lo, c[ic].hi);
          Message("");
        END
      END;
      Message("Estimating tail...");
      c[nChunks] := ComputeTailEstimate(
        tail, nPrimes, o.maxPrime, 
        loPrime, sieveChunk, o.debug
      );
      Message("Total tail:");
      PrintRange("c", c[nChunks].lo, c[nChunks].hi);
      (* Now add the chunk ranges: *)
      s.lo := 0.0d0; s.hi := 0.0d0;
      FOR ic := nChunks-1 TO 0 BY -1 DO
        RDN(); s.lo := s.lo + c[ic].lo;
        RUP(); s.hi := s.hi + c[ic].hi;
      END;
      Message("Series without tail:");
      PrintRange("s", s.lo, s.hi);
      (* Add the tail: *)
      RDN(); s.lo := s.lo + c[nChunks].lo;
      RUP(); s.hi := s.hi + c[nChunks].hi;
      Message("Total series:");
      PrintRange("s", s.lo, s.hi);
      Message("");
    END
  END ComputeConstant;
  
PROCEDURE GenPrimes(maxPrime: NAT): REF NATS =
  VAR p, m, nPrimes: NAT;
  BEGIN
    Message("Generating primes up to " & FI(maxPrime));
    WITH
      isPrime = NEW(REF BOOLS, maxPrime+1)^
    DO
      (* Eratosthenes's sieve: *)
      isPrime[0] := FALSE; isPrime[1] := FALSE; isPrime[2] := TRUE;
      FOR i := 3 TO maxPrime DO isPrime[i] := (i MOD 2 = 1) END;
      p := 3; m := 9;
      WHILE m <= maxPrime DO
        IF isPrime[p] THEN
          REPEAT isPrime[m] := FALSE; m := m + p UNTIL m > maxPrime
        END;
        p := p + 2;
        m := p*p
      END;
      (* Count primes: *)
      nPrimes := 1;
      FOR p := 3 TO maxPrime BY 2 DO
        IF isPrime[p] THEN INC(nPrimes) END
      END;
      (* Gather prime number table: *)
      WITH
        rp = NEW(REF NATS, nPrimes),
        p = rp^
      DO
        p[0] := 2; nPrimes := 1; 
        FOR k := 3 TO maxPrime BY 2 DO
          IF isPrime[k] THEN p[nPrimes] := k; INC(nPrimes) END
        END;
        PrintPrimeData(p);
        RETURN rp
      END
    END
  END GenPrimes;
  
PROCEDURE ComputeSieveChunk(
    READONLY loPrime: NATS;   (* Smapp prime number table *)
    minP, maxP: NAT;          (* Range of prime numbers desired *)
    VAR s: BOOLS;             (* The sieve chunk *)
  ) =
  (*
    Sets "s[k] := TRUE" if and only if "minP + k" is prime,
    for "k" in the range "[0..maxP-minP]".
    The table "loPrime" must contain all primes whose
    square is not greater than "maxP".
  *)
  VAR k, m: NAT;
  BEGIN
    (* Initialize *)
    FOR i := minP TO MAX(2, minP)-1 DO s[i - minP] := FALSE END;
    FOR i := MAX(2, minP) TO maxP DO s[i - minP] := TRUE END;
    (* Eliminate multiples of small primes *)
    k := 0; 
    LOOP
      WITH 
        p = loPrime[k],
        dMax = maxP DIV p
      DO
        IF dMax < p THEN EXIT END;
        m := MAX((minP + p - 1) DIV p, p) * p;
        WHILE m <= maxP DO s[m - minP] := FALSE; m := m + p END;
      END;
      INC(k);
    END
  END ComputeSieveChunk;
  
PROCEDURE ComputeChunkSum(
    term: TermProc;        (* The term function. *)
    VAR nPrimes: NAT;      (* IN: "#{primes < minP}"; OUT: "#{primes < maxP}". *)
    minPrime: NAT;         (* Lower bound for primes to include in sum. *)
    maxPrime: NAT;         (* Upper bound for primes to include in sum. *)
    minP, maxP: NAT;       (* Range of sieve chunk. *)
    READONLY s: BOOLS;     (* The sieve for integers in "[minP..maxP]". *)
    debug: BOOL;           (* TRUE to print each term. *)
  ): Interval = 
  (* 
    Computes the partial sum of "term(i, p[i])" for 
    the primes "p[i]" in the range "[minP..maxP]", given the
    sieve for that range.
  *)
  VAR lo, hi: NUM; i, nTerms, nChunkPrimes: NAT;
  BEGIN
    lo := 0.0d0; hi := 0.0d0;
    (* Count primes: *)
    nChunkPrimes := 0;
    FOR k := 0 TO maxP - minP DO IF s[k] THEN INC(nChunkPrimes) END END;
    nPrimes := nPrimes + nChunkPrimes;
    (* Add terms, in reverse order: *)
    WITH 
      minT = MAX(minP, minPrime), 
      maxT = MIN(maxP, maxPrime)
    DO
      i := nPrimes;
      FOR p := maxP TO maxT+1 BY -1 DO IF s[p-minP] THEN DEC(i) END END;
      nTerms := 0;
      FOR p := maxT TO minT BY -1 DO
        IF s[p-minP] THEN
          WITH t = term(i, p) DO
            IF debug THEN PrintTerm(i, p, t) END;
            RDN(); lo := lo + t.lo;
            RUP(); hi := hi + t.hi;
          END;
          DEC(i); INC(nTerms)
        END
      END
    END;
    Message(FI(nChunkPrimes) & " primes in chunk and " & FI(nTerms) & " terms added");
    RETURN Interval{lo, hi}
  END ComputeChunkSum;

PROCEDURE ComputeTailEstimate(
    tail: TailProc; 
    nPrimes: NAT;            (* Number of primes included in actual sum. *)
    maxPrime: NAT;           (* Max prime included in actual sum *)
    READONLY loPrime: NATS;  (* Table of small primes *)
    VAR s: BOOLS;            (* Work area for "r"-subprime sieve *)
    debug: BOOL;             (* TRUE to print each tail term *)
  ): Interval =
  VAR r, L, K, b, p, j: NAT;
      lo, hi: NUM;
  BEGIN
    (* Compute the period length of the subprimes: *)
    L := 1; r := 0;
    WHILE L < NUMBER(s) DIV loPrime[r] DO L := L*loPrime[r]; INC(r) END;
    (* Compute the sieve for the first period greater than maxPrime: *)
    b := (maxPrime + 1) MOD L; (* s[a] corresponds to integers (b+a) + L*i *)
    FOR a := 0 TO L-1 DO s[a] := TRUE END;
    FOR t := 0 TO r-1 DO
      p := loPrime[t];
      (* find first a such that b+a is a multiple of p: *)
      FOR a := (-b) MOD p TO L-1 BY p DO s[a] := FALSE  END;
    END;
    (* Count primes per period: *)
    K := 0; FOR a := 0 TO L-1 DO IF s[a] THEN INC(K) END END;
    (* Add tail terms: *)
    lo := 0.0d0; hi := 0.0d0;
    j := K;
    FOR a := L-1 TO 0 BY -1 DO
      IF s[a] THEN 
        WITH t = tail(nPrimes+j, maxPrime+1+a, K, L) DO
          IF debug THEN PrintTailTerm(nPrimes+j, maxPrime+1+a, K, L, t) END;
          RDN(); lo := lo + t.lo;
          RUP(); hi := hi + t.hi;
        END;
        DEC(j);
      END
    END;
    RETURN Interval{lo, hi}
  END ComputeTailEstimate;
  
PROCEDURE RUP() = 
  <* FATAL FloatMode.Failure *>
  BEGIN 
    FloatMode.SetRounding(FloatMode.RoundingMode.TowardPlusInfinity)
  END RUP;
  
PROCEDURE RDN() = 
  <* FATAL FloatMode.Failure *>
  BEGIN 
    FloatMode.SetRounding(FloatMode.RoundingMode.TowardMinusInfinity)
  END RDN;
  
PROCEDURE ISqrt(n: NAT): NAT =
  VAR r, m: NAT := 0;
  BEGIN
    IF n = 0 THEN
      RETURN 0
    ELSE
      r := 1;
      REPEAT m := r; r := (r + n DIV r) DIV 2 UNTIL m = r;
      RETURN r
    END;
  END ISqrt;

PROCEDURE PrintPrimeData(READONLY prime: NATS) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    WITH N = NUMBER(prime) DO
      Wr.PutText(stdout, "prime[" & Fmt.Int(N-1) & "] = " & Fmt.Int(prime[N-1]) & "\n"); 
    END;
    Wr.Flush(stdout)
  END PrintPrimeData;

PROCEDURE PrintTerm(i, p: NAT; READONLY t: Interval) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(stdout, "term(" & FI(i) & "," & FI(p) & ") = \n");
    PrintRange("t", t.lo, t.hi)
  END PrintTerm;

PROCEDURE PrintTailTerm(numP, maxP, K, L: NAT; READONLY t: Interval) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(stdout, 
      "tail term(" & FI(numP) & "," & FI(maxP) & 
      "," & FI(K) & "," & FI(L) & ") = \n"
    );
    PrintRange("t", t.lo, t.hi)
  END PrintTailTerm;

PROCEDURE PrintRange(tag: TEXT; lo, hi: NUM) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(stdout, tag & ".lo = " & FLR(lo, 30) & "\n");
    Wr.PutText(stdout, tag & ".hi = " & FLR(hi, 30) & "\n");
    Wr.Flush(stdout)
  END PrintRange;

PROCEDURE Message(txt: TEXT) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(stdout, txt & "\n"); 
    Wr.Flush(stdout)
  END Message;

PROCEDURE FLR(x: NUM; w: NAT): TEXT =
  BEGIN
    RETURN Fmt.LongReal(x, style := Fmt.Style.Fix, prec := w)
  END FLR;
  
PROCEDURE FI(x: NAT): TEXT =
  BEGIN
    RETURN Fmt.Int(x)
  END FI;
  
PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>

  CONST
    MaxMaxPrime = 1024*1024*1024 - 1; (* 2^30-1 *)
    MaxChunkSize = 16*1024*1024;      (* 16 * 2^20 *)

  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO

      TRY
        IF pp.keywordPresent("-minPrime") THEN
          o.minPrime := pp.getNextInt(1, MaxMaxPrime)
        ELSE
          o.minPrime := 1
        END;

        pp.getKeyword("-maxPrime");
        o.maxPrime := pp.getNextInt(o.minPrime, MaxMaxPrime);
        
        IF pp.keywordPresent("-chunkSize") THEN
          o.chunkSize := pp.getNextInt(1, MaxChunkSize)
        ELSE
          o.chunkSize := o.maxPrime
        END;
        
        o.debug := pp.keywordPresent("-debug");
        
        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: PrimeConstants \\\n");
          Wr.PutText(stderr, "  -maxPrime NNN [ -chunkSize NNN ] \n");
          Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;
  
(* ======= INVERSE PRIME SQUARE SERIES *)

CONST PSqrForm = "1/p^2";

PROCEDURE PSqrTerm(<*UNUSED*> i: NAT; p: NAT): Interval =
  VAR lo, hi: NUM;
  BEGIN
    WITH fp = FLOAT(p, NUM) DO
      RDN(); lo := (1.0d0/fp)/fp;
      RUP(); hi := (1.0d0/fp)/fp;
    END;
    RETURN Interval{lo, hi}
  END PSqrTerm;

PROCEDURE PSqrTail(
    <*UNUSED*> numP: NAT; 
    maxP: NAT;
    <*UNUSED*> K: NAT;
    L: NAT;
  ): Interval =
  VAR lo, hi: NUM;
  BEGIN
    (* 
      Lower boudn is zero.
      Upper bound is 
        int((maxP+L*x)^(-2), x=-1..oo) = 
        int((u)^(-2)/L, u=maxP-L..oo) =
        (maxP-L)^(-1)/(L)
    *)
    WITH 
      fp = FLOAT(maxP-L, NUM),
      fl = FLOAT(L, NUM)
    DO
      lo := 0.0d0;
      RUP(); hi := (1.0d0/fp)/fl;
    END;
    RETURN Interval{lo, hi}
  END PSqrTail;

(* ======= INVERSE PRIME CUBE SERIES *)

CONST PCubForm = "1/p^3";

PROCEDURE PCubTerm(<*UNUSED*> i: NAT; p: NAT): Interval =
  VAR lo, hi: NUM;
  BEGIN
    WITH fp = FLOAT(p, NUM) DO
      RDN(); lo := ((1.0d0/fp)/fp)/fp;
      RUP(); hi := ((1.0d0/fp)/fp)/fp;
    END;
    RETURN Interval{lo, hi}
  END PCubTerm;

PROCEDURE PCubTail(
    <*UNUSED*> numP: NAT; 
    maxP: NAT;
    <*UNUSED*> K: NAT;
    L: NAT;
  ): Interval =
  VAR lo, hi: NUM;
  BEGIN
    (* 
      Lower boudn is zero.
      Upper bound is 
        int((maxP+L*x)^(-3), x=-1..oo) = 
        int((u)^(-3)/L, u=maxP-L..oo) =
        (maxP-L)^(-2)/(2*L)
    *)
    WITH 
      fp = FLOAT(maxP-L, NUM),
      fl = FLOAT(2*L, NUM)
    DO
      lo := 0.0d0;
      RUP(); hi := ((1.0d0/fp)/fp)/fl;
    END;
    RETURN Interval{lo, hi}
  END PCubTail;

(* ======= INVERSE PRIME QARTIC SERIES *)

CONST PQtrForm = "1/p^4";

PROCEDURE PQtrTerm(<*UNUSED*> i: NAT; p: NAT): Interval =
  VAR lo, hi: NUM;
  BEGIN
    WITH fp = FLOAT(p, NUM) DO
      RDN(); lo := (1.0d0/fp); lo := lo*lo; lo := lo*lo;
      RUP(); hi := (1.0d0/fp); hi := hi*hi; hi := hi*hi;
    END;
    RETURN Interval{lo, hi}
  END PQtrTerm;

PROCEDURE PQtrTail(
    <*UNUSED*> numP: NAT; 
    maxP: NAT;
    <*UNUSED*> K: NAT;
    L: NAT;
  ): Interval =
  VAR lo, hi: NUM;
  BEGIN
    (* 
      Lower boudn is zero.
      Upper bound is 
        int((maxP+L*x)^(-4), x=-1..oo) = 
        int((u)^(-4)/L, u=maxP-L..oo) =
        (maxP-L)^(-3)/(3*L)
    *)
    WITH 
      fp = FLOAT(maxP-L, NUM),
      fl = FLOAT(3*L, NUM)
    DO
      lo := 0.0d0;
      RUP(); hi := (((1.0d0/fp)/fp)/fp)/fl;
    END;
    RETURN Interval{lo, hi}
  END PQtrTail;

(* ======= INVERSE PRIME QUINTIC SERIES *)

CONST PQuiForm = "1/p^5";

PROCEDURE PQuiTerm(<*UNUSED*> i: NAT; p: NAT): Interval =
  VAR lo, hi: NUM;
  BEGIN
    WITH fp = FLOAT(p, NUM) DO
      RDN(); lo := (1.0d0/fp); lo := lo*lo; lo := lo*lo; lo := lo/fp;
      RUP(); hi := (1.0d0/fp); hi := hi*hi; hi := hi*hi; hi := hi/fp;
    END;
    RETURN Interval{lo, hi}
  END PQuiTerm;

PROCEDURE PQuiTail(
    <*UNUSED*> numP: NAT; 
    maxP: NAT;
    <*UNUSED*> K: NAT;
    L: NAT;
  ): Interval =
  VAR lo, hi: NUM;
  BEGIN
    (* 
      Lower boudn is zero.
      Upper bound is 
        int((maxP+L*x)^(-5), x=-1..oo) = 
        int((u)^(-5)/L, u=maxP-L..oo) =
        (maxP-L)^(-4)/(4*L)
    *)
    WITH 
      fp = FLOAT(maxP-L, NUM),
      fl = FLOAT(4*L, NUM)
    DO
      lo := 0.0d0;
      RUP(); hi := (1.0d0/fp); hi := hi*hi; hi := hi*hi; hi := hi/fl;
    END;
    RETURN Interval{lo, hi}
  END PQuiTail;

(* ======= INVERSE PRIME SEXTIC SERIES *)

CONST PSexForm = "1/p^6";

PROCEDURE PSexTerm(<*UNUSED*> i: NAT; p: NAT): Interval =
  VAR lo, hi: NUM;
  BEGIN
    WITH fp = FLOAT(p, NUM) DO
      RDN(); lo := (1.0d0/fp); lo := lo*lo; lo := lo/fp; lo := lo*lo;
      RUP(); hi := (1.0d0/fp); hi := hi*hi; hi := hi/fp; hi := hi*hi;
    END;
    RETURN Interval{lo, hi}
  END PSexTerm;

PROCEDURE PSexTail(
    <*UNUSED*> numP: NAT; 
    maxP: NAT;
    <*UNUSED*> K: NAT;
    L: NAT;
  ): Interval =
  VAR lo, hi: NUM;
  BEGIN
    (* 
      Lower boudn is zero.
      Upper bound is 
        int((maxP+L*x)^(-6), x=-1..oo) = 
        int((u)^(-6)/L, u=maxP-L..oo) =
        (maxP-L)^(-5)/(5*L)
    *)
    WITH 
      fp = FLOAT(maxP-L, NUM),
      fl = FLOAT(5*L, NUM)
    DO
      lo := 0.0d0;
      RUP(); hi := (1.0d0/fp); hi := hi*hi; hi := hi*hi; hi := (hi/fp)/fl;
    END;
    RETURN Interval{lo, hi}
  END PSexTail;

(* ======= INVERSE PRIME POWER SERIES FOR TESTING *)

CONST PPowForm = "1/2^(-p)";

PROCEDURE PPowTerm(<*UNUSED*> i: NAT; p: NAT): Interval =
  VAR lo, hi: NUM;
  <* FATAL FloatMode.Trap *>
  BEGIN
    IF p > 200 THEN
      lo := 0.0d0; 
      RUP(); hi := LongFloat.Scalb(1.0d0, -200)
    ELSE
      RDN(); lo := LongFloat.Scalb(1.0d0, -p);
      RUP(); hi := LongFloat.Scalb(1.0d0, -p);
    END;
    RETURN Interval{lo, hi}
  END PPowTerm;

PROCEDURE PPowTail(
    <*UNUSED*> numP: NAT; 
    maxP: NAT;
    <*UNUSED*> K: NAT;
    L: NAT;
  ): Interval =
  VAR lo, hi: NUM;
  <* FATAL FloatMode.Trap *>
  BEGIN
    (* 
      Lower bound is zero.
      Upper bound is
        sum(2^-(maxP+L*i), i=0..oo) = 
        2^(-maxP)*sum((2^-L)^i, i=0..oo) <=
        2^(-maxP)*(1 + 2^-(L-1))
    *)
    lo := 0.0d0; 
    (* Upper bound for (1 + 2^-(L-1)): *)
    IF L > 200 THEN
      RUP(); hi := 1.0d0 + LongFloat.Scalb(1.0d0, -199);
    ELSE
      RUP(); hi := 1.0d0 + LongFloat.Scalb(1.0d0, -L+1);
    END;
    (* Now scale it: *)
    IF maxP > 200 THEN
      RUP(); hi := LongFloat.Scalb(hi, -200)
    ELSE
      RUP(); hi := LongFloat.Scalb(hi, -maxP);
    END;
    RETURN Interval{lo, hi}
  END PPowTail;

(* ======= TRIVIAL SERIES FOR TESTING *)

CONST TestForm = "1/3 if p = 3 else 0";

PROCEDURE TestTerm(<*UNUSED*> i: NAT; p: NAT): Interval =
  VAR lo, hi: NUM;
  BEGIN
    IF p = 3 THEN
      WITH fp = FLOAT(p, NUM) DO
        RDN(); lo := 1.0d0/fp;
        RUP(); hi := 1.0d0/fp
      END
    ELSE
      lo := 0.0d0; hi := 0.0d0
    END;
    RETURN Interval{lo, hi}
  END TestTerm;

PROCEDURE TestTail(
    <*UNUSED*> numP: NAT; 
    <*UNUSED*> maxP: NAT;
    <*UNUSED*> K: NAT;
    <*UNUSED*> L: NAT;
  ): Interval =
  VAR lo, hi: NUM;
  BEGIN
    lo := 0.0d0; hi := 0.0d0;
    RETURN Interval{lo, hi}
  END TestTail;

(* =================================================== *)

BEGIN
  DoIt();
END PrimeConstants.