Initial population.

1 files changed, 425 insertions, 0 deletions

diff --git a/src/algebra/contfrac.spad.pamphlet b/src/algebra/contfrac.spad.pamphlet
new file mode 100644
index 00000000..2261d76c
--- /dev/null
+++ b/src/algebra/contfrac.spad.pamphlet
@@ -0,0 +1,425 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra contfrac.spad}
+\author{Stephen M. Watt, Clifton J. Williamson}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{domain CONTFRAC ContinuedFraction}
+<<domain CONTFRAC ContinuedFraction>>=
+)abbrev domain CONTFRAC ContinuedFraction
+++ Author: Stephen M. Watt
+++ Date Created: January 1987
+++ Change History:
+++   11 April   1990
+++    7 October 1991 -- SMW: Treat whole part specially.  Added comments.
+++ Basic Operations:
+++   (Field), (Algebra),
+++   approximants, complete, continuedFraction, convergents, denominators,
+++   extend, numerators, partialDenominators, partialNumerators,
+++   partialQuotients, reducedContinuedFraction, reducedForm, wholePart
+++ Related Constructors:
+++ Also See: Fraction
+++ AMS Classifications: 11A55 11J70 11K50 11Y65 30B70 40A15
+++ Keywords: continued fraction, convergent
+++ References:
+++ Description:  \spadtype{ContinuedFraction} implements general
+++   continued fractions.  This version is not restricted to simple,
+++   finite fractions and uses the \spadtype{Stream} as a
+++   representation.  The arithmetic functions assume that the
+++   approximants alternate below/above the convergence point.
+++   This is enforced by ensuring the partial numerators and partial
+++   denominators are greater than 0 in the Euclidean domain view of \spad{R}
+++   (i.e. \spad{sizeLess?(0, x)}). 
+
+
+ContinuedFraction(R): Exports == Implementation where
+  R :     EuclideanDomain
+  Q   ==> Fraction R
+  MT  ==> MoebiusTransform Q
+  OUT ==> OutputForm
+
+  Exports ==> Join(Algebra R,Algebra Q,Field) with
+    continuedFraction:        Q -> %
+      ++ continuedFraction(r) converts the fraction \spadvar{r} with
+      ++ components of type \spad{R} to a continued fraction over
+      ++ \spad{R}.
+
+    continuedFraction:        (R, Stream R, Stream R) -> %
+      ++ continuedFraction(b0,a,b) constructs a continued fraction in
+      ++ the following way:  if \spad{a = [a1,a2,...]} and \spad{b =
+      ++ [b1,b2,...]} then the result is the continued fraction
+      ++ \spad{b0 + a1/(b1 + a2/(b2 + ...))}.
+
+    reducedContinuedFraction: (R, Stream R) -> %
+      ++ reducedContinuedFraction(b0,b) constructs a continued
+      ++ fraction in the following way:  if \spad{b = [b1,b2,...]}
+      ++ then the result is the continued fraction \spad{b0 + 1/(b1 +
+      ++ 1/(b2 + ...))}.  That is, the result is the same as
+      ++ \spad{continuedFraction(b0,[1,1,1,...],[b1,b2,b3,...])}.
+
+    partialNumerators:   % -> Stream R
+      ++ partialNumerators(x) extracts the numerators in \spadvar{x}.
+      ++ That is, if \spad{x = continuedFraction(b0, [a1,a2,a3,...],
+      ++ [b1,b2,b3,...])}, then \spad{partialNumerators(x) =
+      ++ [a1,a2,a3,...]}.
+
+    partialDenominators: % -> Stream R
+      ++ partialDenominators(x) extracts the denominators in
+      ++ \spadvar{x}.  That is, if \spad{x = continuedFraction(b0,
+      ++ [a1,a2,a3,...], [b1,b2,b3,...])}, then
+      ++ \spad{partialDenominators(x) = [b1,b2,b3,...]}.
+
+    partialQuotients:    % -> Stream R
+      ++ partialQuotients(x) extracts the partial quotients in
+      ++ \spadvar{x}.  That is, if \spad{x = continuedFraction(b0,
+      ++ [a1,a2,a3,...], [b1,b2,b3,...])}, then
+      ++ \spad{partialQuotients(x) = [b0,b1,b2,b3,...]}.
+
+    wholePart:           % -> R
+      ++ wholePart(x) extracts the whole part of \spadvar{x}.  That
+      ++ is, if \spad{x = continuedFraction(b0, [a1,a2,a3,...],
+      ++ [b1,b2,b3,...])}, then \spad{wholePart(x) = b0}.
+
+    reducedForm:         % -> %
+      ++ reducedForm(x) puts the continued fraction \spadvar{x} in
+      ++ reduced form, i.e.  the function returns an equivalent
+      ++ continued fraction of the form
+      ++ \spad{continuedFraction(b0,[1,1,1,...],[b1,b2,b3,...])}.
+
+    approximants:        % -> Stream Q
+      ++ approximants(x) returns the stream of approximants of the
+      ++ continued fraction \spadvar{x}. If the continued fraction is
+      ++ finite, then the stream will be infinite and periodic with
+      ++ period 1.
+
+    convergents:         % -> Stream Q
+      ++ convergents(x) returns the stream of the convergents of the
+      ++ continued fraction \spadvar{x}. If the continued fraction is
+      ++ finite, then the stream will be finite.
+
+    numerators:          % -> Stream R
+      ++ numerators(x) returns the stream of numerators of the
+      ++ approximants of the continued fraction \spadvar{x}. If the
+      ++ continued fraction is finite, then the stream will be finite.
+
+    denominators:        % -> Stream R
+      ++ denominators(x) returns the stream of denominators of the
+      ++ approximants of the continued fraction \spadvar{x}. If the
+      ++ continued fraction is finite, then the stream will be finite.
+
+    extend:              (%,Integer) -> %
+      ++ extend(x,n) causes the first \spadvar{n} entries in the
+      ++ continued fraction \spadvar{x} to be computed.  Normally
+      ++ entries are only computed as needed.
+
+    complete:            % -> %
+      ++ complete(x) causes all entries in \spadvar{x} to be computed.
+      ++ Normally entries are only computed as needed.  If \spadvar{x}
+      ++ is an infinite continued fraction, a user-initiated interrupt is
+      ++ necessary to stop the computation.
+
+  Implementation ==> add
+
+ -- isOrdered  ==> R is Integer
+    isOrdered  ==> R has OrderedRing and R has multiplicativeValuation
+    canReduce? ==> isOrdered or R has additiveValuation
+
+    Rec ==> Record(num: R, den: R)
+    Str ==> Stream Rec
+    Rep :=  Record(value: Record(whole: R, fract: Str), reduced?: Boolean)
+
+    import Str
+
+    genFromSequence:     Stream Q -> %
+    genReducedForm:      (Q, Stream Q, MT)    -> Stream Rec
+    genFractionA:        (Stream R,Stream R)  -> Stream Rec
+    genFractionB:        (Stream R,Stream R)  -> Stream Rec
+    genNumDen:           (R,R, Stream Rec)    -> Stream R
+
+    genApproximants:     (R,R,R,R,Stream Rec) -> Stream Q
+    genConvergents:      (R,R,R,R,Stream Rec) -> Stream Q
+    iGenApproximants:    (R,R,R,R,Stream Rec) -> Stream Q
+    iGenConvergents:     (R,R,R,R,Stream Rec) -> Stream Q
+
+    reducedForm c == 
+        c.reduced? => c
+        explicitlyFinite? c.value.fract =>
+                      continuedFraction last complete convergents c
+        canReduce? => genFromSequence approximants c
+        error "Reduced form not defined for this continued fraction."
+
+    eucWhole(a: Q): R == numer a quo denom a
+
+    eucWhole0(a: Q): R ==
+        isOrdered =>
+            n := numer a
+            d := denom a
+            q := n quo d
+            r := n - q*d
+            if r < 0 then q := q - 1
+            q
+        eucWhole a
+
+    x = y ==
+        x := reducedForm x
+        y := reducedForm y
+
+        x.value.whole ^= y.value.whole => false
+
+        xl := x.value.fract; yl := y.value.fract
+
+        while not empty? xl and not empty? yl repeat
+            frst.xl.den ^= frst.yl.den => return false
+            xl := rst xl; yl := rst yl
+        empty? xl and empty? yl
+
+    continuedFraction q == q :: %
+
+    if isOrdered then
+        continuedFraction(wh,nums,dens) == [[wh,genFractionA(nums,dens)],false]
+
+        genFractionA(nums,dens) ==
+            empty? nums or empty? dens => empty()
+            n := frst nums
+            d := frst dens
+            n < 0 => error "Numerators must be greater than 0."
+            d < 0 => error "Denominators must be greater than 0."
+            concat([n,d]$Rec, delay genFractionA(rst nums,rst dens))
+    else
+        continuedFraction(wh,nums,dens) == [[wh,genFractionB(nums,dens)],false]
+
+        genFractionB(nums,dens) ==
+            empty? nums or empty? dens => empty()
+            n := frst nums
+            d := frst dens
+            concat([n,d]$Rec, delay genFractionB(rst nums,rst dens))
+
+    reducedContinuedFraction(wh,dens) ==
+        continuedFraction(wh, repeating [1], dens)
+
+    coerce(n:Integer):% == [[n::R,empty()], true]
+    coerce(r:R):%       == [[r,   empty()], true]
+
+    coerce(a: Q): % ==
+      wh := eucWhole0 a
+      fr := a - wh::Q
+      zero? fr => [[wh, empty()], true]
+
+      l : List Rec := empty()
+      n := numer fr
+      d := denom fr
+      while not zero? d repeat
+        qr := divide(n,d)
+        l  := concat([1,qr.quotient],l)
+        n  := d
+        d  := qr.remainder
+      [[wh, construct rest reverse_! l], true]
+
+    characteristic() == characteristic()$Q
+
+
+    genFromSequence apps ==
+        lo := first apps; apps := rst apps
+        hi := first apps; apps := rst apps
+        while eucWhole0 lo ^= eucWhole0 hi repeat
+            lo := first apps; apps := rst apps
+            hi := first apps; apps := rst apps
+        wh := eucWhole0 lo
+        [[wh, genReducedForm(wh::Q, apps, moebius(1,0,0,1))], canReduce?]
+
+    genReducedForm(wh0, apps, mt) ==
+        lo: Q := first apps - wh0; apps := rst apps
+        hi: Q := first apps - wh0; apps := rst apps
+        lo = hi and zero? eval(mt, lo) => empty()
+        mt  := recip mt
+        wlo := eucWhole eval(mt, lo)
+        whi := eucWhole eval(mt, hi)
+        while wlo ^= whi repeat
+            wlo := eucWhole eval(mt, first apps - wh0); apps := rst apps
+            whi := eucWhole eval(mt, first apps - wh0); apps := rst apps
+        concat([1,wlo], delay genReducedForm(wh0, apps, shift(mt, -wlo::Q)))
+
+    wholePart c           == c.value.whole
+    partialNumerators c   == map(#1.num, c.value.fract)$StreamFunctions2(Rec,R)
+    partialDenominators c == map(#1.den, c.value.fract)$StreamFunctions2(Rec,R)
+    partialQuotients c    == concat(c.value.whole, partialDenominators c)
+
+    approximants c ==
+      empty? c.value.fract => repeating [c.value.whole::Q]
+      genApproximants(1,0,c.value.whole,1,c.value.fract)
+    convergents c ==
+      empty? c.value.fract => concat(c.value.whole::Q, empty())
+      genConvergents (1,0,c.value.whole,1,c.value.fract)
+    numerators c ==
+      empty? c.value.fract => concat(c.value.whole, empty())
+      genNumDen(1,c.value.whole,c.value.fract)
+    denominators c ==
+      genNumDen(0,1,c.value.fract)
+
+    extend(x,n) == (extend(x.value.fract,n); x)
+    complete(x) == (complete(x.value.fract); x)
+
+    iGenApproximants(pm2,qm2,pm1,qm1,fr) == delay
+      nd := frst fr
+      pm := nd.num*pm2 + nd.den*pm1
+      qm := nd.num*qm2 + nd.den*qm1
+      genApproximants(pm1,qm1,pm,qm,rst fr)
+
+    genApproximants(pm2,qm2,pm1,qm1,fr) ==
+      empty? fr => repeating [pm1/qm1]
+      concat(pm1/qm1,iGenApproximants(pm2,qm2,pm1,qm1,fr))
+
+    iGenConvergents(pm2,qm2,pm1,qm1,fr) == delay
+      nd := frst fr
+      pm := nd.num*pm2 + nd.den*pm1
+      qm := nd.num*qm2 + nd.den*qm1
+      genConvergents(pm1,qm1,pm,qm,rst fr)
+
+    genConvergents(pm2,qm2,pm1,qm1,fr) ==
+      empty? fr => concat(pm1/qm1, empty())
+      concat(pm1/qm1,iGenConvergents(pm2,qm2,pm1,qm1,fr))
+
+    genNumDen(m2,m1,fr) ==
+      empty? fr => concat(m1,empty())
+      concat(m1,delay genNumDen(m1,m2*frst(fr).num + m1*frst(fr).den,rst fr))
+
+    gen  ==> genFromSequence
+    apx  ==> approximants
+
+    c, d: %
+    a: R
+    q: Q
+    n: Integer
+
+    0 == (0$R) :: %
+    1 == (1$R) :: %
+
+    c + d   == genFromSequence map(#1 + #2, apx c, apx d)
+    c - d   == genFromSequence map(#1 - #2, apx c, rest apx d)
+    - c     == genFromSequence map(   - #1, rest apx c)
+    c * d   == genFromSequence map(#1 * #2, apx c, apx d)
+    a * d   == genFromSequence map( a * #1, apx d)
+    q * d   == genFromSequence map( q * #1, apx d)
+    n * d   == genFromSequence map( n * #1, apx d)
+    c / d   == genFromSequence map(#1 / #2, apx c, rest apx d)
+    recip c ==(c = 0 => "failed";
+	       genFromSequence map( 1 / #1, rest apx c))
+
+    showAll?: () -> Boolean
+    showAll?() ==
+      NULL(_$streamsShowAll$Lisp)$Lisp => false
+      true
+
+    zagRec(t:Rec):OUT == zag(t.num :: OUT,t.den :: OUT)
+
+    coerce(c:%): OUT ==
+      wh := c.value.whole
+      fr := c.value.fract
+      empty? fr => wh :: OUT
+      count : NonNegativeInteger := _$streamCount$Lisp
+      l : List OUT := empty()
+      for n in 1..count while not empty? fr repeat
+        l  := concat(zagRec frst fr,l)
+        fr := rst fr
+      if showAll?() then
+        for n in (count + 1).. while explicitEntries? fr repeat
+          l  := concat(zagRec frst fr,l)
+          fr := rst fr
+      if not explicitlyEmpty? fr then l := concat("..." :: OUT,l)
+      l := reverse_! l
+      e := reduce("+",l)
+      zero? wh => e
+      (wh :: OUT) + e
+
+@
+\section{package NCNTFRAC NumericContinuedFraction}
+<<package NCNTFRAC NumericContinuedFraction>>=
+)abbrev package NCNTFRAC NumericContinuedFraction
+++ Author: Clifton J. Williamson
+++ Date Created: 12 April 1990
+++ Change History:
+++ Basic Operations: continuedFraction
+++ Related Constructors: ContinuedFraction, Float
+++ Also See: Fraction
+++ AMS Classifications: 11J70 11A55 11K50 11Y65 30B70 40A15
+++ Keywords: continued fraction
+++ References:
+++ Description: \spadtype{NumericContinuedFraction} provides functions
+++   for converting floating point numbers to continued fractions.
+
+NumericContinuedFraction(F): Exports == Implementation where
+  F :     FloatingPointSystem
+  CFC ==> ContinuedFraction Integer
+  I   ==> Integer
+  ST  ==> Stream I
+
+  Exports ==> with
+    continuedFraction: F -> CFC
+      ++ continuedFraction(f) converts the floating point number
+      ++ \spad{f} to a reduced continued fraction.
+
+  Implementation ==> add
+
+    cfc: F -> ST
+    cfc(a) == delay
+      aa := wholePart a
+      zero?(b := a - (aa :: F)) => concat(aa,empty()$ST)
+      concat(aa,cfc inv b)
+
+    continuedFraction a ==
+      aa := wholePart a
+      zero?(b := a - (aa :: F)) =>
+        reducedContinuedFraction(aa,empty()$ST) 
+      if negative? b then (aa := aa - 1; b := b + 1)
+      reducedContinuedFraction(aa,cfc inv b) 
+
+@
+\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 CONTFRAC ContinuedFraction>>
+<<package NCNTFRAC NumericContinuedFraction>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}