Subject: [clisp-list] floating point problems with clisp 2.30

From: "Vadim V. Zhytnikov" <[email-masked]>
Date: Sun Oct 6 01:34:01 2002

Hi!

It seems that there is some floating point bug in
recently released clisp 2.30.  I've been trying to
compile Maxima CVS with clisp 2.30 and suddenly
I've got large number of

WARNING:
Floating point operation combines numbers of different precision.
See ANSI CL 12.1.4.4 and the CLISP impnotes for details.
The result's actual precision is controlled by
*FLOATING-POINT-CONTAGION-ANSI*.
To shut off this warning, set *WARN-ON-FLOATING-POINT-CONTAGION* to NIL.

warnings in several Maxima modules floating point
evaluation of elliptic and Bessel functions.
Of course I can suppress these warnings by
setting *WARN-ON-FLOATING-POINT-CONTAGION* to NIL
but I'm perfectly sure that they are spurious
and caused by some clisp error because

1. Although Maxima can be compiled it now gives
wrong value for Bessel and elliptic functions
for both values of *FLOATING-POINT-CONTAGION-ANSI*
(T or NIL). Result is not totally meaningless but
in both cases we get right digits only for single
precision while double is expected.

2. With clisp 2.29 everything is just fine.

I tried to apply recent patch from CVS
to lisparit.d but the problem still
persist.

Best wishes,

-- 
      Vadim V. Zhytnikov

       <[email-masked]>
      <[email-masked]>
      <[email-masked]>
     <[email-masked]>

---

From: Sam Steingold <[email-masked]>
Date: Sun Oct 6 06:33:35 2002

> * In message <[email-masked]>
> * On the subject of "[clisp-list] floating point problems with clisp 2.30"
> * Sent on Sun, 06 Oct 2002 12:33:26 +0300
> * Honorable "Vadim V. Zhytnikov" <[email-masked]> writes:

> It seems that there is some floating point bug in recently released
> clisp 2.30.  I've been trying to compile Maxima CVS with clisp 2.30
> and suddenly I've got large number of
> WARNING:

if you are not "combining numbers of different precision", this is a bug.
please isolate it to specific funcalls (by setting *BREAK-ON-SIGNALS*
to T and looking at the backtrace) and report it to <clisp-list>.

Please investigate every single case when you get the warning.

Thanks.

-- 
Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux
<http://www.camera.org> <http://www.iris.org.il> <http://www.memri.org/>
<http://www.mideasttruth.com/> <http://www.palestine-central.com/links.html>
There are two ways to write error-free programs; only the third one works.

---

From: "Vadim V. Zhytnikov" <[email-masked]>
Date: Sun Oct 6 11:55:02 2002

>Please investigate every single case when you get the warning.
Here is one piece of the problem which explains
some excessive warnings

[1]> (expt 5.0d0 0.3d0)

WARNING:
1.6206565966927624d0
[2]>

Obviously warning is not required here.

Unfortunately this doesn't explain wrong numeric values
which I get for some special functions in Maxima with clisp 2.30
since the value above is correct.

Best wishes,

-- 
      Vadim V. Zhytnikov

       <[email-masked]>
      <[email-masked]>
      <[email-masked]>
     <[email-masked]>

---

From: Sam Steingold <[email-masked]>
Date: Sun Oct 6 15:35:02 2002

> * In message <[email-masked]>
> * On the subject of "Re: [clisp-list] floating point problems with clisp 2.30"
> * Sent on Sun, 06 Oct 2002 22:53:50 +0300

thanks.
this is fixed by the appended patch.

> Unfortunately this doesn't explain wrong numeric values which I get
> for some special functions in Maxima with clisp 2.30 since the value
> above is correct.

please send _all_ cases when you get the above warning.

-- 
Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux
<http://www.camera.org> <http://www.iris.org.il> <http://www.memri.org/>
<http://www.mideasttruth.com/> <http://www.palestine-central.com/links.html>
Yeah, yeah, I love cats too... wanna trade recipes?

