As I mentioned before, when ratepsilon is smaller than machine precision,
it treats it as if it were roughly machine precision. This happens
"naturally" in the algorithm, see rat3e.lisp / ration1 (which calculates in
floating point).
And please stop saying "soon enough". rationalize is doing exactly what it
is supposed to do, returning the binary fraction corresponding to the
floating-point number.
-s
On Tue, Jul 3, 2012 at 1:52 PM, Henry Baker <hbaker1 at pipeline.com> wrote:
> Once again, I ask you, how does "rat" know to stop when "rationalize"
> doesn't stop soon enough?
>
> I set ratepsilon far beyond the ability of single precision float.
>
> 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) sqrt2:sqrt(2),numer;
> (%o1) 1.414213562373095
> (%i2) ratepsilon:1.0e-20;
> (%o2) 9.9999999999999995E-21
> (%i3) rat(sqrt2);
> rat: replaced 1.414213562373095 by 131836323/93222358 = 1.414213562373095
> 131836323
> (%o3)/R/ ---------
> 93222358
> (%i4) cf(%);
> (%o4) [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
> (%i5) rationalize(sqrt2);
> 6369051672525773
> (%o5) ----------------
> 4503599627370496
> (%i6) cf(%);
> (%o6) [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]
> (%i7)
>
> At 10:43 AM 7/3/2012, Stavros Macrakis wrote:
> >Not bogus.
> >
> >float(sqrt(2)) is not the same thing as sqrt(2). rationalize has no way
> of knowing whether 1.414213562373095 was calculated via sqrt(2.0) or
> 6369051672525773.0/4503599627370496 or 6369051672525774.0/4503599627370497
> or 1414213562373095.0 * 10^-15, or ...
> >
> >Remember, rationalize is the exact value of the floating-point number,
> which is always of the form A*2^B. It is not usually the minimal rational
> which rounds to that floating-point number.
> >
> > -s
> >
> >On Tue, Jul 3, 2012 at 12:02 PM, Henry Baker <hbaker1 at pipeline.com>
> wrote:
> >[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]
> >
> >The fact that these later numbers stop being 2, 2, 2, 2, ... means that
> rationalize
> >"went too far" and started producing bogus (i.e., not sqrt(2)) results.
> >
> >At 08:39 AM 7/3/2012, Stavros Macrakis wrote:
> >>In some sense, rat is stupider than this. If you set ratepsilon smaller
> than the float precision, it won't change the result. In theory, you might
> want the simplest rational number within (e.g.) 10^-20 of the exact number
> represented by the float, but rat epsilon = 1e-20 won't give it to you.
> >>
> >>Not sure what you're getting at with the sqrt(2) example.
> >>
> >> -s
> >>
> >>On Tue, Jul 3, 2012 at 11:01 AM, Henry Baker <hbaker1 at pipeline.com>
> wrote:
> >>But "rat" seems to be smarter than this, because setting ratepsilon to
> >>a smaller number doesn't change what "rat" returns, as in my last
> >>example, below.
> >>
> >>At 07:33 AM 7/3/2012, Stavros Macrakis wrote:
> >>>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.
>
>