From: Raymond Toy <toy@rt...>  20020827 21:15:36

In 2.28 and 2.29 on Solaris I get: [1]> (expt 5.0 2.0) #C(25.000002 4.3711393E6) [2]> (expt 5.0d0 2.0d0) #C(24.999999999999996d0 6.122503540205582d15) Somewhat surprising results since the correct answers are real numbers. This also happens on the x86 version. FWIW, this breaks some things with maxima. Ray 
From: Sam Steingold <sds@gn...>  20020828 13:51:42

> * In message <4nznv82ec1.fsf@...> > * On the subject of "expt bug" > * Sent on 27 Aug 2002 17:15:26 0400 > * Honorable Raymond Toy <toy@...> writes: > > In 2.28 and 2.29 on Solaris I get: > > [1]> (expt 5.0 2.0) > #C(25.000002 4.3711393E6) > [2]> (expt 5.0d0 2.0d0) > #C(24.999999999999996d0 6.122503540205582d15) > > Somewhat surprising results since the correct answers are real > numbers. > > This also happens on the x86 version. > > FWIW, this breaks some things with maxima. please try the appended patch.  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.palestinecentral.com/links.html>; Do not tell me what to do and I will not tell you where to go.  comptran.d.~1.9.~ Wed Jul 17 13:46:22 2002 +++ comptran.d Wed Aug 28 09:49:22 2002 @@ 254,6 +254,15 @@ var object x; var object y; { + if (N_floatp(y)) { + pushSTACK(x); pushSTACK(y); /* save x & y */ + R_round_I_R(y); + if (R_zerop(STACK_0)) /* y is actually an integer */ + y = STACK_1; + else y = STACK_2; + x = STACK_3; + skipSTACK(4); + } if (N_realp(y) && R_rationalp(y)) { # y rational if (RA_integerp(y)) { 
From: Raymond Toy <toy@rt...>  20020828 14:04:49

>>>>> "Sam" == Sam Steingold <sds@...> writes: Sam> please try the appended patch. This works beautifully! Thanks! Ray P.S. Did you ever get around to fixing the problem where (tan (/ pi 2)) was quite a bit different from (/ (cos (/ pi 2)))? Just curious. You said you knew where the problem was but I don't remember any resolution. Not a big deal. 
From: Sam Steingold <sds@gn...>  20020828 15:33:03

> * In message <4n8z2r2i6h.fsf@...> > * On the subject of "Re: expt bug" > * Sent on 28 Aug 2002 10:04:38 0400 > * Honorable Raymond Toy <toy@...> writes: > > >>>>> "Sam" == Sam Steingold <sds@...> writes: > Sam> please try the appended patch. > This works beautifully! Thanks! you are welcome. > P.S. Did you ever get around to fixing the problem where > (tan (/ pi 2)) was quite a bit different from (/ (cos (/ pi 2)))? > Just curious. You said you knew where the problem was but I don't > remember any resolution. Not a big deal. The problem there is that all trigonometric functions are computed from the MacLaurent series for (sin(x)/x)^2. This means taking SQRT for SIN and TAN (and many others), which creates the problem you reported. I think that this decision (to use (sin(x)/x)^2) was made in the early 1990ies when the project started, and it persisted into the CLN. Fixing the problem you reported requires replacing (sin(x)/x)^2 with something more palatable (it is not clear what exactly). Switching to CORDIC is a very serious project. Any volunteers?  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.palestinecentral.com/links.html>; A year spent in artificial intelligence is enough to make one believe in God. 
From: Raymond Toy <toy@rt...>  20020828 16:01:42

