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

)abbrev domain DFLOAT DoubleFloat
++ Author: Michael Monagan
++ Date Created: January 1988
++ 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{DoubleFloat} is intended to make accessible
++ hardware floating point arithmetic in Axiom, either native double
++ precision, or IEEE. On most machines, there will be hardware support for
++ the arithmetic operations: ++ +, *, / and possibly also the
++ sqrt operation.
++ The operations exp, log, sin, cos, atan are normally coded in
++ software based on minimax polynomial/rational approximations.
++
++ Some general comments about the accuracy of the operations:
++ the operations +, *, / and sqrt are expected to be fully accurate.
++ The operations exp, log, sin, cos and atan are not expected to be
++ fully accurate.  In particular, sin and cos
++ will lose all precision for large arguments.
++
++ The Float domain provides an alternative to the DoubleFloat domain.
++ It provides an arbitrary precision model of floating point arithmetic.
++ This means that accuracy problems like those above are eliminated
++ by increasing the working precision where necessary.  \spadtype{Float}
++ provides some special functions such as erf, the error function
++ in addition to the elementary functions.  The disadvantage of Float is that
++ it is much more expensive than small floats when the latter can be used.

DoubleFloat() : SIG == CODE where

  SIG ==> Join(FloatingPointSystem, DifferentialRing, OpenMath,
               TranscendentalFunctionCategory, SpecialFunctionCategory, _
               ConvertibleTo InputForm) with

    _/ : (%, Integer) -> %
      ++ x / i computes the division from x by an integer i.

    _*_* : (%,%) -> %
      ++ x ** y returns the yth power of x (equal to \spad{exp(y log x)}).

    exp1 : () -> %
      ++ exp1() returns the natural log base \spad{2.718281828...}.

    hash : % -> Integer
      ++ hash(x) returns the hash key for x

    log2 : % -> %
      ++ log2(x) computes the logarithm with base 2 for x.

    log10 : % -> %
      ++ log10(x) computes the logarithm with base 10 for x.

    atan : (%,%) -> %
      ++ atan(x,y) computes the arc tangent from x with phase y.

    Gamma : % -> %
      ++ Gamma(x) is the Euler Gamma function.

    Beta : (%,%) -> %
      ++ Beta(x,y) is \spad{Gamma(x) * Gamma(y)/Gamma(x+y)}.

    doubleFloatFormat : String -> String
      ++doubleFloatFormat changes the output format for doublefloats 
      ++using lisp format strings

    rationalApproximation : (%, NonNegativeInteger) -> Fraction Integer
      ++rationalApproximation(f, n) computes a rational approximation
      ++r to f with relative error \spad{< 10**(-n)}.

    rationalApproximation : (%, NonNegativeInteger, NonNegativeInteger) -> _
          Fraction Integer
      ++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)}).

    machineFraction : % -> Fraction Integer
      ++machineFraction(x) returns a bit-exact fraction of the machine
      ++floating point number using the common lisp integer-decode-float 
      ++function. See Steele, ISBN 0-13-152414-3 p354
      ++This function can be used to print results which do not depend
      ++on binary-to-decimal conversions
      ++
      ++X a:DFLOAT:=-1.0/3.0
      ++X machineFraction a

    integerDecode : % -> List Integer
      ++integerDecode(x) returns the multiple values of the  common 
      ++lisp integer-decode-float function. 
      ++See Steele, ISBN 0-13-152414-3 p354. This function can be used
      ++to ensure that the results are bit-exact and do not depend on
      ++the binary-to-decimal conversions.
      ++
      ++X a:DFLOAT:=-1.0/3.0
      ++X integerDecode a

  CODE ==> add

   format: String := "~G"

   MER ==> Record(MANTISSA:Integer,EXPONENT:Integer)

   manexp: % -> MER

   doubleFloatFormat(s:String): String ==
     ss: String := format
     format := s
     ss

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

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

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

   OMwrite(dev: OpenMathDevice, x: %, wholeObj: Boolean): Void ==
     if wholeObj then
       OMputObject(dev)
     OMputFloat(dev, convert x)
     if wholeObj then
       OMputEndObject(dev)

   checkComplex(x:%):% == C_-TO_-R(x)$Lisp
   -- In AKCL we used to have to make the arguments to ASIN ACOS ACOSH ATANH
   -- complex to get the correct behaviour.
   --makeComplex(x: %):% == COMPLEX(x, 0$%)$Lisp

   machineFraction(df:%):Fraction(Integer) ==
     numer:Integer:=INTEGER_-DECODE_-FLOAT_-NUMERATOR(df)$Lisp
     denom:Integer:=INTEGER_-DECODE_-FLOAT_-DENOMINATOR(df)$Lisp
     sign:Integer:=INTEGER_-DECODE_-FLOAT_-SIGN(df)$Lisp
     sign*numer/denom

   integerDecode(df:%):List(Integer) ==
     numer:Integer:=INTEGER_-DECODE_-FLOAT_-NUMERATOR(df)$Lisp
     exp:Integer:=INTEGER_-DECODE_-FLOAT_-EXPONENT(df)$Lisp
     sign:Integer:=INTEGER_-DECODE_-FLOAT_-SIGN(df)$Lisp
     [numer,exp,sign]

   base()           == FLOAT_-RADIX(0$%)$Lisp

   mantissa x       == manexp(x).MANTISSA

   exponent x       == manexp(x).EXPONENT

   precision()      == FLOAT_-DIGITS(0$%)$Lisp

   bits()           ==
     base() = 2 => precision()
     base() = 16 => 4*precision()
     wholePart(precision()*log2(base()::%))::PositiveInteger

   max()            == MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp

   min()            == MOST_-NEGATIVE_-DOUBLE_-FLOAT$Lisp

   order(a) == precision() + exponent a - 1

   0                == FLOAT(0$Lisp,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp

   1                == FLOAT(1$Lisp,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp
   -- rational approximation to e accurate to 23 digits

   exp1()  == FLOAT(534625820200,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp / _
              FLOAT(196677847971,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp

   pi()    == FLOAT(PI$Lisp,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp

   coerce(x:%):OutputForm == 
     x >= 0 => message(FORMAT(NIL$Lisp,format,x)$Lisp @ String)
     - (message(FORMAT(NIL$Lisp,format,-x)$Lisp @ String))

   convert(x:%):InputForm == convert(x pretend DoubleFloat)$InputForm

   x < y            == DFLESSTHAN(x,y)$Lisp

   - x              == DFUNARYMINUS(x)$Lisp

   x + y            == DFADD(x,y)$Lisp

   x:% - y:%        == DFSUBTRACT(x,y)$Lisp

   x:% * y:%        == DFMULTIPLY(x,y)$Lisp

   i:Integer * x:%  == DFINTEGERMULTIPLY(i,x)$Lisp

   max(x,y)         == DFMAX(x,y)$Lisp

   min(x,y)         == DFMIN(x,y)$Lisp

   x = y            == DFEQL(x,y)$Lisp

   x:% / i:Integer  == DFINTEGERDIVIDE(x,i)$Lisp

   sqrt x           == checkComplex DFSQRT(x)$Lisp

   log10 x          == checkComplex DFLOG(x,10)$Lisp

   x:% ** i:Integer == DFINTEGEREXPT(x,i)$Lisp

   x:% ** y:%       == checkComplex DFEXPT(x,y)$Lisp

   coerce(i:Integer):% == FLOAT(i,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp

   exp x            == DFEXP(x)$Lisp

   log x            == checkComplex DFLOGE(x)$Lisp

   log2 x           == checkComplex DFLOG(x,2)$Lisp

   sin x            == DFSIN(x)$Lisp

   cos x            == DFCOS(x)$Lisp

   tan x            == DFTAN(x)$Lisp

   cot x            == COT(x)$Lisp

   sec x            == SEC(x)$Lisp

   csc x            == CSC(x)$Lisp

   asin x           == checkComplex DFASIN(x)$Lisp -- can be complex

   acos x           == checkComplex DFACOS(x)$Lisp -- can be complex

   atan x           == DFATAN(x)$Lisp

   acsc x           == checkComplex ACSC(x)$Lisp

   acot x           == ACOT(x)$Lisp

   asec x           == checkComplex ASEC(x)$Lisp

   sinh x           == SINH(x)$Lisp

   cosh x           == COSH(x)$Lisp

   tanh x           == TANH(x)$Lisp

   csch x           == CSCH(x)$Lisp

   coth x           == COTH(x)$Lisp

   sech x           == SECH(x)$Lisp

   asinh x          == DFASINH(x)$Lisp

   acosh x          == checkComplex DFACOSH(x)$Lisp -- can be complex

   atanh x          == checkComplex DFATANH(x)$Lisp -- can be complex

   acsch x          == ACSCH(x)$Lisp

   acoth x          == checkComplex ACOTH(x)$Lisp

   asech x          == checkComplex ASECH(x)$Lisp

   x:% / y:%        == DFDIVIDE(x,y)$Lisp

   negative? x      == DFMINUSP(x)$Lisp

   zero? x          == ZEROP(x)$Lisp

   hash x           == SXHASH(x)$Lisp

   recip(x)         == (zero? x => "failed"; 1 / x)

   differentiate x  == 0

   SFSFUN           ==> DoubleFloatSpecialFunctions()

   sfx              ==> x pretend DoubleFloat

   sfy              ==> y pretend DoubleFloat

   airyAi x         == airyAi(sfx)$SFSFUN pretend %

   airyBi x         == airyBi(sfx)$SFSFUN pretend %

   besselI(x,y)     == besselI(sfx,sfy)$SFSFUN pretend %

   besselJ(x,y)     == besselJ(sfx,sfy)$SFSFUN pretend %

   besselK(x,y)     == besselK(sfx,sfy)$SFSFUN pretend %

   besselY(x,y)     == besselY(sfx,sfy)$SFSFUN pretend %

   Beta(x,y)        == Beta(sfx,sfy)$SFSFUN pretend %

   digamma x        == digamma(sfx)$SFSFUN pretend %

   Gamma x          == Gamma(sfx)$SFSFUN pretend %

   polygamma(x,y)   ==
       if (n := retractIfCan(x@%)@Union(Integer, "failed")) case Integer _
          and n >= 0
       then polygamma(n::Integer::NonNegativeInteger,sfy)$SFSFUN pretend %
       else error "polygamma: first argument should be a nonnegative integer"

   wholePart x            == TRUNCATE(x)$Lisp

   float(ma,ex,b)   == ma*(b::%)**ex

   convert(x:%):DoubleFloat == x pretend DoubleFloat

   convert(x:%):Float == convert(x pretend DoubleFloat)$Float

   rationalApproximation(x, d) == rationalApproximation(x, d, 10)

   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

   retract(x:%):Fraction(Integer) ==
     rationalApproximation(x, (precision() - 1)::NonNegativeInteger, base())

   retractIfCan(x:%):Union(Fraction Integer, "failed") ==
     rationalApproximation(x, (precision() - 1)::NonNegativeInteger, base())

   retract(x:%):Integer ==
     x = ((n := wholePart x)::%) => n
     error "Not an integer"

   retractIfCan(x:%):Union(Integer, "failed") ==
     x = ((n := wholePart x)::%) => n
     "failed"

   sign(x) == retract FLOAT_-SIGN(x,1)$Lisp

   abs x   == FLOAT_-SIGN(1,x)$Lisp

   manexp(x) ==
      zero? x => [0,0]
      s := sign x; x := abs x
      if x > max()$% then return [s*mantissa(max())+1,exponent max()]
      me:Record(man:%,exp:Integer) := MANEXP(x)$Lisp 
      two53:= base()**precision()
      [s*wholePart(two53 * me.man ),me.exp-precision()]

   rationalApproximation(f,d,b) ==
      -- this algorithm expresses f as n / d where d = BASE ** k
      -- then all arithmetic operations are done over the integers
      (nu, ex) := manexp f
      BASE := base()
      ex >= 0 => (nu * BASE ** (ex::NonNegativeInteger))::Fraction(Integer)
      de :Integer := BASE**((-ex)::NonNegativeInteger)
      b < 2 => error "base must be > 1"
      tol := b**d
      s := nu; t := de
      p0:Integer := 0; p1:Integer := 1; q0:Integer := 1; q1:Integer := 0
      repeat
         (q,r) := divide(s, t)
         p2 := q*p1+p0
         q2 := q*q1+q0
         r = 0 or tol*abs(nu*q2-de*p2) < de*abs(p2) => return(p2/q2)
         (p0,p1) := (p1,p2)
         (q0,q1) := (q1,q2)
         (s,t) := (t,r)

   x:% ** r:Fraction Integer ==
      zero? x =>
         zero? r => error "0**0 is undefined"
         negative? r => error "division by 0"
         0
      zero? r or (x = 1) => 1
      (r = 1) => x
      n := numer r
      d := denom r
      negative? x =>
         odd? d =>
            odd? n => return -((-x)**r)
            return ((-x)**r)
         error "negative root"
      d = 2 => sqrt(x) ** n
      x ** (n::% / d::%)