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.