Hello,
On Thu, 27 Apr 2023, Jakub Jelinek wrote:
> The first really large error I see is for sinl with
> x/2gx &val
> 0x748160ed90d9425b 0xefd8b811d6293294
> i.e.
> 1.5926552660973502228303666578452949e+253
> with most significant double being
> 1.5926552660973502e+253
> and low double
> -5.9963639272208416e+230
Such large numbers will always be a problem with the range reduction step
in sin/cos. With double-double the possible mantissage length can be much
larger than 106, and the range reduction needs to be precise to at
least those many bits to give anything sensible.
Realistically I think you can only expect reasonably exact results for
double-double on operations that require range-reductions for
(a) "smallish" values. Where the low double is (say) <= 2^16 * PI, or
(b) where the low double is consecutive to the high double, i.e. the
overall mantissa size (including the implicit zeros in the middle) is
less than 107 (or whatever the current precision for the
range-reduction step on large values is)
> given is
> -0.4025472157704263326278375983156912
> and expected (mpfr computed)
> -0.46994008859023245970759964236618727
> But if I try on x86_64:
> #define _GNU_SOURCE
> #include <math.h>
>
> int
> main ()
> {
> _Float128 f, f2, f3, f4;
> double d, d2;
> f = 1.5926552660973502228303666578452949e+253f128;
> d = 1.5926552660973502e+253;
> f2 = d;
> f2 += -5.9963639272208416e+230;
> f3 = sinf128 (f);
> f4 = sinf128 (f2);
> d2 = sin (d);
> return 0;
> }
> where I think f2 is what matches most closely the 106 bit precision value,
> (gdb) p f
> $7 = 1.5926552660973502228303666578452949e+253
> (gdb) p f2
> $8 = 1.59265526609735022283036665784527174e+253
> (gdb) p f3
> $9 = -0.277062522218693980443596385112227247
> (gdb) p f4
> $10 = -0.402547215770426332627837598315693221
> and f4 is much closer to the given than to expected.
Sure, but that's because f2 is only "close" to the double-double exact
value of (1.5926552660973502e+253 + -5.9963639272208416e+230) relatively,
not absolutely. As you already wrote the mantissa to represent the exact
value (which double-double and mpfr can!) is larger than 106 bits. The
necessary error of rounding to represent it in f128 is small in ULPs, but
very very large in magnitude. Large magnitude changes of input value to
sin/cos essentially put the value into completely different quadrants and
positions within those quadrants, and hence the result under such rounding
in input can be _wildly_ off.
E.g. imagine a double-double representing (2^107 * PI + PI/2) exactly
(assume PI is the 53-bit representation of pi, that's why I say
"exactly"). The correct result of sin() on that is 1. The result on the
nearest f128 input value (2^107 * PI) will be 0. So you really can't
compare f128 arithmetic with double-double one when the mantissas are too
far away.
So, maybe you want to only let your tester test "good" double-double
values, i.e. those that can be interpreted as a about-106-bit number where
(high-exp - low-exp) <= about 53.
(Or just not care about the similarities of cos() on double-double to a
random number generator :) )
Ciao,
Michael.
