path: root/src/algebra/fortmac.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra fortmac.spad}
\author{Mike Dewar}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain MINT MachineInteger}
<<domain MINT MachineInteger>>=
)abbrev domain MINT MachineInteger
++ Author: Mike Dewar
++ Date Created:  December 1993
++ Date Last Updated:
++ Basic Operations:
++ Related Domains:
++ Also See: FortranExpression, FortranMachineTypeCategory, MachineFloat,
++  MachineComplex
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description: A domain which models the integer representation
++ used by machines in the AXIOM-NAG link.
MachineInteger(): Exports == Implementation where

  S ==> String

  Exports ==> Join(FortranMachineTypeCategory,IntegerNumberSystem) with
    maxint : PositiveInteger -> PositiveInteger
     ++ maxint(u) sets the maximum integer in the model to u
    maxint : () -> PositiveInteger
     ++ maxint() returns the maximum integer in the model
    coerce : Expression Integer -> Expression $
     ++ coerce(x) returns x with coefficients in the domain

  Implementation  ==> Integer add

    MAXINT : PositiveInteger := 2**32

    maxint():PositiveInteger == MAXINT

    maxint(new:PositiveInteger):PositiveInteger ==
      old := MAXINT
      MAXINT := new
      old

    coerce(u:Expression Integer):Expression($) ==
      map(coerce,u)$ExpressionFunctions2(Integer,$)

    coerce(u:Integer):$ ==
      import S
      abs(u) > MAXINT => 
        message: S := concat [convert(u)@S,"  > MAXINT(",convert(MAXINT)@S,")"]
        error message
      per u

    retract(u:$):Integer == rep u

    retractIfCan(u:$):Union(Integer,"failed") == rep u

@
\section{domain MFLOAT MachineFloat}
<<domain MFLOAT MachineFloat>>=
)abbrev domain MFLOAT MachineFloat
++ Author: Mike Dewar
++ Date Created:  December 1993
++ Date Last Updated:
++ Basic Operations:
++ Related Domains:
++ Also See: FortranExpression, FortranMachineTypeCategory, MachineInteger,
++  MachineComplex
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description: A domain which models the floating point representation
++ used by machines in the AXIOM-NAG link.
MachineFloat(): Exports == Implementation where

  PI  ==> PositiveInteger
  NNI ==> NonNegativeInteger
  F   ==> Float
  I   ==> Integer
  S   ==> String
  FI  ==> Fraction Integer
  SUP ==> SparseUnivariatePolynomial
  SF  ==> DoubleFloat
 
  Exports ==> 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
    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) \undocumented{}

  Implementation ==> 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 pretend 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 pretend I)::FI**exp
      sign   : I  := sign f
      f      : FI := abs f
      newMan : I  := wholePart f
      zero? f => [0,0]$Rep
      BB     : FI := (B pretend 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 pretend 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 pretend I)**((e1-e2) pretend NNI)+m2,e2]$Rep
      e2 > e1 + POWER +2 => 
               zero? f2 => f1
               f2
      normalise [m2*(B pretend 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

@
\section{domain MCMPLX MachineComplex}
<<domain MCMPLX MachineComplex>>=
)abbrev domain MCMPLX MachineComplex
++ Date Created:  December 1993
++ Date Last Updated:
++ Basic Operations:
++ Related Domains:
++ Also See: FortranExpression, FortranMachineTypeCategory, MachineInteger,
++  MachineFloat
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description: A domain which models the complex number representation
++ used by machines in the AXIOM-NAG link.
MachineComplex():Exports == Implementation where

  Exports ==> Join (FortranMachineTypeCategory,
                    ComplexCategory(MachineFloat)) with
    coerce : Complex Float -> $
      ++ coerce(u) transforms u into a MachineComplex
    coerce : Complex Integer -> $
      ++ coerce(u) transforms u into a MachineComplex
    coerce : Complex MachineFloat -> $
      ++ coerce(u) transforms u into a MachineComplex
    coerce : Complex MachineInteger -> $
      ++ coerce(u) transforms u into a MachineComplex
    coerce : $ -> Complex Float
      ++ coerce(u) transforms u into a COmplex Float

  Implementation ==> Complex MachineFloat add

    coerce(u:Complex Float):$ == 
      complex(real(u)::MachineFloat,imag(u)::MachineFloat)

    coerce(u:Complex Integer):$ ==
      complex(real(u)::MachineFloat,imag(u)::MachineFloat)

    coerce(u:Complex MachineInteger):$ ==
      complex(real(u)::MachineFloat,imag(u)::MachineFloat)

    coerce(u:Complex MachineFloat):$ == 
      complex(real(u),imag(u))

    coerce(u:$):Complex Float ==
      complex(real(u)::Float,imag(u)::Float)

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--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.
@
<<*>>=
<<license>>

<<domain MINT MachineInteger>>
<<domain MFLOAT MachineFloat>>
<<domain MCMPLX MachineComplex>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}