diff options
Diffstat (limited to 'src/algebra/contfrac.spad.pamphlet')
| -rw-r--r-- | src/algebra/contfrac.spad.pamphlet | 425 |
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} |