Not a problem, exactly. I'm curious as to how "rat" is smart enough to
stop the continued fraction process before it starts giving bogus
results, whereas "rationalize" isn't smart enough to stop while it's
ahead.
BTW, cf(sqrt(2)) should be the infinite repeating continued fraction
[1,2,2,2,2,2,...]
Maxima 5.25.1 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) sqrt(2),numer;
(%o1) 1.414213562373095
(%i2) ratepsilon;
(%o2) 2.0E-8
(%i3) cf(rat(%o1));
rat: replaced 1.414213562373095 by 8119/5741 = 1.414213551646055
(%o3) [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
(%i4) ratepsilon:1.0e-15;
(%o4) 1.0000000000000001E-15
(%i5) cf(rat(%o1));
rat: replaced 1.414213562373095 by 22619537/15994428 = 1.414213562373097
(%o5) [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
(%i6) ratepsilon:1.0e-16;
(%o6) 9.9999999999999998E-17
(%i7) cf(rat(%o1));
rat: replaced 1.414213562373095 by 131836323/93222358 = 1.414213562373095
(%o7) [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
(%i8) ratepsilon:1.0e-17;
(%o8) 1.0000000000000001E-17
(%i9) cf(rat(%o1));
rat: replaced 1.414213562373095 by 131836323/93222358 = 1.414213562373095
(%o9) [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
(%i10) rationalize(%o1);
6369051672525773
(%o10) ----------------
4503599627370496
(%i11) cf(%);
(%o11) [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1,
1, 2, 7, 1, 2, 33, 2, 7, 5, 2, 1, 1, 16, 2]
(%i12) ratepsilon:1.0e-20;
(%o12) 9.9999999999999995E-21
(%i13) cf(rat(%o1));
rat: replaced 1.414213562373095 by 131836323/93222358 = 1.414213562373095
(%o13) [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
(%i14)
At 06:30 AM 7/3/2012, Stavros Macrakis wrote:
>Not sure I understand what the issue is. Can you give an example?
>
>On Mon, Jul 2, 2012 at 2:54 PM, Henry Baker <hbaker1 at pipeline.com> wrote:
>Over the weekend, I was playing with my toy version of complex-rationalize.
>
>The problem in programming this function isn't the complex numbers themselves,
>but what the "stopping criterion" for the continued fraction should be.
>
>In particular, I notice that the Lisp (& probably the Maxima) version of
>"rationalize" doesn't do a very good job of stopping. For example,
>the CF expansion of sqrt(2)=[1,2,2,2,2,2,...], so if the CF expansion
>starts spitting out non-2's, then Houston, we have a Problem. The
>Common Lisp implementations I have (GCL, SBCL) both seem to go too far.
>
>According to the Common Lisp manual, "rationalize" may return any rational
>number for which the floating-point number is the best available approximation
>of its format; in doing this it attempts to keep both numerator an denominator
>small.
>
>The Common Lisp manual further goes on to stipulate that
>
>(float (rationalize x) x) == x
>
>So, one might conclude that the problem isn't rationalize itself, but the
>implementation of sqrt(2).
>
>Nevertheless, I think that "rationalize" should have a second argument
>to indicate the precision of the result -- e.g., |x-x'|/|x|<epsilon,
>where x is the number being rationalized, and x' is the rational
>approximation.
>
>But even implementing this stopping criterion is fraught with problems.
>
>It is possible to implement a recursive Euclidean algorithm which
>is _tail-recursive_, i.e., it builds the product of 2x2 matrices
>while recursing down, rather than while returning back up (this
>works due to the associativity of matrix multiplication). When
>the stopping criterion is reached, the matrix (and hence the
>rational approximation) has already been built and merely needs
>to be returned.
>
>But even when this is done, the programming task is hampered by
>the fact that Common Lisp rationals won't allow "1/0" (one/zero).
>This is a pity, because the 1) IEEE floating point numbers already
>allow such an unsigned infinity; and 2) the most perspicuous code
>allows rationals of the form "1/0" to exist, because error'ing
>out forces too many, and too early, tests to avoid them.