>>>>> "Sam" == Sam Steingold <sds@...> writes: >> * In message <4n8z2r2i6h.fsf@...> >> P.S. Did you ever get around to fixing the problem where >> (tan (/ pi 2)) was quite a bit different from (/ (cos (/ pi 2)))? >> Just curious. You said you knew where the problem was but I don't >> remember any resolution. Not a big deal. Sam> The problem there is that all trigonometric functions are computed from Sam> the MacLaurent series for (sin(x)/x)^2. This means taking SQRT for SIN Sam> and TAN (and many others), which creates the problem you reported. Sam> I think that this decision (to use (sin(x)/x)^2) was made in the early Sam> 1990ies when the project started, and it persisted into the CLN. Too bad. I wonder how this decision was made. Bruno? Sam> Fixing the problem you reported requires replacing (sin(x)/x)^2 with Sam> something more palatable (it is not clear what exactly). Sam> Switching to CORDIC is a very serious project. I think I used to have some CORDIC code in Lisp, but I doubt I can find it anymore. Anyway CORDIC is pretty simple to do, but getting all the tables recomputed when the precision changes can be a pain. Ray 
From: Bruno Haible <haible@il...>  20020828 16:01:06

Raymond Toy writes: > > In 2.28 and 2.29 on Solaris I get: > > > > [1]> (expt 5.0 2.0) > > #C(25.000002 4.3711393E6) > > [2]> (expt 5.0d0 2.0d0) > > #C(24.999999999999996d0 6.122503540205582d15) > > > > Somewhat surprising results since the correct answers are real > > numbers. Raymond, you know what a floatingpoint number stands for: an interval of real numbers. 2.0 stands for the interval [2  eps, 2 + 2*eps], where eps is singlefloatnegativeepsilon. Now, (expt 5.0 1.999999) is a complex number, (expt 5.0 2.000001) is a complex number. You can therefore not expect (expt 5.0 2.0) to be a real number. But I could agree if you change the result from #C(25.000002 4.3711393E6) to a less biased #C(25.000002 0.0) > > FWIW, this breaks some things with maxima. Then fix Maxima to use an exact integer as exponent. Sam writes: > please try the appended patch. The patch is wrong for above reasons. Bruno 
From: Bruno Haible <haible@il...>  20020828 16:06:27

Sam writes: > The problem there is that all trigonometric functions are computed from > the MacLaurent series for (sin(x)/x)^2. This means taking SQRT for SIN > and TAN (and many others), which creates the problem you reported. There is nothing wrong with (sin(x)/x)^2. The real problem is that the series is computed around 0, pi, 2*pi, etc.  not expecting rounding problems around pi/2. The real fix would compute a different series around pi/2, 3*pi/2, 5*pi/2 etc. Bruno 
From: Sam Steingold <sds@gn...>  20020828 16:12:27

> * In message <15724.62520.469499.305994@...> > * On the subject of "Re: expt bug" > * Sent on Wed, 28 Aug 2002 18:03:04 +0200 (CEST) > * Honorable Bruno Haible <haible@...> writes: > > But I could agree if you change the result from > #C(25.000002 4.3711393E6) > to a less biased > #C(25.000002 0.0) > > Sam writes: > > please try the appended patch. > > The patch is wrong for above reasons. the patch optimizes the (expt 5.0 2.0) computation by computing (expt 5.0 2) instead (since they are mathematically identical).  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.palestinecentral.com/links.html>; To a Lisp hacker, XML is Sexpressions with extra cruft. 
From: Bruno Haible <haible@il...>  20020828 17:17:38

Sam writes: > > But I could agree if you change the result from > > #C(25.000002 4.3711393E6) > > to a less biased > > #C(25.000002 0.0) > > > > Sam writes: > > > please try the appended patch. > > > > The patch is wrong for above reasons. > > the patch optimizes the (expt 5.0 2.0) computation by computing > (expt 5.0 2) instead (since they are mathematically identical). No they are not identical. For the same reason that (expt 5.0 2) must return 25.0 and not the exact number 25, for the same reason the imaginary part of (expt 5.0 2.0) must be the inexact 0.0 and not the exact 0. Note that the complex canonicalization rule in clisp reads #c(x 0) = x but #c(x 0.0) /= x Bruno 
From: Raymond Toy <toy@rt...>  20020828 16:17:13