--- realelem.d.~1.17.~	Thu Sep 12 04:14:11 2002
+++ realelem.d	Sun Oct  6 18:16:53 2002
@@ -250,7 +250,7 @@
                    return F_to_SF(x),
                    return F_to_FF(x),
                    return F_to_DF(x),
-                   return F_to_LF(x,I_to_UL(O(LF_digits))),
+                   return F_to_LF(x,lfloat_length(TheLfloat(y))),
                    pushSTACK(x), x = popSTACK());
 }
 local object RA_R_float_F (object x, object y)
@@ -259,7 +259,7 @@
                    return RA_to_SF(x),
                    return RA_to_FF(x),
                    return RA_to_DF(x),
-                   return RA_to_LF(x,I_to_UL(O(LF_digits))),
+                   return RA_to_LF(x,lfloat_length(TheLfloat(y))),
                    pushSTACK(x), x = popSTACK());
 }
 local object R_R_float_F (object x, object y)

---

From: Sam Steingold <[email-masked]>
Date: Mon Oct 7 08:19:07 2002

> for some special functions in Maxima with clisp 2.30 since the value

note that BESSEL &co call explicit (FLOAT * 1.0) which can easily
result in floating point contagion.

again, _please_ supply me with an test example where a CLISP built-in
function produces incorrect results (either WRT precision, of value or
whatever).

thanks.

-- 
Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux
<http://www.camera.org> <http://www.iris.org.il> <http://www.memri.org/>
<http://www.mideasttruth.com/> <http://www.palestine-central.com/links.html>
Murphy's Law was probably named after the wrong guy.

---

From: Raymond Toy <[email-masked]>
Date: Mon Oct 7 09:21:07 2002

>>>>> "Sam" == Sam Steingold <[email-masked]> writes:

    >> for some special functions in Maxima with clisp 2.30 since the value

    Sam> note that BESSEL &co call explicit (FLOAT * 1.0) which can easily
    Sam> result in floating point contagion.

Can you point out where that is?  I can't find any such things in the
current sources.

Also, be aware that maxima defines (float x) to be (lisp::float x 1d0)
somwhere.  Very confusing until I realized that.

Ray

---

From: Sam Steingold <[email-masked]>
Date: Mon Oct 7 09:38:05 2002

> * In message <[email-masked]>
> * Sent on 07 Oct 2002 12:20:32 -0400
> * Honorable Raymond Toy <[email-masked]> writes:
>

I meant (FLOAT N) and (FLOAT (F* 2 N)) in (DEFUN BESSEL) in bessel.lisp

> somwhere.  Very confusing until I realized that.

yuk.
this goes together with (re)defining TYPE-ERROR, LOOP-FINISH, IF &c.
someone's gotta do a cleanup.

