--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--Copyright (C) 2007-2010, Gabriel Dos Reis.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

)abbrev domain FLOAT Float

B ==> Boolean
I ==> Integer
S ==> String
PI ==> PositiveInteger
RN ==> Fraction Integer
SF ==> DoubleFloat
N ==> NonNegativeInteger

++ Author: Michael Monagan
++ Date Created:
++   December 1987
++ Change History:
++   19 Jun 1990
++ Basic Operations: outputFloating, outputFixed, outputGeneral, outputSpacing,
++   atan, convert, exp1, log2, log10, normalize, rationalApproximation,
++   relerror, shift, / , **
++ Keywords: float, floating point, number
++ Description: \spadtype{Float} implements arbitrary precision floating
++ point arithmetic.
++ The number of significant digits of each operation can be set
++ to an arbitrary value (the default is 20 decimal digits).
++ The operation \spad{float(mantissa,exponent,\spadfunFrom{base}{FloatingPointSystem})} for integer
++ \spad{mantissa}, \spad{exponent} specifies the number
++ \spad{mantissa * \spadfunFrom{base}{FloatingPointSystem} ** exponent}
++ The underlying representation for floats is binary
++ not decimal. The implications of this are described below.
++
++ The model adopted is that arithmetic operations are rounded to
++ to nearest unit in the last place, that is, accurate to within
++ \spad{2**(-\spadfunFrom{bits}{FloatingPointSystem})}.
++ Also, the elementary functions and constants are
++ accurate to one unit in the last place.
++ A float is represented as a record of two integers, the mantissa
++ and the exponent.  The \spadfunFrom{base}{FloatingPointSystem}
++ of the representation is binary, hence
++ a \spad{Record(m:mantissa,e:exponent)} represents the number \spad{m * 2 ** e}.
++ Though it is not assumed that the underlying integers are represented
++ with a binary \spadfunFrom{base}{FloatingPointSystem},
++ the code will be most efficient when this is the
++ the case (this is true in most implementations of Lisp).
++ The decision to choose the \spadfunFrom{base}{FloatingPointSystem} to be
++ binary has some unfortunate
++ consequences.  First, decimal numbers like 0.3 cannot be represented
++ exactly.  Second, there is a further loss of accuracy during
++ conversion to decimal for output.  To compensate for this, if d
++ digits of precision are specified, \spad{1 + ceiling(log2 d)} bits are used.
++ Two numbers that are displayed identically may therefore be
++ not equal.  On the other hand, a significant efficiency loss would
++ be incurred if we chose to use a decimal \spadfunFrom{base}{FloatingPointSystem} when the underlying
++ integer base is binary.
++
++ Algorithms used:
++ For the elementary functions, the general approach is to apply
++ identities so that the taylor series can be used, and, so
++ that it will converge within \spad{O( sqrt n )} steps.  For example,
++ using the identity \spad{exp(x) = exp(x/2)**2}, we can compute
++ \spad{exp(1/3)} to n digits of precision as follows.  We have
++ \spad{exp(1/3) = exp(2 ** (-sqrt s) / 3) ** (2 ** sqrt s)}.
++ The taylor series will converge in less than sqrt n steps and the
++ exponentiation requires sqrt n multiplications for a total of
++ \spad{2 sqrt n} multiplications.  Assuming integer multiplication costs
++ \spad{O( n**2 )} the overall running time is \spad{O( sqrt(n) n**2 )}.
++ This approach is the best known approach for precisions up to
++ about 10,000 digits at which point the methods of Brent
++ which are \spad{O( log(n) n**2 )} become competitive.  Note also that
++ summing the terms of the taylor series for the elementary
++ functions is done using integer operations.  This avoids the
++ overhead of floating point operations and results in efficient
++ code at low precisions.  This implementation makes no attempt
++ to reuse storage, relying on the underlying system to do
++ \spadgloss{garbage collection}.  I estimate that the efficiency of this
++ package at low precisions could be improved by a factor of 2
++ if in-place operations were available.
++
++ Running times: in the following, n is the number of bits of precision
++      \spad{*}, \spad{/}, \spad{sqrt}, \spad{pi}, \spad{exp1}, \spad{log2}, \spad{log10}: \spad{ O( n**2 )}
++      \spad{exp}, \spad{log}, \spad{sin}, \spad{atan}:  \spad{ O( sqrt(n) n**2 )}
++ The other elementary functions are coded in terms of the ones above.
Float():
 Join(FloatingPointSystem, DifferentialRing, ConvertibleTo String,_
  CoercibleTo DoubleFloat, TranscendentalFunctionCategory, _
    ConvertibleTo InputForm,ConvertibleFrom SF) with
   /  : (%, I) -> %
      ++ x / i computes the division from x by an integer i.
   **: (%, %) -> %
      ++ x ** y computes \spad{exp(y log x)} where \spad{x >= 0}.
   normalize: % -> %
      ++ normalize(x) normalizes x at current precision.
   relerror : (%, %) -> I
      ++ relerror(x,y) computes the absolute value of \spad{x - y} divided by
      ++ y, when \spad{y \~= 0}.
   shift: (%, I) -> %
      ++ shift(x,n) adds n to the exponent of float x.
   rationalApproximation: (%, N) -> RN
     ++ rationalApproximation(f, n) computes a rational approximation
     ++ r to f with relative error \spad{< 10**(-n)}.
   rationalApproximation: (%, N, N) -> RN
     ++ rationalApproximation(f, n, b) computes a rational
     ++ approximation r to f with relative error \spad{< b**(-n)}, that is
     ++ \spad{|(r-f)/f| < b**(-n)}.
   log2 : () -> %
      ++ log2() returns \spad{ln 2}, i.e. \spad{0.6931471805...}.
   log10: () -> %
      ++ log10() returns \spad{ln 10}: \spad{2.3025809299...}.
   exp1 : () -> %
      ++ exp1() returns  exp 1: \spad{2.7182818284...}.
   atan : (%,%) -> %
      ++ atan(x,y) computes the arc tangent from x with phase y.
   log2 : % -> %
      ++ log2(x) computes the logarithm for x to base 2.
   log10: % -> %
      ++ log10(x) computes the logarithm for x to base 10.
   outputFloating: () -> Void
      ++ outputFloating() sets the output mode to floating (scientific) notation, i.e.
      ++ \spad{mantissa * 10 exponent} is displayed as  \spad{0.mantissa E exponent}.
   outputFloating: N -> Void
      ++ outputFloating(n) sets the output mode to floating (scientific) notation
      ++ with n significant digits displayed after the decimal point.
   outputFixed: () -> Void
      ++ outputFixed() sets the output mode to fixed point notation;
      ++ the output will contain a decimal point.
   outputFixed: N -> Void
      ++ outputFixed(n) sets the output mode to fixed point notation,
      ++ with n digits displayed after the decimal point.
   outputGeneral: () -> Void
      ++ outputGeneral() sets the output mode (default mode) to general
      ++ notation; numbers will be displayed in either fixed or floating
      ++ (scientific) notation depending on the magnitude.
   outputGeneral: N -> Void
      ++ outputGeneral(n) sets the output mode to general notation
      ++ with n significant digits displayed.
   outputSpacing: N -> Void
      ++ outputSpacing(n) inserts a space after n (default 10) digits on output;
      ++ outputSpacing(0) means no spaces are inserted.
   arbitraryPrecision
   arbitraryExponent
  == add
   BASE ==> 2
   BITS:Reference(PI) := ref 68 -- 20 digits
   LENGTH ==> %ilength$Foreign(Builtin)
   ISQRT ==> approxSqrt$IntegerRoots(I)
   Rep := Record( mantissa:I, exponent:I )
   StoredConstant ==> Record( precision:PI, value:% )
   UCA ==> Record( unit:%, coef:%, associate:% )
   inc ==> increasePrecision
   dec ==> decreasePrecision

   -- local utility operations
   shift2 : (I,I) -> I           -- WSP: fix bug in shift
   times : (%,%) -> %            -- multiply x and y with no rounding
   itimes: (I,%) -> %            -- multiply by a small integer
   chop: (%,PI) -> %             -- chop x at p bits of precision
   dvide: (%,%) -> %             -- divide x by y with no rounding
   square: (%,I) -> %            -- repeated squaring with chopping
   power: (%,I) -> %             -- x ** n with chopping
   plus: (%,%) -> %              -- addition with no rounding
   sub: (%,%) -> %               -- subtraction with no rounding
   negate: % -> %                -- negation with no rounding
   ceillog10base2: PI -> PI      -- rational approximation
   floorln2: PI -> PI            -- rational approximation

   atanSeries: % -> %            -- atan(x) by taylor series |x| < 1/2
   atanInverse: I -> %           -- atan(1/n) for n an integer > 1
   expInverse: I -> %            -- exp(1/n) for n an integer
   expSeries: % -> %             -- exp(x) by taylor series  |x| < 1/2
   logSeries: % -> %             -- log(x) by taylor series 1/2 < x < 2
   sinSeries: % -> %             -- sin(x) by taylor series |x| < 1/2
   cosSeries: % -> %             -- cos(x) by taylor series |x| < 1/2
   piRamanujan: () -> %          -- pi using Ramanujans series

   shift2(x,y) == sign(x)*shift(sign(x)*x,y)

   asin x ==
      zero? x => 0
      negative? x => -asin(-x)
      one? x => pi()/2
      x > 1 => error "asin: argument > 1 in magnitude"
      inc 5; r := atan(x/sqrt(sub(1,times(x,x)))); dec 5
      normalize r

   acos x ==
      zero? x => pi()/2
      negative? x => (inc 3; r := pi()-acos(-x); dec 3; normalize r)
      one? x => 0
      x > 1 => error "acos: argument > 1 in magnitude"
      inc 5; r := atan(sqrt(sub(1,times(x,x)))/x); dec 5
      normalize r

   atan(x,y) ==
      x = 0 =>
         positive? y => pi()/2
         negative? y => -pi()/2
         0
      -- Only count on first quadrant being on principal branch.
      theta := atan abs(y/x)
      if negative? x then theta := pi() - theta
      if negative? y then theta := - theta
      theta

   atan x ==
      zero? x => 0
      negative? x => -atan(-x)
      if x > 1 then
         inc 4
         r := if zero? fractionPart x and x < [bits(),0] then atanInverse wholePart x
                 else atan(1/x)
         r := pi()/2 - r
         dec 4
         return normalize r
      -- make |x| < O( 2**(-sqrt p) ) < 1/2 to speed series convergence
      -- by using the formula  atan(x) = 2*atan(x/(1+sqrt(1+x**2)))
      k := ISQRT (bits()-100)::I quo 5
      k := max(0,2 + k + order x)
      inc(2*k)
      for i in 1..k repeat x := x/(1+sqrt(1+x*x))
      t := atanSeries x
      dec(2*k)
      t := shift(t,k)
      normalize t

   atanSeries x ==
      -- atan(x) = x (1 - x**2/3 + x**4/5 - x**6/7 + ...)  |x| < 1
      p := bits() + LENGTH bits() + 2
      s:I := d:I := shift(1,p)
      y := times(x,x)
      t := m := - shift2(y.mantissa,y.exponent+p)
      for i in 3.. by 2 while t ~= 0 repeat
         s := s + t quo i
         t := (m * t) quo d
      x * [s,-p]

   atanInverse n ==
      -- compute atan(1/n) for an integer n > 1
      -- atan n = 1/n - 1/n**3/3 + 1/n**5/4 - ...
      --   pi = 16 atan(1/5) - 4 atan(1/239)
      n2 := -n*n
      e:I := bits() + LENGTH bits() + LENGTH n + 1
      s:I := shift(1,e) quo n
      t:I := s quo n2
      for k in 3.. by 2 while t ~= 0 repeat
         s := s + t quo k
         t := t quo n2
      normalize [s,-e]

   sin x ==
      s := sign x; x := abs x; p := bits(); inc 4
      if x > [6,0] then (inc p; x := 2*pi()*fractionPart(x/pi()/2); bits p)
      if x > [3,0] then (inc p; s := -s; x := x - pi(); bits p)
      if x > [3,-1] then (inc p; x := pi() - x; dec p)
      -- make |x| < O( 2**(-sqrt p) ) < 1/2 to speed series convergence
      -- by using the formula  sin(3*x/3) = 3 sin(x/3) - 4 sin(x/3)**3
      -- the running time is O( sqrt p M(p) ) assuming |x| < 1
      k := ISQRT (bits()-100)::I quo 4
      k := max(0,2 + k + order x)
      if positive? k then (inc k; x := x / 3**k::N)
      r := sinSeries x
      for i in 1..k repeat r := itimes(3,r)-shift(r**3,2)
      bits p
      s * r

   sinSeries x ==
      -- sin(x) = x (1 - x**2/3! + x**4/5! - x**6/7! + ... |x| < 1/2
      p := bits() + LENGTH bits() + 2
      y := times(x,x)
      s:I := d:I := shift(1,p)
      m:I := - shift2(y.mantissa,y.exponent+p)
      t:I := m quo 6
      for i in 4.. by 2 while t ~= 0 repeat
         s := s + t
         t := (m * t) quo (i*(i+1))
         t := t quo d
      x * [s,-p]

   cos x ==
     s:I := 1; x := abs x; p := bits(); inc 4
     if x > [6,0] then (inc p; x := 2*pi()*fractionPart(x/pi()/2); dec p)
     if x > [3,0] then (inc p; s := -s; x := x-pi(); dec p)
     if x > [1,0] then
         -- take care of the accuracy problem near pi/2
         inc p; x := pi()/2-x; bits p; x := normalize x
         return (s * sin x)
      -- make |x| < O( 2**(-sqrt p) ) < 1/2 to speed series convergence
      -- by using the formula  cos(2*x/2) = 2 cos(x/2)**2 - 1
      -- the running time is O( sqrt p M(p) ) assuming |x| < 1
     k := ISQRT (bits()-100)::I quo 3
     k := max(0,2 + k + order x)
      -- need to increase precision by more than k, otherwise recursion
      -- causes loss of accuracy.
      -- Michael Monagan suggests adding a factor of log(k)
     if positive? k then (inc(k+length(k)**2); x := shift(x,-k))
     r := cosSeries x
     for i in 1..k repeat r := shift(r*r,1)-1
     bits p
     s * r



   cosSeries x ==
      -- cos(x) = 1 - x**2/2! + x**4/4! - x**6/6! + ... |x| < 1/2
      p := bits() + LENGTH bits() + 1
      y := times(x,x)
      s:I := d:I := shift(1,p)
      m:I := - shift2(y.mantissa,y.exponent+p)
      t:I := m quo 2
      for i in 3.. by 2 while t ~= 0 repeat
         s := s + t
         t := (m * t) quo (i*(i+1))
         t := t quo d
      normalize [s,-p]

   tan x ==
      s := sign x; x := abs x; p := bits(); inc 6
      if x > [3,0] then (inc p; x := pi()*fractionPart(x/pi()); dec p)
      if x > [3,-1] then (inc p; x := pi()-x; s := -s; dec p)
      if x > 1 then (c := cos x; t := sqrt(1-c*c)/c)
      else (c := sin x; t := c/sqrt(1-c*c))
      bits p
      s * t

   P:StoredConstant := [1,[1,2]]
   pi() ==
      -- We use Ramanujan's identity to compute pi.
      -- The running time is quadratic in the precision.
      -- This is about twice as fast as Machin's identity on Lisp/VM
      --   pi = 16 atan(1/5) - 4 atan(1/239)
      bits() <= P.precision => normalize P.value
      (P := [bits(), piRamanujan()]) value

   piRamanujan() ==
      -- Ramanujans identity for 1/pi
      -- Reference: Shanks and Wrench, Math Comp, 1962
      -- "Calculation of pi to 100,000 Decimals".
      n := bits() + LENGTH bits() + 11
      t:I := shift(1,n) quo 882
      d:I := 4*882**2
      s:I := 0
      for i in 2.. by 2 for j in 1123.. by 21460 while t ~= 0 repeat
         s := s + j*t
         m := -(i-1)*(2*i-1)*(2*i-3)
         t := (m*t) quo (d*i**3)
      1 / [s,-n-2]

   sinh x ==
      zero? x => 0
      lost:I := max(- order x,0)
      2*lost > bits() => x
      inc(5+lost); e := exp x; s := (e-1/e)/2; dec(5+lost)
      normalize s

   cosh x ==
      (inc 5; e := exp x; c := (e+1/e)/2; dec 5; normalize c)

   tanh x ==
      zero? x => 0
      lost:I := max(- order x,0)
      2*lost > bits() => x
      inc(6+lost); e := exp x; e := e*e; t := (e-1)/(e+1); dec(6+lost)
      normalize t

   asinh x ==
      p := min(0,order x)
      if zero? x or 2*p < -bits() then return x
      inc(5-p); r := log(x+sqrt(1+x*x)); dec(5-p)
      normalize r

   acosh x ==
      if x < 1 then error "invalid argument to acosh"
      inc 5; r := log(x+sqrt(sub(times(x,x),1))); dec 5
      normalize r

   atanh x ==
      if x > 1 or x < -1 then error "invalid argument to atanh"
      p := min(0,order x)
      if zero? x or 2*p < -bits() then return x
      inc(5-p); r := log((x+1)/(1-x))/2; dec(5-p)
      normalize r

   log x ==
      negative? x => error "negative log"
      zero? x => error "log 0 generated"
      p := bits(); inc 5
      -- apply  log(x) = n log 2 + log(x/2**n)  so that  1/2 < x < 2
      if negative?(n := order x) then n := n+1
      l: % :=
        n = 0 => 0
        x := shift(x,-n)
        n * log2()
      -- speed the series convergence by finding m and k such that
      -- | exp(m/2**k) x - 1 |  <  1 / 2 ** O(sqrt p)
      -- write  log(exp(m/2**k) x) as m/2**k + log x
      k := ISQRT (p-100)::I quo 3
      if k > 1 then
         k := max(1,k+order(x-1))
         inc k
         ek := expInverse (2**k::N)
         dec(p quo 2); m := order square(x,k); inc(p quo 2)
         m := (6847196937 * m) quo 9878417065   -- m := m log 2
         x := x * ek ** (-m)
         l := l + [m,-k]
      l := l + logSeries x
      bits p
      normalize l

   logSeries x ==
      -- log(x) = 2 y (1 + y**2/3 + y**4/5 ...)  for  y = (x-1) / (x+1)
      -- given 1/2 < x < 2 on input we have -1/3 < y < 1/3
      p := bits() + (g := LENGTH bits() + 3)
      inc g; y := (x-1)/(x+1); dec g
      s:I := d:I := shift(1,p)
      z := times(y,y)
      t := m := shift2(z.mantissa,z.exponent+p)
      for i in 3.. by 2 while t ~= 0 repeat
         s := s + t quo i
         t := m * t quo d
      y * [s,1-p]

   L2:StoredConstant := [1,1]
   log2() ==
      --  log x  =  2 * sum( ((x-1)/(x+1))**(2*k+1)/(2*k+1), k=1.. )
      --  log 2  =  2 * sum( 1/9**k / (2*k+1), k=0..n ) / 3
      n := bits() :: N
      n <= L2.precision => normalize L2.value
      n := n + LENGTH n + 3  -- guard bits
      s:I := shift(1,n+1) quo 3
      t:I := s quo 9
      for k in 3.. by 2 while t ~= 0 repeat
         s := s + t quo k
         t := t quo 9
      L2 := [bits(),[s,-n]]
      normalize L2.value

   L10:StoredConstant := [1,[1,1]]
   log10() ==
      --  log x  =  2 * sum( ((x-1)/(x+1))**(2*k+1)/(2*k+1), k=0.. )
      --  log 5/4  =  2 * sum( 1/81**k / (2*k+1), k=0.. ) / 9
      n := bits() :: N
      n <= L10.precision => normalize L10.value
      n := n + LENGTH n + 5  -- guard bits
      s:I := shift(1,n+1) quo 9
      t:I := s quo 81
      for k in 3.. by 2 while t ~= 0 repeat
         s := s + t quo k
         t := t quo 81
      -- We have log 10 = log 5 + log 2 and log 5/4 = log 5 - 2 log 2
      inc 2; L10 := [bits(),[s,-n] + 3*log2()]; dec 2
      normalize L10.value

   log2(x) == (inc 2; r := log(x)/log2(); dec 2; normalize r)
   log10(x) == (inc 2; r := log(x)/log10(); dec 2; normalize r)

   exp(x) ==
      -- exp(n+x) = exp(1)**n exp(x) for n such that |x| < 1
      p := bits(); inc 5; e1:% := 1
      if (n := wholePart x) ~= 0 then
         inc LENGTH n; e1 := exp1() ** n; dec LENGTH n
         x := fractionPart x
      if zero? x then (bits p; return normalize e1)
      -- make |x| < O( 2**(-sqrt p) ) < 1/2 to speed series convergence
      -- by repeated use of the formula exp(2*x/2) = exp(x/2)**2
      -- results in an overall running time of O( sqrt p M(p) )
      k := ISQRT (p-100)::I quo 3
      k := max(0,2 + k + order x)
      if positive? k then (inc k; x := shift(x,-k))
      e := expSeries x
      if positive? k then e := square(e,k)
      bits p
      e * e1

   expSeries x ==
      -- exp(x) = 1 + x + x**2/2 + ... + x**i/i!  valid for all x
      p := bits() + LENGTH bits() + 1
      s:I := d:I := shift(1,p)
      t:I := n:I := shift2(x.mantissa,x.exponent+p)
      for i in 2.. while t ~= 0 repeat
         s := s + t
         t := (n * t) quo i
         t := t quo d
      normalize [s,-p]

   expInverse k ==
      -- computes exp(1/k) via continued fraction
      p0:I := 2*k+1; p1:I := 6*k*p0+1
      q0:I := 2*k-1; q1:I := 6*k*q0+1
      for i in 10*k.. by 4*k while 2 * LENGTH p0 < bits() repeat
         (p0,p1) := (p1,i*p1+p0)
         (q0,q1) := (q1,i*q1+q0)
      dvide([p1,0],[q1,0])

   E:StoredConstant := [1,[1,1]]
   exp1() ==
      if bits() > E.precision then E := [bits(),expInverse 1]
      normalize E.value

   sqrt x ==
      negative? x => error "negative sqrt"
      m := x.mantissa; e := x.exponent
      l := LENGTH m
      p := 2 * bits() - l + 2
      if odd?(e-l) then p := p - 1
      i := shift2(x.mantissa,p)
      -- ISQRT uses a variable precision newton iteration
      i := ISQRT i
      normalize [i,(e-p) quo 2]

   bits() == deref BITS
   bits(n) == (t := bits(); setref(BITS,n); t)
   precision() == bits()
   precision(n) == bits(n)
   increasePrecision n == (b := bits(); bits((b + n)::PI); b)
   decreasePrecision n == (b := bits(); bits((b - n)::PI); b)
   ceillog10base2 n == ((13301 * n + 4003) quo 4004) :: PI
   digits() == max(1,4004 * (bits()-1) quo 13301)::PI
   digits(n) == (t := digits(); bits (1 + ceillog10base2 n); t)

   order(a) == LENGTH a.mantissa + a.exponent - 1
   relerror(a,b) == order((a-b)/b)
   0 == [0,0]
   1 == [1,0]
   base() == BASE
   mantissa x == x.mantissa
   exponent x == x.exponent
   one? a == a = 1
   zero? a == zero?(a.mantissa)
   negative? a == negative?(a.mantissa)
   positive? a == positive?(a.mantissa)

   chop(x,p) ==
      e : I := LENGTH x.mantissa - p
      if positive? e then x := [shift2(x.mantissa,-e),x.exponent+e]
      x
   float(m,e) == normalize [m,e]
   float(m,e,b) ==
      m = 0 => 0
      inc 4; r := m * [b,0] ** e; dec 4
      normalize r
   normalize x ==
      m := x.mantissa
      m = 0 => 0
      e : I := LENGTH m - bits()
      if positive? e then
         y := shift2(m,1-e)
         if odd? y then
            y := (if positive? y then y+1 else y-1) quo 2
            if LENGTH y > bits() then
               y := y quo 2
               e := e+1
         else y := y quo 2
         x := [y,x.exponent+e]
      x
   shift(x:%,n:I) == [x.mantissa,x.exponent+n]

   x = y ==
      order x = order y and sign x = sign y and zero? (x - y)
   x < y ==
      y.mantissa = 0 => negative? x.mantissa
      x.mantissa = 0 => positive? y.mantissa
      negative? x and positive? y => true
      negative? y and positive? x => false
      order x < order y => positive? x
      order x > order y => negative? x
      negative? (x-y)

   abs x == if negative? x then -x else normalize x
   ceiling x ==
      if negative? x then return (-floor(-x))
      if zero? fractionPart x then x else truncate x + 1
   wholePart x == shift2(x.mantissa,x.exponent)
   floor x == if negative? x then -ceiling(-x) else truncate x
   round x == (half := [sign x,-1]; truncate(x + half))
   sign x == if negative? x.mantissa then -1 else 1
   truncate x ==
      if x.exponent >= 0 then return x
      normalize [shift2(x.mantissa,x.exponent),0]
   recip(x) == if x=0 then "failed" else 1/x
   differentiate x == 0

   - x == normalize negate x
   negate x == [-x.mantissa,x.exponent]
   x + y == normalize plus(x,y)
   x - y == normalize plus(x,negate y)
   sub(x,y) == plus(x,negate y)
   plus(x,y) ==
      mx := x.mantissa; my := y.mantissa
      mx = 0 => y
      my = 0 => x
      ex := x.exponent; ey := y.exponent
      ex = ey => [mx+my,ex]
      de := ex + LENGTH mx - ey - LENGTH my
      de > bits()+1 => x
      de < -(bits()+1) => y
      if ex < ey then (mx,my,ex,ey) := (my,mx,ey,ex)
      mw := my + shift2(mx,ex-ey)
      [mw,ey]

   x:% * y:% == normalize times (x,y)
   x:I * y:% ==
      if LENGTH x > bits() then normalize [x,0] * y
      else normalize [x * y.mantissa,y.exponent]
   x:% / y:% == normalize dvide(x,y)
   x:% / y:I ==
      if LENGTH y > bits() then x / normalize [y,0] else x / [y,0]
   inv x == 1 / x

   times(x:%,y:%) == [x.mantissa * y.mantissa, x.exponent + y.exponent]
   itimes(n:I,y:%) == [n * y.mantissa,y.exponent]

   dvide(x,y) ==
      ew := LENGTH y.mantissa - LENGTH x.mantissa + bits() + 1
      mw := shift2(x.mantissa,ew) quo y.mantissa
      ew := x.exponent - y.exponent - ew
      [mw,ew]

   square(x,n) ==
      ma := x.mantissa; ex := x.exponent
      for k in 1..n repeat
         ma := ma * ma; ex := ex + ex
         l:I := bits()::I - LENGTH ma
         ma := shift2(ma,l); ex := ex - l
      [ma,ex]

   power(x,n) ==
      y:% := 1; z:% := x
      repeat
         if odd? n then y := chop( times(y,z), bits() )
         if (n := n quo 2) = 0 then return y
         z := chop( times(z,z), bits() )

   x:% ** y:% ==
      x = 0 =>
         y = 0 => error "0**0 is undefined"
         negative? y => error "division by 0"
         0
      y = 0 => 1
      y = 1 => x
      x = 1 => 1
      p := abs order y + 5
      inc p; r := exp(y*log(x)); dec p
      normalize r

   x:% ** r:RN ==
      x = 0 =>
         r = 0 => error "0**0 is undefined"
         negative? r => error "division by 0"
         0
      r = 0 => 1
      r = 1 => x
      x = 1 => 1
      n := numer r
      d := denom r
      negative? x =>
         odd? d =>
            odd? n => return -((-x)**r)
            return ((-x)**r)
         error "negative root"
      if d = 2 then
         inc LENGTH n; y := sqrt(x); y := y**n; dec LENGTH n
         return normalize y
      y := [n,0]/[d,0]
      x ** y

   x:% ** n:I ==
      x = 0 =>
         n = 0 => error "0**0 is undefined"
         negative? n => error "division by 0"
         0
      n = 0 => 1
      n = 1 => x
      x = 1 => 1
      p := bits()
      bits(p + LENGTH n + 2)
      y := power(x,abs n)
      if negative? n then y := dvide(1,y)
      bits p
      normalize y

   -- Utility routines for conversion to decimal
   ceilLength10: I -> I
   chop10: (%,I) -> %
   convert10:(%,I) -> %
   floorLength10: I -> I
   length10: I -> I
   normalize10: (%,I) -> %
   quotient10: (%,%,I) -> %
   power10: (%,I,I) -> %
   times10: (%,%,I) -> %

   convert10(x,d) ==
      m := x.mantissa; e := x.exponent
      --!! deal with bits here
      b := bits(); (q,r) := divide(abs e, b)
      b := 2**b::N; r := 2**r::N
      -- compute 2**e = b**q * r
      h := power10([b,0],q,d+5)
      h := chop10([r*h.mantissa,h.exponent],d+5)
      if negative? e then h := quotient10([m,0],h,d)
      else times10([m,0],h,d)

   ceilLength10 n == 146 * LENGTH n quo 485 + 1
   floorLength10 n == 643 *  LENGTH n quo 2136
   length10 n ==
      ln := LENGTH(n:=abs n)
      upper := 76573 * ln quo 254370
      lower := 21306 * (ln-1) quo 70777
      upper = lower => upper + 1
      n := n quo (10**lower::N)
      while n >= 10 repeat
         n:= n quo 10
         lower := lower + 1
      lower + 1

   chop10(x,p) ==
      e : I := floorLength10 x.mantissa - p
      if positive? e then x := [x.mantissa quo 10**e::N,x.exponent+e]
      x
   normalize10(x,p) ==
      ma := x.mantissa
      ex := x.exponent
      e : I := length10 ma - p
      if positive? e then
         ma := ma quo 10**(e-1)::N
         ex := ex + e
         (ma,r) := divide(ma, 10)
         if r > 4 then
            ma := ma + 1
            if ma = 10**p::N then (ma := 1; ex := ex + p)
      [ma,ex]
   times10(x,y,p) == normalize10(times(x,y),p)
   quotient10(x,y,p) ==
      ew := floorLength10 y.mantissa - ceilLength10 x.mantissa + p + 2
      if negative? ew then ew := 0
      mw := (x.mantissa * 10**ew::N) quo y.mantissa
      ew := x.exponent - y.exponent - ew
      normalize10([mw,ew],p)
   power10(x,n,d) ==
      x = 0 => 0
      n = 0 => 1
      n = 1 => x
      x = 1 => 1
      p:I := d + LENGTH n + 1
      e:I := n
      y:% := 1
      z:% := x
      repeat
         if odd? e then y := chop10(times(y,z),p)
         if (e := e quo 2) = 0 then return y
         z := chop10(times(z,z),p)

   --------------------------------
   -- Output routines for Floats --
   --------------------------------
   zero ==> char("0")
   separator ==> space()$Character

   SPACING : Reference(N) := ref 10
   OUTMODE : Reference(S) := ref "general"
   OUTPREC : Reference(I) := ref(-1)

   fixed : % -> S
   floating : % -> S
   general : % -> S

   padFromLeft(s:S):S ==
      zero? deref SPACING => s
      n:I := #s - 1
      t := new( (n + 1 + n quo deref SPACING) :: N , separator )
      for i in 0..n for j in minIndex t .. repeat
         t.j := s.(i + minIndex s)
         if (i+1) rem deref SPACING = 0 then j := j+1
      t
   padFromRight(s:S):S ==
      deref SPACING = 0 => s
      n:I := #s - 1
      t := new( (n + 1 + n quo deref SPACING) :: N , separator )
      for i in n..0 by -1 for j in maxIndex t .. by -1 repeat
         t.j := s.(i + minIndex s)
         if (n-i+1) rem deref SPACING = 0 then j := j-1
      t

   fixed f ==
      zero? f => "0.0"
      zero? exponent f =>
        padFromRight concat(string mantissa f, ".0")
      negative? f => concat("-", fixed abs f)
      d := if deref OUTPREC = -1 then digits()::I else deref OUTPREC
--    g := convert10(abs f,digits); m := g.mantissa; e := g.exponent
      g := convert10(abs f,d); m := g.mantissa; e := g.exponent
      if deref OUTPREC ~= -1 then
         -- round g to OUTPREC digits after the decimal point
         l := length10 m
         if -e > deref OUTPREC and -e < 2*digits()::I then
            g := normalize10(g,l+e+deref OUTPREC)
            m := g.mantissa; e := g.exponent
      s := string m; n := #s; o := e+n
      p := if deref OUTPREC = -1 then n::I else deref OUTPREC
      t:S
      if e >= 0 then
         s := concat(s, new(e::N, zero))
         t := ""
      else if o <= 0 then
         t := concat(new((-o)::N,zero), s)
         s := "0"
      else
         t := s(o + minIndex s .. n + minIndex s - 1)
         s := s(minIndex s .. o + minIndex s - 1)
      n := #t
      if deref OUTPREC = -1 then
         t := rightTrim(t,zero)
         if t = "" then t := "0"
      else if n > p then t := t(minIndex t .. p + minIndex t- 1)
                    else t := concat(t, new((p-n)::N,zero))
      concat(padFromRight s, concat(".", padFromLeft t))

   floating f ==
      zero? f => "0.0"
      negative? f => concat("-", floating abs f)
      t:S := if zero? deref SPACING then "E" else " E "
      zero? exponent f =>
        s := string mantissa f
        concat ["0.", padFromLeft s, t, string(#s)]
      -- base conversion to decimal rounded to the requested precision
      d := if deref OUTPREC = -1 then digits()::I else deref OUTPREC
      g := convert10(f,d); m := g.mantissa; e := g.exponent
      -- I'm assuming that length10 m = # s given n > 0
      s := string m; n := #s; o := e+n
      s := padFromLeft s
      concat ["0.", s, t, string o]

   general(f) ==
      zero? f => "0.0"
      negative? f => concat("-", general abs f)
      d := if deref OUTPREC = -1 then digits()::I else deref OUTPREC
      zero? exponent f =>
        d := d + 1
        s := string mantissa f
        deref OUTPREC ~= -1 and (e := #s) > d =>
          t:S := if zero? deref SPACING then "E" else " E "
          concat ["0.", padFromLeft s, t, string e]
        padFromRight concat(s, ".0")
      -- base conversion to decimal rounded to the requested precision
      g := convert10(f,d); m := g.mantissa; e := g.exponent
      -- I'm assuming that length10 m = # s given n > 0
      s := string m; n := #s; o := n + e
      -- Note: at least one digit is displayed after the decimal point
      -- and trailing zeroes after the decimal point are dropped
      if positive? o and o <= max(n,d) then
         -- fixed format: add trailing zeroes before the decimal point
         if o > n then s := concat(s, new((o-n)::N,zero))
         t := rightTrim(s(o + minIndex s .. n + minIndex s - 1), zero)
         if t = "" then t := "0" else t := padFromLeft t
         s := padFromRight s(minIndex s .. o + minIndex s - 1)
         concat(s, concat(".", t))
      else if o <= 0 and o >= -5 then
         -- fixed format: up to 5 leading zeroes after the decimal point
         concat("0.",padFromLeft concat(new((-o)::N,zero),rightTrim(s,zero)))
      else
         -- print using E format written  0.mantissa E exponent
         t := padFromLeft rightTrim(s,zero)
         s := if zero? deref SPACING then "E" else " E "
         concat ["0.", t, s, string(e+n)]

   outputSpacing n == setref(SPACING,n)
   outputFixed() == (setref(OUTMODE,"fixed"); setref(OUTPREC, -1))
   outputFixed n == (setref(OUTMODE,"fixed"); setref(OUTPREC,n::I))
   outputGeneral() == (setref(OUTMODE,"general"); setref(OUTPREC,-1))
   outputGeneral n == (setref(OUTMODE,"general"); setref(OUTPREC,n::I))
   outputFloating() == (setref(OUTMODE,"floating"); setref(OUTPREC,-1))
   outputFloating n == (setref(OUTMODE,"floating"); setref(OUTPREC,n::I))

   convert(f):S ==
      b:Integer :=
        deref OUTPREC = -1 and not zero? f =>
          bits(length(abs mantissa f)::PositiveInteger)
        0
      s :=
        deref OUTMODE = "fixed" => fixed f
        deref OUTMODE = "floating" => floating f
        deref OUTMODE = "general" => general f
        empty()$String
      if positive? b then bits(b::PositiveInteger)
      s = empty()$String => error "bad output mode"
      s

   coerce(f):OutputForm ==
     f >= 0 => message(convert(f)@S)
     - (coerce(-f)@OutputForm)

   convert(f):InputForm ==
     convert [convert('float), convert mantissa f,
              convert exponent f, convert base()]$List(InputForm)

   -- Conversion routines
   convert(x:%):Float == x pretend Float
   convert(x:%):SF == makeSF(x.mantissa,x.exponent)$Lisp
   coerce(x:%):SF == convert(x)@SF
   convert(sf:SF):% == float(mantissa sf,exponent sf,base()$SF)

   retract(f:%):RN == rationalApproximation(f,(bits()-1)::N,BASE)

   retractIfCan(f:%):Union(RN, "failed") ==
     rationalApproximation(f,(bits()-1)::N,BASE)

   retract(f:%):I ==
     (f = (n := wholePart f)::%) => n
     error "Not an integer"

   retractIfCan(f:%):Union(I, "failed") ==
     (f = (n := wholePart f)::%) => n
     "failed"

   rationalApproximation(f,d) == rationalApproximation(f,d,10)

   rationalApproximation(f,d,b) ==
      t: Integer
      nu := f.mantissa; ex := f.exponent
      if ex >= 0 then return ((nu*BASE**(ex::N))/1)
      de := BASE**((-ex)::N)
      if b < 2 then error "base must be > 1"
      tol := b**d
      s := nu; t := de
      p0,p1,q0,q1 : Integer
      p0 := 0; p1 := 1; q0 := 1; q1 := 0
      repeat
         (q,r) := divide(s, t)
         p2 := q*p1+p0
         q2 := q*q1+q0
         if r = 0 or tol*abs(nu*q2-de*p2) < de*abs(p2) then return (p2/q2)
         (p0,p1) := (p1,p2)
         (q0,q1) := (q1,q2)
         (s,t) := (t,r)

