\documentclass{article} \usepackage{open-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} <>= )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 negative? r 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 negative? n => error "Numerators must be greater than 0." negative? d => 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 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: Integer * 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} <>= )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} <>= --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. @ <<*>>= <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}