at any rate, I would greatly appreciate some test cases
(based on today's CVS).

-- 
Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux
<http://www.camera.org> <http://www.iris.org.il> <http://www.memri.org/>
<http://www.mideasttruth.com/> <http://www.palestine-central.com/links.html>
main(a){printf(a,34,a="main(a){printf(a,34,a=%c%s%c,34);}",34);}

---

From: Raymond Toy <[email-masked]>
Date: Mon Oct 7 09:47:05 2002

>> Can you point out where that is?  I can't find any such things in the

    Sam> I meant (FLOAT N) and (FLOAT (F* 2 N)) in (DEFUN BESSEL) in bessel.lisp

I think these are ok because maxima redefines float. (Well, not
really.  This float is defined in the maxima package, not the lisp
package.)

    >> somwhere.  Very confusing until I realized that.

    Sam> yuk.
    Sam> this goes together with (re)defining TYPE-ERROR, LOOP-FINISH, IF &c.
    Sam> someone's gotta do a cleanup.

I think that's slated for a later release.  There are lots of silly
issues like this such as macros redefined as functions and vice versa,
wrong number of args in various places (maybe), and so on.

Ray

---

From: Raymond Toy <[email-masked]>
Date: Mon Oct 7 09:49:05 2002

Oh, what version are you looking at?  I can't find (defun bessel ...)
in bessel.lisp in my copy of the CVS sources.

Ray

---

From: Sam Steingold <[email-masked]>
Date: Mon Oct 7 09:57:02 2002

> * In message <[email-masked]>
> * Sent on 07 Oct 2002 12:48:47 -0400

RCS file: /home/cvs/maxima/maxima/src/bessel.lisp,v
Working file: bessel.lisp
head: 1.2
line 477
Root: :pserver:[email-masked]:/home/cvs/maxima

-- 
Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux
<http://www.camera.org> <http://www.iris.org.il> <http://www.memri.org/>
<http://www.mideasttruth.com/> <http://www.palestine-central.com/links.html>
Your mouse has moved - WinNT has to be restarted for this to take effect.

---

From: Sam Steingold <[email-masked]>
Date: Mon Oct 7 09:58:07 2002

> * In message <[email-masked]>
> * Sent on 07 Oct 2002 12:46:40 -0400

you will go a long way if you look at the CLISP warnings and eliminate
them.  Or give me CVS write access and I will do that.

-- 
Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux
<http://www.camera.org> <http://www.iris.org.il> <http://www.memri.org/>
<http://www.mideasttruth.com/> <http://www.palestine-central.com/links.html>
Never succeed from the first try - if you do, nobody will think it was hard.

---

From: Raymond Toy <[email-masked]>
Date: Mon Oct 7 10:09:08 2002

>> in bessel.lisp in my copy of the CVS sources.

    Sam> RCS file: /home/cvs/maxima/maxima/src/bessel.lisp,v
    Sam> Working file: bessel.lisp
    Sam> head: 1.2
    Sam> line 477
    Sam> Root: :pserver:[email-masked]:/home/cvs/maxima

Ok, that explains it.  The maxima code has moved to
maxima.sourceforge.net.  A lot of the issues you've mentioned,
however, still exist there.

I think the fixes you're interested in are meant for something after
the imminent 5.9.0 release.  5.9.0 is just an "update" to 5.6 so that
clisp, gcl, and cmucl can build maxima.

Ray

---

From: "Vadim V. Zhytnikov" <[email-masked]>
Date: Mon Oct 7 13:17:07 2002

>function produces incorrect results (either WRT precision, of value or
Unfortunately I can't provide any simple example yet but below are results
which make me think that something is wrong with clisp 2.30 or at least
something is wrong in combination Maxima + clisp 2.30.
All results are obtained on different lisps with Maxima CVS (SourceForge):
-------------------------------------------------
CMUCL 18d
(C1) gamma(1.1);

(D1)                     .9513507698668733
-------------------------------------------------
GCL 2.5.0 CVS

(C1) gamma(1.1);

(D1)                    0.95135076986687
-------------------------------------------------
CLISP 2.29
(C1) gamma(1.1);
(D1)                     .9513507698668733
-------------------------------------------------
CLISP 2.30
(C1) gamma(1.1);

WARNING:
<snip>
(D1)                    0.9513507
(C2) :lisp (setq custom::*floating-point-contagion-ansi* t)
T
(C2) gamma(1.1);

WARNING:
<snip>
(D2)                          .9513507355077885
------------------------------------------------------

Notice that all lisps including clisp 2.29 (and 2.28, and 2.27)
give exactly the same answer and only clisp 2.30
result differs.   Actual computation is performed by
gammafloat and gamma-lanczos functions in csimp2.lisp.

I'm glad that Raymond pointed out that Maxima redefines
float.  It solves some mystery to me.  Originally I extracted
gamma-lanczos into separate file to play with and obtained
quite strange results which differ from ones obtained from
Maxima.  Now I understand why.

-- 
     Vadim V. Zhytnikov

      <[email-masked]> 
     <[email-masked]>
     <[email-masked]>
    <[email-masked]>

---

From: Sam Steingold <[email-masked]>
Date: Mon Oct 7 14:55:05 2002

> * In message <[email-masked]>
> * Sent on Tue, 08 Oct 2002 00:15:18 +0300
> * Honorable "Vadim V. Zhytnikov" <[email-masked]> writes:
>
> Sam Steingold write:
> 
> >function produces incorrect results (either WRT precision, of value
> >or whatever).

> Unfortunately I can't provide any simple example yet

try (setq *break-on-warnings* t)
please do try to make a stand-alone (non-maxima) case!

I just fixed some (but not all) problems.

Live Lisp and prosper.

Index: src/comptran.d
===================================================================
RCS file: /cvsroot/clisp/clisp/src/comptran.d,v
retrieving revision 1.13
retrieving revision 1.16
diff -u -w -b -r1.13 -r1.16
--- src/comptran.d	11 Sep 2002 23:15:23 -0000	1.13
+++ src/comptran.d	7 Oct 2002 21:47:37 -0000	1.16
@@ -31,16 +31,27 @@
     return R_exp_R(x,start_p,end_p);
   } else { /* x=a+bi */
     pushSTACK(TheComplex(x)->c_real); /* save a */
-    R_cos_sin_R_R(TheComplex(x)->c_imag,start_p,NULL); /* (cos b), (sin b) */
-    /* stack layout: a, cos(b), sin(b). */
+    pushSTACK(TheComplex(x)->c_imag); /* save b */
+    pushSTACK(R_R_contagion_R(STACK_0,STACK_1));
+    /* since x is complex, the result is a float anyway */
+    if (R_rationalp(STACK_1)) /* b */
+      STACK_1 = RA_R_float_F(STACK_1,STACK_0);
+    if (R_rationalp(STACK_2)) /* a */
+      STACK_2 = RA_R_float_F(STACK_2,STACK_0);
+    if (start_p) {
+      STACK_1 = F_extend2_F(STACK_1); /* b */
+      STACK_2 = F_extend2_F(STACK_2); /* a */
+    }
+    R_cos_sin_R_R(STACK_1,false,NULL); /* (cos b), (sin b) */
+    /* stack layout: a, b, contagion, cos(b), sin(b). */
     /* b != Fixnum_0 ==> sin(b) ~= Fixnum_0. */
-    STACK_2 = R_exp_R(STACK_2,start_p,NULL); /* (exp a) */
-    /* stack layout: exp(a), cos(b), sin(b). */
+    STACK_2 = R_exp_R(STACK_4,false,NULL); /* (exp a) */
+    /* stack layout: a, exp(a), cos(b), sin(b). */
     STACK_0 = R_R_mal_R(STACK_2,STACK_0); /* (* (exp a) (sin b)) != Fixnum_0 */
     STACK_1 = R_R_mal_R(STACK_2,STACK_1); /* (* (exp a) (cos b)) */
     x = R_R_complex_C(F_R_float_F(STACK_1,*end_p),
                       F_R_float_F(STACK_0,*end_p)); /* (complex ... ...) */
-    skipSTACK(3); return x;
+    skipSTACK(5); return x;
   }
 }
 
