BUG in continued fraction, to wit cf(sqrt(31^3));

I agree that the result is wrong, but your etiological conjecture is not
correct:

     cf(3^(3/2)) =>      5
     cf(float(3^(3/2)))  =>   [5, 5, 10, 5, 10]

Unfortunately, the cf routine only correctly calculates the continued
fractions of quadratic irrationals of the simplified form a+b*sqrt(c), for
integers a, b, c and where b<>c.  It is also not correct for e.g. sqrt(3)/5,
(sqrt(3)+1)/5, etc.

This is clearly a bug.  Thanks for reporting it.

              -s

On Wed, Oct 6, 2010 at 09:00, Bernhard Marx <b.w.marx at gmx.de> wrote:

> cf(sqrt(31^3));
> SHOULD at least give something like
> [172,
>
>  1,1,1,1,56,1,14,38,3,2,5,1,5,1,1,4,1,2,4,4,31,6,1,6,1,4,2,1,3,1,10,1,2,1,1,2,1,
>
>  68,3,8,11,2,1,1,2,2,2,7,3,1,7,3,1,2,2,2,3,1,22,4,6,34,2,1,3,2,22,1,1,2,1,9,1,2,
>
>  1,13,15,1,1,1,1,1,1,1,2,10,12,1,2,4,1,1,1,18,1,1,6,1,114,4,1,171,1,4,114,1,6,1,
>
>  1,18,1,1,1,4,2,1,12,10,2,1,1,1,1,1,1,1,15,13,1,2,1,9,1,2,1,1,22,2,3,1,2,34,6,4,
>
>  22,1,3,2,2,2,1,3,7,1,3,7,2,2,2,1,1,2,11,8,3,68,1,2,1,1,2,1,10,1,3,1,2,4,1,6,1,
>  6,31,4,4,2,1,4,1,1,5,1,5,2,3,38,14,1,56,1,1,1,1,344]
>
> and NOT what I get (even with cflength:3)
> [172,
>  1,1,1,1,56,1,14,38,3,2,5,1,5,1,1,4,1,2,5,1,2]
>
> I fear the continued fraction is done rather simplistic, having a float of
> the value and continuously subtracting the floor() and taking the inverse
> of
> the remaining float, which would eat up a certain number of bits of the
> precision of the float PER iteration ...
> I fear this is a BUG, that I stumbled across even for smaller numbers than
> 29791.
>
> If this should be the case, there is an easy way to avoid the problem.
> Let r:m0*sqrt(n)+c0 be the periodic continued fraction to be evaluated.
>  (m0 and c0 are exact rationals (m0 NOT 0), n=a^2+o, n,a,o positive
> integers
>  and o <= 2a)
> f:bfloat(sqrt(n));
>
> a0:floor(bfloat(m0*f+c0));
>        /* But now instead of altering f */
> d: m0^2*n-(a0-c0)^2;
> m1:m0/d;
> c1:(a0-c0)/d;
>        /* exact computations taking of a LIMITED number of bits from the
>           INITIAL precision of f (constant!) in EACH iteration */
>
> a1:floor(bfloat(m1*f+c1));
>
> d: m1^2*n-(a1-c1)^2;
> m2:m1/d;
> c2:(a1-c1)/d;
>        /* ... */
>        /* until [mj, cj] is equal to a pair previously computed */
>
>
>
>
>
> --
> Neu: GMX De-Mail - Einfach wie E-Mail, sicher wie ein Brief!
> Jetzt De-Mail-Adresse reservieren: http://portal.gmx.net/de/go/demail
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>