realroots



On 8/25/06, Stavros Macrakis <macrakis at gmail.com> wrote:

> The problem seems to be in the float-to-rational conversion routine,
> which converts any number < 1.0e-38 to 0 (regardless of the value of
> ratepsilon).  My guess is that that magic number comes from the
> PDP-10's float range (!!!).

Looks like the problem is in PREPFLOAT (src/rat3e.lisp).
Here is a patch. If I don't hear otherwise I'll commit this patch.

--- rat3e.lisp  19 Nov 2005 18:24:11 -0000      1.11
+++ rat3e.lisp  26 Aug 2006 16:29:03 -0000
@@ -972,11 +972,10 @@
                  (cons num den))))))))

 (defun prepfloat (f)
-  (cond ((< (abs f) 1.0e-37) (setq f 0.0))) ;changed 38 to 37 --wfs
   (cond (modulus (merror "Floating point meaningless unless `modulus'
= `false'"))
        ($ratprint (mtell "~&`rat' replaced ~A by" f)))
   (setq f (maxima-rationalize f))
-  (if $ratprint (mtell " ~A//~A = ~A~%"  (car f) (cdr f)
+  (if $ratprint (mtell " ~A/~A = ~A~%"  (car f) (cdr f)
                       (fpcofrat1 (car f) (cdr f))))
   f)


(I've also taken the liberty of also changing "x//y" to "x/y" ....)

Here are some results I get (Maxima 5.9.3cvs / Clisp 2.38):

rat (1e-45);
 => 1/1000000000000000088213614053064226407018659840

rat (%pi * 1e-45), numer;
 => 1/318309886183790666016233356636961235628195840

?least\-positive\-double\-float;
 => 2.2250738585072E-308

rat (?least\-positive\-double\-float);
 => 1/44942328371557897693232629769725618340449424473557664318357520289433168951375240783177119330601884005280028469967848339414697442203604155623211857659868531094441973356216371319075554900311523529863270738021251442209537670585615720368478277635206809290837627671146574559986811484619929076208839082406056034304

HTH
Robert