@@ -166,7 +177,7 @@
         pushSTACK(a); pushSTACK(b);
         STACK_1 = N_log_N(STACK_1,true,&STACK_1); /* (log a) */
         STACK_0 = N_log_N(STACK_0,true,&STACK_0); /* (log b) */
-        a = N_N_durch_N(a,b); /* divide */
+        a = N_N_durch_N(STACK_1,STACK_0); /* divide */
         skipSTACK(2); return a;
       }
     }
@@ -502,10 +513,11 @@
     STACK_3 = R_R_contagion_R(STACK_3,STACK_2);
     STACK_2 = x = R_exp_R(R_minus_R(STACK_2),false,NULL); /* (exp (- b)) */
     /* stack layout: exp(-b), cos(a), sin(a). */
-    STACK_0 = R_R_mal_R(x,STACK_0); /* (* (exp (- b)) (sin a)) */
+    STACK_0 = R_R_mal_R(x,STACK_0); /* (* (exp (- b)) (sin a)), maybe 0 */
     x = R_R_mal_R(STACK_2,STACK_1); /* (* (exp (- b)) (cos a)) */
     x = R_R_complex_N(F_R_float_F(x,STACK_3), /* (complex ... ...) */
-                      F_R_float_F(STACK_0,STACK_3));
+                      eq(Fixnum_0,STACK_0) ? Fixnum_0
+                      : F_R_float_F(STACK_0,STACK_3));
     skipSTACK(4); return x;
   }
 }

