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.