>>>>> "Bruno" == Bruno Haible <haible@...> writes: Bruno> Raymond Toy writes: >> > In 2.28 and 2.29 on Solaris I get: >> > >> > [1]> (expt 5.0 2.0) >> > #C(25.000002 4.3711393E6) >> > [2]> (expt 5.0d0 2.0d0) >> > #C(24.999999999999996d0 6.122503540205582d15) >> > >> > Somewhat surprising results since the correct answers are real >> > numbers. Bruno> Raymond, you know what a floatingpoint number stands for: an interval Bruno> of real numbers. 2.0 stands for the interval [2  eps, 2 + 2*eps], Bruno> where eps is singlefloatnegativeepsilon. We've had these discussions before and I think this is were we've always disagreed. You want to treat all FP numbers as intervals. I want to treat them as exact rationals. Pretending that every FP number is an interval seems somewhat wrong, because when I write 2.0, maybe it's really from the interval [1.5, 2.5). You can't know thatonly I know that. I think computer should compute things as accurately as possible, without trying to guess what the numbers might have meant. Bruno> Now, (expt 5.0 1.999999) is a complex number, Bruno> (expt 5.0 2.000001) is a complex number. You can therefore not expect Bruno> (expt 5.0 2.0) to be a real number. Well, of course I can, because my 2.0 is the exact rational 2.0, represented as a floatingpoint number. :) Bruno> But I could agree if you change the result from Bruno> #C(25.000002 4.3711393E6) Bruno> to a less biased Bruno> #C(25.000002 0.0) I do not disagree with your reasoning, here, however. By why, not go for broke and produce the even less biased result of #c(25.0 0.0). :) Yes, I'll see what I can do with maxima about this. Ray 
From: Raymond Toy <toy@rt...>  20020828 16:21:59

>>>>> "Bruno" == Bruno Haible <haible@...> writes: Bruno> Sam writes: >> The problem there is that all trigonometric functions are computed from >> the MacLaurent series for (sin(x)/x)^2. This means taking SQRT for SIN >> and TAN (and many others), which creates the problem you reported. Bruno> There is nothing wrong with (sin(x)/x)^2. The real problem is that Bruno> the series is computed around 0, pi, 2*pi, etc.  not expecting Bruno> rounding problems around pi/2. The real fix would compute a different Bruno> series around pi/2, 3*pi/2, 5*pi/2 etc. You actually have many series for (sin(x)/x)^2 expanded at all those different integer multiples of pi? In any case, I don't really care about the rounding problems, I was just surprised that tan(x) was not closer to cos(x) for x near pi/2. Ray 
From: Bruno Haible <haible@il...>  20020828 17:31:48

Raymond Toy writes: > You actually have many series for (sin(x)/x)^2 expanded at all those > different integer multiples of pi? Of course not. I use (mod x pi) and develop around 0. Bruno 
From: Sam Steingold <sds@gn...>  20020828 16:22:30

> * In message <15724.62832.742961.452125@...> > * On the subject of "Re: expt bug" > * Sent on Wed, 28 Aug 2002 18:08:16 +0200 (CEST) > * Honorable Bruno Haible <haible@...> writes: > > Sam writes: > > > The problem there is that all trigonometric functions are computed from > > the MacLaurent series for (sin(x)/x)^2. This means taking SQRT for SIN > > and TAN (and many others), which creates the problem you reported. > > There is nothing wrong with (sin(x)/x)^2. The real problem is that > the series is computed around 0, pi, 2*pi, etc.  not expecting > rounding problems around pi/2. The real fix would compute a different > series around pi/2, 3*pi/2, 5*pi/2 etc. sure, computing (sin(x)/x)^2 at y=pi/2x is the right decision here. are you going to do this?  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.palestinecentral.com/links.html>; Lisp is a way of life. C is a way of death. 
From: Sam Steingold <sds@gn...>  20020828 16:30:56

