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

)abbrev domain MFLOAT MachineFloat
++ Author: Mike Dewar
++ Date Created:  December 1993
++ Description: 
++ A domain which models the floating point representation
++ used by machines in the AXIOM-NAG link.

MachineFloat() : SIG == CODE where

  PI  ==> PositiveInteger
  NNI ==> NonNegativeInteger
  F   ==> Float
  I   ==> Integer
  S   ==> String
  FI  ==> Fraction Integer
  SUP ==> SparseUnivariatePolynomial
  SF  ==> DoubleFloat
 
  SIG ==> Join(FloatingPointSystem,FortranMachineTypeCategory,Field,
            RetractableTo(Float),RetractableTo(Fraction(Integer)),
              CharacteristicZero) with

    precision : PI -> PI
      ++ precision(p) sets the number of digits in the model to p

    precision : () -> PI
      ++ precision() returns the number of digits in the model 

    base : PI -> PI
      ++ base(b) sets the base of the model to b

    base : () -> PI
      ++ base() returns the base of the model

    maximumExponent : I -> I
      ++ maximumExponent(e) sets the maximum exponent in the model to e

    maximumExponent : () -> I
      ++ maximumExponent() returns the maximum exponent in the model

    minimumExponent : I -> I
      ++ minimumExponent(e) sets the minimum exponent in the model to e

    minimumExponent : () -> I
      ++ minimumExponent() returns the minimum exponent in the model

    coerce : $ -> F
      ++ coerce(u) transforms a MachineFloat to a standard Float

    coerce : MachineInteger -> $
      ++ coerce(u) transforms a MachineInteger into a MachineFloat

    mantissa : $ -> I
      ++ mantissa(u) returns the mantissa of u

    exponent : $ -> I
      ++ exponent(u) returns the exponent of u

    changeBase : (I,I,PI) -> $
      ++ changeBase(exp,man,base) is not documented

  CODE ==> add

    import F
    import FI

    Rep := Record(mantissa:I,exponent:I)

    -- Parameters of the Floating Point Representation
    P : PI := 16      -- Precision
    B : PI := 2       -- Base
    EMIN : I := -1021 -- Minimum Exponent
    EMAX : I :=  1024 -- Maximum Exponent

    -- Useful constants
    POWER : PI := 53  -- The maximum power of B which will yield P
                      -- decimal digits.
    MMAX  : PI := B**POWER 

    -- locals
    locRound:(FI)->I
    checkExponent:($)->$
    normalise:($)->$
    newPower:(PI,PI)->Void
 
    retractIfCan(u:$):Union(FI,"failed") == 
      mantissa(u)*(B/1)**(exponent(u))

    wholePart(u:$):Integer ==
       man:I:=mantissa u
       exp:I:=exponent u
       f:=
           positive? exp => man*B**(exp pretend PI)
           zero? exp => man
           wholePart(man/B**((-exp) pretend PI))

    normalise(u:$):$ ==
      -- We want the largest possible mantissa, to ensure a canonical
      -- representation.
      exp : I  := exponent u
      man : I  := mantissa u
      BB  : I  := B @ I
      sgn : I := sign man ; man := abs man
      zero? man => [0,0]$Rep
      if man < MMAX then 
        while man < MMAX repeat
          exp := exp - 1
          man := man * BB
      if man > MMAX then
        q1:FI:= man/1
        BBF:FI:=BB/1
        while wholePart(q1) > MMAX repeat
          q1:= q1 / BBF
          exp:=exp + 1
        man := locRound(q1)  
      positive?(sgn) => checkExponent [man,exp]$Rep
      checkExponent [-man,exp]$Rep

    mantissa(u:$):I == elt(u,mantissa)$Rep

    exponent(u:$):I == elt(u,exponent)$Rep

    newPower(base:PI,prec:PI):Void ==
      power   : PI := 1
      target  : PI := 10**prec
      current : PI := base
      while (current := current*base) < target repeat power := power+1
      POWER := power
      MMAX  := B**POWER
      void()

    changeBase(exp:I,man:I,base:PI):$ ==
      newExp : I  := 0
      f      : FI := man*(base @ I)::FI**exp
      sign   : I  := sign f
      f      : FI := abs f
      newMan : I  := wholePart f
      zero? f => [0,0]$Rep
      BB     : FI := (B @ I)::FI
      if newMan < MMAX then
        while newMan < MMAX repeat
          newExp := newExp - 1
          f := f*BB
          newMan := wholePart f
      if newMan > MMAX then
        while newMan > MMAX repeat
          newExp := newExp + 1
          f := f/BB
          newMan := wholePart f
      [sign*newMan,newExp]$Rep

    checkExponent(u:$):$ ==
      exponent(u) < EMIN or exponent(u) > EMAX =>
        message :S := concat(["Exponent out of range: ",
                              convert(EMIN)@S, "..", convert(EMAX)@S])$S
        error message
      u
    
    coerce(u:$):OutputForm == 
      coerce(u::F)

    coerce(u:MachineInteger):$ ==
      checkExponent changeBase(0,retract(u)@Integer,10)

    coerce(u:$):F ==
      oldDigits : PI := digits(P)$F
      r : F := float(mantissa u,exponent u,B)$Float
      digits(oldDigits)$F
      r
    
    coerce(u:F):$ ==
      checkExponent changeBase(exponent(u)$F,mantissa(u)$F,base()$F)

    coerce(u:I):$ ==
       checkExponent changeBase(0,u,10)

    coerce(u:FI):$ == (numer u)::$/(denom u)::$

    retract(u:$):FI == 
      value : Union(FI,"failed") := retractIfCan(u)
      value case "failed" => error "Cannot retract to a Fraction Integer"
      value::FI

    retract(u:$):F == u::F

    retractIfCan(u:$):Union(F,"failed") == u::F::Union(F,"failed")

    retractIfCan(u:$):Union(I,"failed") ==
      value:FI := mantissa(u)*(B @ I)::FI**exponent(u)
      zero? fractionPart(value) => wholePart(value)::Union(I,"failed")
      "failed"::Union(I,"failed")

    retract(u:$):I ==
      result : Union(I,"failed") := retractIfCan u
      result = "failed" => error "Not an Integer"
      result::I

    precision(p: PI):PI ==
      old : PI := P
      newPower(B,p)
      P := p
      old

    precision():PI == P

    base(b:PI):PI ==
      old : PI := b
      newPower(b,P)
      B := b
      old

    base():PI == B

    maximumExponent(u:I):I ==
      old : I := EMAX
      EMAX := u
      old

    maximumExponent():I == EMAX

    minimumExponent(u:I):I ==
      old : I := EMIN
      EMIN := u
      old

    minimumExponent():I == EMIN

    0 == [0,0]$Rep

    1 == changeBase(0,1,10)

    zero?(u:$):Boolean == u=[0,0]$Rep

    f1:$
    f2:$

    locRound(x:FI):I ==
      abs(fractionPart(x)) >= 1/2 => wholePart(x)+sign(x)
      wholePart(x)

    recip f1 ==
      zero? f1 => "failed"
      normalise [ locRound(B**(2*POWER)/mantissa f1),-(exponent f1 + 2*POWER)]
    
    f1 * f2 == 
      normalise [mantissa(f1)*mantissa(f2),exponent(f1)+exponent(f2)]$Rep
    
    f1 **(p:FI) ==
      ((f1::F)**p)::%

