From: Nikodemus Siivola <nikodemus@ra...>  20110324 06:07:41

Lutz Euler is the man  even though he's not subscribed, he used his powers to telepathically read the mail and figured out what's going on! Cheers,  Nikodemus  Forwarded message  From: Lutz Euler <lutz.euler@...> Date: 24 March 2011 01:49 Subject: Re: sbcl1.0.45 and maxima5.23.1 testsuite To: Nikodemus Siivola <nikodemus@...>, "Andrey G. Grozin" <A.G.Grozin@...> Cc: sbclhelp@... Nikodemus wrote: > On 17 January 2011 18:15, Andrey G. Grozin <A.G.Grozin@...> wrote: > >> Now there is a new testsuite failure with sbcl. It's rtest16, problem 385. >> It's about floatingpoint calculation of zeta(%i+3) (%i is the imaginary >> unit). sbcl produced a resulr with error of order 3*10^(9) instead of >> expected 10^(15). Earlier versions of sbcl did not have this failure. >> Looks like a regression. All the other lisps listed above also don't have >> this particular failure (ecl has 2 failures in rtest8, this is another >> story). Looks like something is calculated in single precision when it >> should be done in double. > > Do you recall what's the latest version in which that test passed? (It > fails with 1.0.44 as well.) > > Alternatively, if a Maxima wizard can reduce a testcase that doesn't > involve Maxima  just shows an SBCL operation that returns the wrong > value  what would help too. I am not a maxima wizard but I believe I have found the cause of the problem: SBCL's "expt" with integer or rational base and complex doublefloat exponent uses an intermediate singlefloat value for its calculations thus, while yielding a complex doublefloat result as required, this result is exact only to singlefloat precision. This behaviour is the same already in SBCL 0.8.12, which is the oldest version I happen to have available readily compiled. The details: Maxima's calculation of the zeta function with the argument of the test case is done using the lisp function "floatzeta", defined in maxima5.23.1/src/combin.lisp. The formula contains a sum of several terms each containing a quotient of something by "expt" of a small integer and the argument to the zeta function, the latter converted to complex doublefloat. These "expt"s are calculated directly using SBCL's "expt". Here is a trace: Maxima 5.23.1 http://maxima.sourceforge.net using Lisp SBCL 1.0.46.42 (%i1) to_lisp(); Type (tomaxima) to restart, ($quit) to quit Maxima. MAXIMA> (trace expt) (EXPT) MAXIMA> (floatzeta #$(3+%i)$) 0: (EXPT 4 3.0) 0: EXPT returned 64.0 0: (EXPT 2 #C(2.0 1.0)) 0: EXPT returned #C(0.19230972430417587 0.15974031883619208) 0: (EXPT #C(7.2421875 1.0) #C(1.25 0.5)) 0: EXPT returned #C(4.418528303829868 10.318276280305101) 0: (EXPT 1 #C(3.0 1.0)) 0: EXPT returned #C(1.0 0.0) 0: (EXPT 1 0) 0: EXPT returned 1 0: (EXPT 1 0) 0: EXPT returned 1 0: (EXPT 21 #C(3.0 1.0)) 0: EXPT returned #C(9217.405235813078 897.5555997887233) 0: (EXPT 2 #C(3.0 1.0)) 0: EXPT returned #C(6.1539112363389945 5.11169025143816) ... ... [many more calls to (expt <small int> #C(3.0 1.0))] ... #C(1.1072144060921958 0.14829086482826961) The test case expects the following result: #C(1.10721440843141  0.1482908671781754) As a simple reduced test case without maxima we have: (expt 2 #C(3.0d0 1.0d0)) SBCL gives: #C(6.1539112363389945d0 5.11169025143816d0) Clisp (2.48) gives: #C(6.153911210911777d0 5.111690210509078d0) and Clisp with more precision (setf (system::longfloatdigits) 500): (expt 2 #C(3.0l0 1.0l0)) #C(6.15391121091177701262663994929016561152716795929571082535897...l0 5.11169021050907840920026329171761427405784430264463783596469...l0) The definition of "expt" is in sbcl/src/code/irrat.lisp. (defun expt (base power) ... (numberdispatch ((base number) (power number)) ... (((foreach fixnum (or bignum ratio) singlefloat doublefloat) complex) (if (and (zerop base) (plusp (realpart power))) (* base power) (exp (* power (log base))))) ...)) The calculation of (log base), base being an integer, yields a singlefloat as required by the standard. For the purpose of "expt" it may be a bug to use (log base) as the standard says in section 12.1.4.1 "Rule of Float and Rational Contagion": "When rationals and floats are combined by a numerical function, the rational is first converted to a float of the same format." For a quick test I replaced the calculation in the SBCL source as follows, (exp (* power (log (coerce base 'doublefloat)))) then compiled SBCL and maxima. This SBCL calculates: * (expt 2 #C(3.0d0 1.0d0)) #C(6.153911210911775d0 5.111690210509077d0) This seems to differ from Clisp's result only in the last bit of the real and imaginary part each. The zeta value calculated by maxima is then: (%i1) zeta (3 + %i), numer=true; (%o1) 1.107214408431409  .1482908671781753 %i The resulting maxima passes its test suite with no unexpected errors. In case it matters, the above is under x8664, except when using SBCL 0.8.12, which is x86. If "expt" indeed needs to be modified a correct fix seems a bit more complicated as currently there are two clauses dealing with complex exponents, both not distinguishing between single and double floats. The standard seems to require to look at all argument types and to use singlefloat calculations if all arguments are integer or rational or singlefloat, be it real or complex, and to use doublefloat calculations if at least one argument is doublefloat real or complex? Also, the current implementation of "expt" allows (with a complex exponent) bignums or ratios as the base that are too big or too small to be represented as a singlefloat, provided their logarithm can (and the exponent is such that the result fits in its type). If one doesn't want to lose this property, a function for log of ratios and integers returning a doublefloat would need to be defined. This may be easy as "log" already contains the needed logic, only unfortunately always coercing the result to singlefloat. On the other hand I believe the standard does not require this, so hopefully nobody depends on this accidental property. And, when the exponent is real, bignum/ratio bases too large and too small don't work anyway, currently. On the third hand, as "log" (intentionally, seemingly) deals with bignums and ratios that don't fit in a float, maybe "expt" should, too, even in the case of a real exponent? Moreover, "log" is implemented the way it is to be more exact in the case that its argument is a ratio very near to one. This carries over to "expt": The original SBCL calculates: * (expt (/ (expt 2 64) (1+ (expt 2 64))) #c(10 10)) #C(1.0 5.421011e19) whereas with the above hack we have the less exact #C(1.0d0 0.0d0) (Sorry, can't resist: Did you notice that "log" returns zero for rational arguments very near to one if the integer length of numerator and denominator differs by one? (log (/ (expt 2 64) (1 (expt 2 64)))) > 0.0 whereas (log (/ (1+ (expt 2 64)) (expt 2 64))) > 5.421011e20 which result would be correct for the former ratio, too.) Greetings Lutz 