Initial population.

1 files changed, 461 insertions, 0 deletions

diff --git a/src/algebra/fortmac.spad.pamphlet b/src/algebra/fortmac.spad.pamphlet
new file mode 100644
index 00000000..8c23f9ac
--- /dev/null
+++ b/src/algebra/fortmac.spad.pamphlet
@@ -0,0 +1,461 @@
+\documentclass{article}
+\usepackage{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
+      u pretend $
+
+    retract(u:$):Integer == u pretend Integer
+
+    retractIfCan(u:$):Union(Integer,"failed") == u pretend Integer
+
+@
+\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
+    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) \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}