[Maxima] ratepsilon

Stavros Macrakis macrakis at alum.mit.edu
Tue Jul 3 10:39:59 CDT 2012

```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
> >
> >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)
> >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.
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20120703/f4484725/attachment.html>
```