Rat and rationalize have quite different definitions:
rat(fl) produces the simplest rational within ratepsilon of fl
rationalize(fl) converts fl to its *exact* rational representation = M *
2^N. That is usually *not *going to be the simplest rational such that
float(rationalize(fl)) == fl.
(%i1) rationalize(1.0/3);
(%o1) 6004799503160661/18014398509481984
(%i2) factor(%);
(%o2) 3^3*7*19*73*87211*262657 / 2^54
(%i3) cf(%o1);
(%o3) [0,3,6004799503160661]
But float(1/3) == 1.0/3
On Tue, Jul 3, 2012 at 10:06 AM, Henry Baker <hbaker1 at pipeline.com> wrote:
> 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.
>
>