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

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