continued fraction algorithm

Dear all,

I need help understanding an algorithm used in axiom for computing with
continued fractions. More precisely, I would like to work with continued
fractions over univariate polynomials with coefficients being rational numbers.

In the introduction to the axiom package we can read

++ 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)}).

but that would make any reduced continued fraction in UP(x, FRAC INT) illegal:
the Euclidean size of any rational number, in particular of 1, is zero in that
domain. (the Euclidean size of such a polynomial is simply its degree)

Below is the most relevant part of the code, where R is (in our case) the
domain of univariate polynomials over the rationals, and Q denotes the
quotientfield thereof.

The main routine is

genFromSequence: Stream Q -> %

which takes a stream of approximants and should return a reduced continued
fraction. In fact, this works well for continued fractions over the integers
(i.e., R=Integer, Q=Fraction Integer).

Unfortunately, in the case I'm interested, the loop in genFromSequence will not
necessarily terminate. For an example that works over the integers but fails
in the domain of univariate polynomials, consider:

(1) -> d := continuedFraction(1, repeating [1],repeating [1])

(1)
1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
1 + +---+ + +---+ + +---+ + +---+ + +---+ + +---+ + +---+ + +---+ + +---+
| 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
+
1 |
+---+ + ...
| 1
Type: ContinuedFraction Integer
(2) -> approximants d

3 5 8 13 21 34 55 89
(2) [1,2,-,-,-,--,--,--,--,--,...]
2 3 5 8 13 21 34 55
Type: Stream Fraction Integer

(3) -> 2*d

(3)
1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
3 + +---+ + +---+ + +---+ + +---+ + +---+ + +---+ + +---+ + +---+ + +---+
| 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4
+
1 |
+---+ + ...
| 4
Type: ContinuedFraction Integer

The final output is computed by multiplying every approximant by 2 and then
recreating the continued fraction, i.e., genFromSequence receives as argument

[2,4,3,10/3,16/5,13/4,...]

It then goes on to compute eucWhole0 from 2, 4, 3, 10/3, where and stops there,
since eucWhole0(3) = eucWhole0(10/3) = 3

If we switch to univariate polynomials however, eucWhole of a number is always
the number itself, so the loop never terminates.

Maybe one should not check eucWhole0(lo) = eucWhole0(hi) but rather whether
these two whole parts are associated to each other. At least, then the loop
would terminate, but I do not know how to go on.

I can imagine that it will be hard to help, but I won't give up hope... Any
hints would be greatly appreciated.

Many thanks in advance,

Martin

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

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)))




Help would be greatly appreciated!

Martin
.