Index: src/lisparit.d
===================================================================
RCS file: /cvsroot/clisp/clisp/src/lisparit.d,v
retrieving revision 1.41
retrieving revision 1.43
diff -u -w -b -r1.41 -r1.43
--- src/lisparit.d	24 Sep 2002 17:55:51 -0000	1.41
+++ src/lisparit.d	7 Oct 2002 21:13:29 -0000	1.43
@@ -980,8 +980,12 @@
 # (EXP number), CLTL S. 203
 LISPFUNN(exp,1) {
   check_number(STACK_0);
-  VALUES1(N_exp_N(STACK_0,true,&STACK_0));
-  skipSTACK(1);
+  if (complexp(STACK_0))
+    pushSTACK(R_R_contagion_R(TheComplex(STACK_0)->c_real,
+                              TheComplex(STACK_0)->c_imag));
+  else pushSTACK(STACK_0);
+  VALUES1(N_exp_N(STACK_1,true,&STACK_0));
+  skipSTACK(2);
 }
 
 LISPFUNN(expt,2)

---

From: Sam Steingold <[email-masked]>
Date: Thu Oct 10 07:46:01 2002

> for some special functions in Maxima with clisp 2.30 since the value

I fixed in the CVS all bugs you have reported so far.
please get the CLISP CVS HEAD and check that everything is OK.
Thanks.

-- 
Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux
<http://www.camera.org> <http://www.iris.org.il> <http://www.memri.org/>
<http://www.mideasttruth.com/> <http://www.palestine-central.com/links.html>
Trespassers will be shot.  Survivors will be prosecuted.

---

From: "Vadim V. Zhytnikov" <[email-masked]>
Date: Sat Oct 12 11:11:05 2002

> please get the CLISP CVS HEAD and check that everything is OK.

Today's CLISP CVS works fine and passes all
tests with Maxima CVS.

Thank you!

-- 
      Vadim V. Zhytnikov

       <[email-masked]>
      <[email-masked]>
      <[email-masked]>
     <[email-masked]>

---

From: Sam Steingold <[email-masked]>
Date: Sat Oct 12 20:00:02 2002

> * In message <[email-masked]>
> * On the subject of "Re: floating point problems with clisp 2.30"
> * Sent on Sat, 12 Oct 2002 22:10:26 +0300

what about the case when some function was given a DOUBLE-FLOAT and
returned a SINGLE-FLOAT (or a SHORT-FLOAT)?

could you please re-check that this problem is gone too?

thanks!

-- 
Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux
<http://www.camera.org> <http://www.iris.org.il> <http://www.memri.org/>
<http://www.mideasttruth.com/> <http://www.palestine-central.com/links.html>
Isn't "Microsoft Works" an advertisement lie?

---

From: "Vadim V. Zhytnikov" <[email-masked]>
Date: Mon Oct 14 12:58:05 2002

> could you please re-check that this problem is gone too?

Yes this problem is fixed as well.
In any case with CLISP CVS 12.10.2002 I can't
reproduce any problem I observed earlier with
clisp 2.30.  Everything seems to be fine.

But I want to stress that I was not able
to reduce these problems to some simple
test outside Maxima (except wrong
warning in expt) and now they all gone.

-- 
      Vadim V. Zhytnikov

       <[email-masked]>
      <[email-masked]>
      <[email-masked]>
     <[email-masked]>