axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / FLOAT.spad

)abbrev domain FLOAT Float

++ Author: Michael Monagan
++ Date Created:  December 1987
++ Change History: 19 Jun 1990
++ References:
++ Corl00 According to Abramowitz and Stegun or arccoth needn't be Uncouth
++ Fate01a A Critique of OpenMath and Thoughts on Encoding Mathematics
++ 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,base)} for integer
++ \spad{mantissa}, \spad{exponent} specifies the number
++ \spad{mantissa * base ** 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**(-bits)}. 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 base 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 base, 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 base 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(10^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 base 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\br
++ \spad{*}, \spad{/}, \spad{sqrt}, \spad{pi}, \spad{exp1}, \spad{log2}, 
++ \spad{log10}: \spad{ O( n**2 )} \br
++ \spad{exp}, \spad{log}, \spad{sin}, \spad{atan}: \spad{O(sqrt(n) n**2)}\br
++ The other elementary functions are coded in terms of the ones above.

Float() : SIG == CODE where

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

  SIG ==> Join(FloatingPointSystem, DifferentialRing, ConvertibleTo String, 
               OpenMath,  CoercibleTo DoubleFloat, 
               TranscendentalFunctionCategory, ConvertibleTo InputForm) 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}, \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.

     convert : SF -> %
       ++ convert(x) converts a \spadtype{DoubleFloat} x to a \spadtype{Float}.

     outputFloating : () -> Void
       ++ outputFloating() sets the output mode to floating (scientific) 
       ++ notation, \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 space after n (default 10) digits on output;
       ++ outputSpacing(0) means no spaces are inserted.

     arbitraryPrecision

     arbitraryExponent

  CODE ==> add

   BASE ==> 2

   BITS:Reference(PI) := ref 68 -- 20 digits

   LENGTH ==> INTEGER_-LENGTH$Lisp

   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

   writeOMFloat(dev: OpenMathDevice, x: %): Void ==
      OMputApp(dev)
      OMputSymbol(dev, "bigfloat1", "bigfloat")
      OMputInteger(dev, mantissa x)
      OMputInteger(dev, 2)
      OMputInteger(dev, exponent x)
      OMputEndApp(dev)

   OMwrite(x: %): String ==
      s: String := ""
      sp := OM_-STRINGTOSTRINGPTR(s)$Lisp
      dev: OpenMathDevice := OMopenString(sp pretend String, OMencodingXML)
      OMputObject(dev)
      writeOMFloat(dev, x)
      OMputEndObject(dev)
      OMclose(dev)
      s := OM_-STRINGPTRTOSTRING(sp)$Lisp pretend String
      s

   OMwrite(x: %, wholeObj: Boolean): String ==
      s: String := ""
      sp := OM_-STRINGTOSTRINGPTR(s)$Lisp
      dev: OpenMathDevice := OMopenString(sp pretend String, OMencodingXML)
      if wholeObj then
         OMputObject(dev)
      writeOMFloat(dev, x)
      if wholeObj then
         OMputEndObject(dev)
      OMclose(dev)
      s := OM_-STRINGPTRTOSTRING(sp)$Lisp pretend String
      s

   OMwrite(dev: OpenMathDevice, x: %): Void ==
      OMputObject(dev)
      writeOMFloat(dev, x)
      OMputEndObject(dev)

   OMwrite(dev: OpenMathDevice, x: %, wholeObj: Boolean): Void ==
      if wholeObj then
         OMputObject(dev)
      writeOMFloat(dev, x)
      if wholeObj then
         OMputEndObject(dev)
   
   shift2(x,y) == sign(x)*shift(sign(x)*x,y)

   asin x ==
      zero? x => 0
      negative? x => -asin(-x)
      (x = 1) => 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)
      (x = 1) => 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 =>
         y > 0 => pi()/2
         y < 0 => -pi()/2
         0
      -- Only count on first quadrant being on principal branch.
      theta := atan abs(y/x)
      if x < 0 then theta := pi() - theta
      if y < 0 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 k > 0 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 k > 0 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 (n := order x) < 0 then n := n+1
      l := if n = 0 then 0 else (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 k > 0 then (inc k; x := shift(x,-k))
      e := expSeries x
      if k > 0 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() == BITS()

   bits(n) == (t := bits(); 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 e > 0 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 e > 0 then
         y := shift2(m,1-e)
         if odd? y then
            y := (if y>0 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 => x.mantissa < 0
      x.mantissa = 0 => y.mantissa > 0
      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 x.mantissa < 0 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"
         y < 0 => error "division by 0"
         y > 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"
         r < 0 => error "division by 0"
         r > 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"
         n < 0 => error "division by 0"
         n > 0 => 0
      n = 0 => 1
      n = 1 => x
      x = 1 => 1
      p := bits()
      bits(p + LENGTH n + 2)
      y := power(x,abs n)
      if n < 0 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 e < 0 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 e > 0 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 e > 0 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 ew < 0 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? SPACING() => s
      n:I := #s - 1
      t := new( (n + 1 + n quo SPACING()) :: N , separator )
      for i in 0..n for j in minIndex t .. repeat
         t.j := s.(i + minIndex s)
         if (i+1) rem SPACING() = 0 then j := j+1
      t
   padFromRight(s:S):S ==
      SPACING() = 0 => s
      n:I := #s - 1
      t := new( (n + 1 + n quo 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 SPACING() = 0 then j := j-1
      t

   fixed f ==
      d := if OUTPREC() = -1 then digits::I else OUTPREC()
      dpos:N:= if (d > 0) then d::N else 1::N
      zero? f =>
        OUTPREC() = -1 => "0.0"
        concat("0",concat(".",padFromLeft new(dpos,zero)))
      zero? exponent f =>
        concat(padFromRight convert(mantissa f)@S,
               concat(".",padFromLeft new(dpos,zero)))
      negative? f => concat("-", fixed abs f)
      bl := LENGTH(f.mantissa) + f.exponent
      dd :=
        OUTPREC() = -1 => d
        bl > 0 => (146*bl) quo 485 + 1 + d
        d
      g := convert10(abs f,dd)
      m := g.mantissa
      e := g.exponent
      if OUTPREC() ^= -1 then
         -- round g to OUTPREC digits after the decimal point
         l := length10 m
         if -e > OUTPREC() and -e < 2*digits::I then
            g := normalize10(g,l+e+OUTPREC())
            m := g.mantissa; e := g.exponent
      s := convert(m)@S; n := #s; o := e+n
      p := if OUTPREC() = -1 then n::I else 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 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? SPACING() then "E" else " E "
      zero? exponent f =>
        s := convert(mantissa f)@S
        concat ["0.", padFromLeft s, t, convert(#s)@S]
      -- base conversion to decimal rounded to the requested precision
      d := if OUTPREC() = -1 then digits::I else OUTPREC()
      g := convert10(f,d); m := g.mantissa; e := g.exponent
      -- I'm assuming that length10 m = # s given n > 0
      s := convert(m)@S; n := #s; o := e+n
      s := padFromLeft s
      concat ["0.", s, t, convert(o)@S]

   general(f) ==
      zero? f => "0.0"
      negative? f => concat("-", general abs f)
      d := if OUTPREC() = -1 then digits::I else OUTPREC()
      zero? exponent f =>
        d := d + 1
        s := convert(mantissa f)@S
        OUTPREC() ^= -1 and (e := #s) > d =>
          t:S := if zero? SPACING() then "E" else " E "
          concat ["0.", padFromLeft s, t, convert(e)@S]
        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 := convert(m)@S; 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 o > 0 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? SPACING() then "E" else " E "
         concat ["0.", t, s, convert(e+n)@S]

   outputSpacing n == SPACING() := n

   outputFixed() == (OUTMODE() := "fixed"; OUTPREC() := -1)

   outputFixed n == (OUTMODE() := "fixed"; OUTPREC() := n::I)

   outputGeneral() == (OUTMODE() := "general"; OUTPREC() := -1)

   outputGeneral n == (OUTMODE() := "general"; OUTPREC() := n::I)

   outputFloating() == (OUTMODE() := "floating"; OUTPREC() := -1)

   outputFloating n == (OUTMODE() := "floating"; OUTPREC() := n::I)

   convert(f):S ==
      b:Integer :=
        OUTPREC() = -1 and not zero? f =>
          bits(length(abs mantissa f)::PositiveInteger)
        0
      s :=
        OUTMODE() = "fixed" => fixed f
        OUTMODE() = "floating" => floating f
        OUTMODE() = "general" => general f
        empty()$String
      if b > 0 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"::Symbol), 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)