> * In message <4nu1lf0xhf.fsf@...> > * On the subject of "Re: expt bug" > * Sent on 28 Aug 2002 12:17:00 0400 > * Honorable Raymond Toy <toy@...> writes: > > >>>>> "Bruno" == Bruno Haible <haible@...> writes: > > Bruno> Raymond Toy writes: > >> > In 2.28 and 2.29 on Solaris I get: > >> > > >> > [1]> (expt 5.0 2.0) > >> > #C(25.000002 4.3711393E6) > >> > [2]> (expt 5.0d0 2.0d0) > >> > #C(24.999999999999996d0 6.122503540205582d15) > >> > > >> > Somewhat surprising results since the correct answers are real > >> > numbers. > > Bruno> Raymond, you know what a floatingpoint number stands for: an interval > Bruno> of real numbers. 2.0 stands for the interval [2  eps, 2 + 2*eps], > Bruno> where eps is singlefloatnegativeepsilon. > > We've had these discussions before and I think this is were we've > always disagreed. You want to treat all FP numbers as intervals. I > want to treat them as exact rationals. this is wrong, even thought it might look good. if floats were exact rationals, then the product of two single floats would have been a double float. the issue here is accuracy vs precision, see <http://clisp.cons.org/impnotes.html#flocont>;  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.palestinecentral.com/links.html>; Of course, I haven't tried it. But it will work.  Isaak Asimov 
From: Raymond Toy <toy@rt...>  20020828 17:15:48

>>>>> "Sam" == Sam Steingold <sds@...> writes: Sam> this is wrong, even thought it might look good. Sam> if floats were exact rationals, then the product of two single floats Sam> would have been a double float. Yes, I was being sloppy with wordsfloats are not exact rationals. I just want things computed as accurately as possible given the float that was given, without trying to pretend to know what the actual accuracy of the float is. Ray 
From: Sam Steingold <sds@gn...>  20020828 17:22:21

> * In message <4nlm6q29c6.fsf@...> > * On the subject of "Re: expt bug" > * Sent on 28 Aug 2002 13:15:37 0400 > * Honorable Raymond Toy <toy@...> writes: > > >>>>> "Sam" == Sam Steingold <sds@...> writes: > > Sam> this is wrong, even thought it might look good. if floats > Sam> were exact rationals, then the product of two single floats > Sam> would have been a double float. > > Yes, I was being sloppy with wordsfloats are not exact rationals. > I just want things computed as accurately as possible given the float > precisely, not accurately! :) CLISP has no idea about accuracy! > that was given, without trying to pretend to know what the actual > accuracy of the float is. exactly! (sorry for nitpicking! :)  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.palestinecentral.com/links.html>; Your mouse has moved  WinNT has to be restarted for this to take effect. 
From: Raymond Toy <toy@rt...>  20020828 17:44:55

>>>>> "Sam" == Sam Steingold <sds@...> writes: >> * In message <4nlm6q29c6.fsf@...> >> * On the subject of "Re: expt bug" >> * Sent on 28 Aug 2002 13:15:37 0400 >> * Honorable Raymond Toy <toy@...> writes: >> >> >>>>> "Sam" == Sam Steingold <sds@...> writes: >> Sam> this is wrong, even thought it might look good. if floats Sam> were exact rationals, then the product of two single floats Sam> would have been a double float. >> >> Yes, I was being sloppy with wordsfloats are not exact rationals. >> I just want things computed as accurately as possible given the float > precisely, not accurately! :) Sam> CLISP has no idea about accuracy! Yes, precisely! >> that was given, without trying to pretend to know what the actual >> accuracy of the float is. Sam> exactly! Sam> (sorry for nitpicking! :) No problem. We should strive to be as precise and accurate as possible. :) Ray 
From: Bruno Haible <haible@il...>  20020828 17:30:51