--inline

    f1 / f2 == 
      zero? f2 => error "division by zero"
      zero? f1 => 0
      f1=f2 => 1
      normalise [locRound(mantissa(f1)*B**(2*POWER)/mantissa(f2)),
         exponent(f1)-(exponent f2 + 2*POWER)]
    
    inv(f1) == 1/f1

    f1 exquo f2 == f1/f2

    divide(f1,f2) == [ f1/f2,0]

    f1 quo f2 == f1/f2

    f1 rem f2 == 0

    u:I * f1 == 
      normalise [u*mantissa(f1),exponent(f1)]$Rep

    f1 = f2 == mantissa(f1)=mantissa(f2) and exponent(f1)=exponent(f2)

    f1 + f2 == 
      m1 : I := mantissa f1
      m2 : I := mantissa f2
      e1 : I := exponent f1
      e2 : I := exponent f2
      e1 > e2 => 
--insignificance
         e1 > e2 + POWER + 2 => 
               zero? f1 => f2
               f1 
         normalise [m1*(B @ I)**((e1-e2) pretend NNI)+m2,e2]$Rep
      e2 > e1 + POWER +2 => 
               zero? f2 => f1
               f2
      normalise [m2*(B @ I)**((e2-e1) pretend NNI)+m1,e1]$Rep

    - f1 == [- mantissa f1,exponent f1]$Rep

    f1 - f2 == f1 + (-f2)

    f1 < f2 == 
      m1 : I := mantissa f1
      m2 : I := mantissa f2
      e1 : I := exponent f1
      e2 : I := exponent f2
      sign(m1) = sign(m2) =>
        e1 < e2 => true
        e1 = e2 and m1 < m2 => true
        false
      sign(m1) = 1 => false
      sign(m1) = 0 and sign(m2) = -1 => false
      true

    characteristic():NNI == 0