On 2/16/07, Andreas Eder <aeder at arcor.de> wrote:
>
> > The code for continued fractions for square roots appears to
> > compute everything with integer arithmetic. It would be ideal if
> > other roots could be handled so nicely; I don't know how to go
> > about that. Be that as it may, cf should handle special cases
> > for which it knows particular algorithms, and fall back on something
> > like cf(bfloat(expr)) in the general case.
>
> Quadratic irrationals have a rather simple cf expansion - it is
> essentially periodic and vice versa. But as far sa I know no athere
> simple scheme for cf expansions is known.
The bfloat method would actually be an improvement in all but the pure
rational and sqrt cases. Currently, cf gets *incorrect results* for
algebraic combinations of sqrts and integers, see below. As for the general
integer^rational (with the new code), it can be extremely slow (see below).
Using bfloat is much more straightforward. Here's some simple code:
cfb(ex,terms) :=
if ratnump(ex) then cf(ex) else
block([res,fpprec:3],
while length(res:cf(bfloat(ex))) < terms+3 do fpprec:fpprec*2,
/* correct rounding */
res: rest( reverse(res), length(res)-terms-1 ),
reverse( if res[1]=1 then cons(res[2]+1,rest(res,2)) else
rest(res,1) ) )$
This handles *arbitrary* expressions, including transcendental functions.
There is one important case where it fails, namely when the result is in
fact rational, e.g. for sin(2)^2+cos(2)^2 or asin(1/sqrt(2))/%pi. I'm close
to done with the code for these cases (which warns the user that the answer
is probably correct, but uses approximation methods), but it is a bit fiddly
and I don't want to distribute it yet. Simply giving up after N iterations
would be easy but a bit crude...
In the meantime, I'd like to know if sees other flaws....
-s
Examples where current code is wrong (and also produces lots of spurious
extra terms):
cf(sqrt(2/3)),cflength:12 =>
[0, 1, 4, 2, 4, 2, 4, 2, 4, 2, 3, 1, 281, 4, 2, 1, 9, 1, 3, 2, 2, 14]
^^^ first incorrect term @ position 11
(would be correct if rounded to 11 terms, but not to 12)
cf((3/5)^(1/3)),cflength:5 =>
[0, 1, 5, 3, 2, 2, 2]
^^^ first incorrect term @ position 4
cf(100000^(1/3)),cflength:5 == cflength(10*100^(1/3)) =>
[46,2,2,2,8,1,1 ...]
^ first incorrect term @ 3
Example where current code is very slow:
cf(21^(9/10)),cflength:3 =>
[15,2,20] -- correct, but takes 423 seconds for 3 terms
cfb(21^(9/10),100) takes 0.34 seconds for 100 terms