Raymond Toy writes: > You want to treat all FP numbers as intervals. Yes. > I want to treat them as exact rationals. This point of view is justified in C or Fortran  which lack builtin RATIONAL number type, but in Lisp and Scheme, where we have the distinction between exact rationals and inexact floatingpoint numbers. > Pretending that every FP number is an interval seems somewhat wrong, > because when I write 2.0, maybe it's really from the interval [1.5, > 2.5). You can't know thatonly I know that. The ideal floating point data type contains the midpoint of the interval and an estimation of the length of the interval. Some computer algebra systems get this right. In Lisp, we only have 4 possible values for the interval length: short/single/double/long float precision. Still the precision of a rational is *way* different, because it is an infinitely short interval. > I think computer should compute things as accurately as possible, > without trying to guess what the numbers might have meant. I think computers should keep track about the accuracy. We have a rule that rationals may be replaced by floats during computations (CLHS 12.1.3.3. Rule of Float Substitutability), but never the converse. > produce the even less biased result of #c(25.0 0.0). :) Why not :) Bruno 
From: Sam Steingold <sds@gn...>  20020828 17:47:30

> * In message <15725.2367.707619.426819@...> > * On the subject of "Re: expt bug" > * Sent on Wed, 28 Aug 2002 19:32:47 +0200 (CEST) > * Honorable Bruno Haible <haible@...> writes: > > > I think computer should compute things as accurately as possible, > > without trying to guess what the numbers might have meant. > > I think computers should keep track about the accuracy. how can they without explicit interval arithmetic? > > produce the even less biased result of #c(25.0 0.0). :) > Why not :) right now: (/ (imagpart (expt 5.0 2.0)) singlefloatepsilon) ==> 73.33554 this means that the computations are __VERY__ impresize. (7400% error!!!)  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.palestinecentral.com/links.html>; Life is a sexually transmitted disease with 100% mortality. 
From: Raymond Toy <toy@rt...>  20020828 17:48:50

>>>>> "Bruno" == Bruno Haible <haible@...> writes: [snipped parts that I agree with] >> I think computer should compute things as accurately as possible, >> without trying to guess what the numbers might have meant. Bruno> I think computers should keep track about the accuracy. Yes. But we don't really have that yet, unless we do interval arithmetic by hand. (I understand interval arithmetic has it's own problems, sometimes producing bounds much wider than necessary, but I don't follow that field at all.) Ray 
From: Bruno Haible <haible@il...>  20020828 17:33:05

Sam writes: > sure, computing (sin(x)/x)^2 at y=pi/2x is the right decision here. > are you going to do this? If you don't beat me at it: yes. Bruno 
From: Sam Steingold <sds@gn...>  20020828 17:52:50

> * In message <15725.2501.307234.693664@...> > * On the subject of "Re: expt bug" > * Sent on Wed, 28 Aug 2002 19:35:01 +0200 (CEST) > * Honorable Bruno Haible <haible@...> writes: > > Sam writes: > > > sure, computing (sin(x)/x)^2 at y=pi/2x is the right decision here. > > are you going to do this? > > If you don't beat me at it: yes. I do not want a race. When are you going to do it? (today? tomorrow? next week? next year?)  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.palestinecentral.com/links.html>; I may be getting older, but I refuse to grow up! 
From: Sam Steingold <sds@gn...>  20020905 22:10:22

> * In message <4nznv82ec1.fsf@...> > * On the subject of "expt bug" > * Sent on 27 Aug 2002 17:15:26 0400 > * Honorable Raymond Toy <toy@...> writes: > > In 2.28 and 2.29 on Solaris I get: > > [1]> (expt 5.0 2.0) > #C(25.000002 4.3711393E6) > [2]> (expt 5.0d0 2.0d0) > #C(24.999999999999996d0 6.122503540205582d15) The reason is that CLISP computes the transcendental functions with higher precision than the argument is but discards the extra digits before returning the value. The problem is that CLISP computes (expt x y) as (exp (* y (log x))) and it discards the extra intermediate digits of (log x). The fix is straightforward but extensive (many functions have to be carefully modified) This is a great opportunity to get into CLISP development. Any volunteers?  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.palestinecentral.com/links.html>; There are 10 kinds of people: those who count in binary and those who do not. 