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.
>
>