From: <lut...@fr...> - 2012-05-28 15:39:28
|
Hi, the mentioned test is tagged as "Expected failure", currently only on x86-64 (but also fails on PPC, see lp#964348 -- I agree with Nikodemus that that is expected). I am getting earnestly annoyed by this, so I investigated. It turns out things are completely different. The test is about range reduction for trigonometric functions on large float arguments, like (sin (expt 2d0 70)). It was introduced with commit ad9090dc91fc922a with a fix for lp#327192 where Gabor had noticed that on x86 the above expression returns 0d0 and Paul committed a change introducing range reduction into the relevant VOPs on x86. In short: The test is bogus. Even on x86 it only doesn't fail due to an (un)lucky choice of test values. Trigonometric functions on large float arguments on x86 still deliver completely wrong results. The x86-64 port gets everything right but the test punishes that. Basically, x86 gets things wrong because it does range reduction (even for arguments below (expt 2 63), it seems) modulo some limited-precision approximation of 2 pi while libm uses an algorithm that delivers results as if the infinitely precise pi was used to reduce mod 2 pi. So what do we want: - Document that x86 is imprecise. - Completely revert the fix for lp#327192 on the account that delivering 0 or some random value are both equally bad (the bug that "tan" crashes can be fixed independently; there is only an "fldz" missing). - Change the test to test for correct results and mark x86 as failing. - Drop the test as nobody would be using arguments this large to trigonometric functions anyway. - Implement a better range reduction on x86 (in Lisp, presumably) and use it for arguments above (expt 2 20) or so. - Convert x86 to use libm for these functions. Any comments appreciated! Some details: Here are some results of (sin (expt 2d0 n)) on SBCL x86, respectively (coerce (sin (expt 2l0 n)) 'double-float) on Clisp with the long-float precision set to 3000 bits (the results of the first expression on SBCL x86-64 are equal to Clisp's results and I assume that both are mathematically correct within one ULP or so): n = 20: SBCL x86: 0.33049314002173596d0 Clisp: 0.3304931400217347d0 n = 40: SBCL x86: -0.4057050128266265d0 Clisp: -0.40570501153282873d0 n = 70: SBCL x86: 0.936042960761405d0 Clisp: -0.9981794021933068d0 So, with n = 70, x86 does not even get the sign of the result right. Here is why the test only accidentally works on x86: For several values it compares for example (sin value) with (sin (mod value (* 2 pi))). It turns out that for the only value large enough to have a chance to make "almost=" in the test fail (in light of the size of the errors seen above), namely (* pi (expt 2d0 70)), the modular reductions by mod (Lisp) and by fprem/fsin (x87) yield the same, namely 0d0: (All following expressions are evaluated on x86) * (sin (* pi (expt 2d0 70))) 0.0d0 * (mod (* pi (expt 2d0 70)) (* 2 pi)) 0.0d0 * (sin (mod (* pi (expt 2d0 70)) (* 2 pi))) 0.0d0 But for some other values these are not equal (or "almost=") anymore (and still both completely wrong, I need not mention): * (let* ((value (* 1.625d0 pi (expt 2d0 70))) (mod (mod value (* 2 pi)))) (values (sin value) mod (sin mod))) 0.1676180272408028d0 0.0d0 0.0d0 * (let* ((value (* 1.6d0 (expt 2d0 70))) (mod (mod value (* 2 pi)))) (values (sin value) mod (sin mod))) 0.5720626950150568d0 262144.0d0 -0.08410702780950105d0 Greetings, Lutz |