[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Axiom-developer] [ContinuedFractions] (new)
From: |
Bill Page |
Subject: |
[Axiom-developer] [ContinuedFractions] (new) |
Date: |
Mon, 30 Jan 2006 19:44:55 -0600 |
Changes http://wiki.axiom-developer.org/ContinuedFractions/diff
--
Description:
'ContinuedFraction' implements general continued fractions.
This version is not restricted to simple, finite fractions and
uses the '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 'R' (i.e. 'sizeLess?(0, x)').
\begin{spad}
)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:
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
\end{spad}
--
forwarded from http://wiki.axiom-developer.org/address@hidden