This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: libm-test.inc: Computing ulps near FP_ZERO.
- From: "Carlos O'Donell" <carlos at redhat dot com>
- To: Rich Felker <dalias at aerifal dot cx>
- Cc: Brooks Moses <brooks_moses at mentor dot com>, GNU C Library <libc-alpha at sourceware dot org>, "Joseph S. Myers" <joseph at codesourcery dot com>, Andreas Jaeger <aj at suse dot com>, Thomas Schwinge <thomas at codesourcery dot com>, David Miller <davem at davemloft dot net>
- Date: Tue, 09 Apr 2013 13:34:47 -0400
- Subject: Re: libm-test.inc: Computing ulps near FP_ZERO.
- References: <5160A4D7 dot 9080802 at codesourcery dot com> <516183BD dot 6060700 at redhat dot com> <5161AADF dot 2040403 at codesourcery dot com> <5161FA1A dot 5060100 at redhat dot com> <51620F35 dot 5080101 at codesourcery dot com> <51634CA1 dot 4050006 at redhat dot com> <51636912 dot 10702 at codesourcery dot com> <5163710F dot 9080500 at mentor dot com> <20130409020954 dot GT20323 at brightrain dot aerifal dot cx> <51643640 dot 5090900 at redhat dot com> <20130409160141 dot GX20323 at brightrain dot aerifal dot cx>
On 04/09/2013 12:01 PM, Rich Felker wrote:
> On Tue, Apr 09, 2013 at 11:39:44AM -0400, Carlos O'Donell wrote:
>> On 04/08/2013 10:09 PM, Rich Felker wrote:
>>> I'm a bit confused about this whole discussion. If the test is
>>> cos(M_PI_2) then according to Wolfram Alpha,
>>> cos(1.5707963267948965579989817342720925807952880859375), where that
>>> argument is the exact value of M_PI_2, evaluates to approximately:
>>>
>>> 6.123233995736765886130329661375001... Ã 10^-17
>>
>> Could you do me a favour and compute:
>>
>> cpow (0.75 + 1.25 i, 0.75 + 1.25 i)?
>
> Wolfram Alpha says:
>
> 0.11750629391447355542027983221042048264448457971691710227141... +
> 0.34655274770833867648302535206041800107792019145641861264503... i
>
> (http://www.wolframalpha.com/input/?i=%280.75%2B1.25i%29%5E%280.75%2B1.25i%29)
Awesome. Thanks, that matches glibc's expectations.
I'd never used Wolfram Alpha before, but it's certainly awesome.
So this:
@@ -5330,7 +5356,19 @@ cos_test (void)
TEST_f_f (cos, M_PI_6l * 2.0, 0.5);
TEST_f_f (cos, M_PI_6l * 4.0, -0.5);
- TEST_f_f (cos, M_PI_2l, 0);
+ /* The value of M_PI_2l is never exactly PI/2, and therefore the
+ answer is never exactly zero. The answer is equal to the error
+ in rounding PI/2 for the type used. Thus the answer is unique
+ to each type. */
+#ifdef TEST_FLOAT
+ TEST_f_f (cos, M_PI_2l, -4.371139000186241438857289400265215e-8L);
+#endif
+#if defined TEST_DOUBLE
+ TEST_f_f (cos, M_PI_2l, 6.123233995736765886130329661375001e-17L);
+#endif
+#if defined TEST_LDOUBLE
+ TEST_f_f (cos, M_PI_2l, -2.50827880633416601177866354016537e-20L);
+#endif
TEST_f_f (cos, 0.75L, 0.731688868873820886311838753000084544L);
And this:
@@ -12133,8 +12182,20 @@ sincos_test (void)
TEST_extra (sincos, plus_infty, qnan_value, qnan_value, INVALID_EXCEPTION);
TEST_extra (sincos, minus_infty, qnan_value, qnan_value, INVALID_EXCEPTION);
TEST_extra (sincos, qnan_value, qnan_value, qnan_value);
+ /* The value of M_PI_2l is never exactly PI/2, and therefore the
+ answer is never exactly zero. The answer is equal to the error
+ in rounding PI/2 for the type used. Thus the answer is unique
+ to each type. */
+#ifdef TEST_FLOAT
+ TEST_extra (sincos, M_PI_2l, 1, -4.371139000186241438857289400265215e-8L);
+#endif
+#if defined TEST_DOUBLE
+ TEST_extra (sincos, M_PI_2l, 1, 6.123233995736765886130329661375001e-17L);
+#endif
+#if defined TEST_LDOUBLE
+ TEST_extra (sincos, M_PI_2l, 1, -2.50827880633416601177866354016537e-20L);
+#endif
- TEST_extra (sincos, M_PI_2l, 1, 0);
TEST_extra (sincos, M_PI_6l, 0.5, 0.86602540378443864676372317075293616L);
TEST_extra (sincos, M_PI_6l*2.0, 0.86602540378443864676372317075293616L, 0.5);
TEST_extra (sincos, 0.75L, 0.681638760023334166733241952779893935L, 0.731688868873820886311838753000084544L);
Fix the cos and sincos regressions.
Sadly the cpow fix of the same form:
@@ -5651,11 +5689,22 @@ cpow_test (void)
TEST_cc_c (cpow, 1, 0, 0, 0, 1.0, 0.0);
TEST_cc_c (cpow, 2, 0, 10, 0, 1024.0, 0.0);
- TEST_cc_c (cpow, M_El, 0, 0, 2 * M_PIl, 1.0, 0.0);
+ /* The value of M_PIl is never exactly PI, and therefore the
+ answer is never exactly zero. The answer is equal to the error
+ in rounding PI for the type used. Thus the answer is unique
+ to each type. The same applies to M_El. */
+#ifdef TEST_FLOAT
+ TEST_cc_c (cpow, M_El, 0, 0, 2 * M_PIl, 0.999999999999999872617861385637398697901594693607738138L, 1.596133694991510372437756115619381943355320422428308993e-8L);
+#endif
+#if defined TEST_DOUBLE
+ TEST_cc_c (cpow, M_El, 0, 0, 2 * M_PIl, 0.99999999999999999999999999999983233080834759314891457L, 5.79084090011816593139976536980817613270685802172233262e-16L);
+#endif
+#if defined TEST_LDOUBLE
+ TEST_cc_c (cpow, M_El, 0, 0, 2 * M_PIl, 0.99999999999999999999999999999999999996691529119356675L, 2.57234168828455786661626197956089780140664404404216932e-19L);
+#endif
TEST_cc_c (cpow, 2, 3, 4, 0, -119.0, -120.0);
Shows that the result is an order of magnitude worse than it should be:
Regenerating ULPs for /home/carlos/build/glibc/math/test-ifloat
testing float (inline functions)
Failure: Test: Imaginary part of: cpow (e + 0 i, 0 + 2 * M_PIl i) == 0.999999999999999872617861385637398697901594693607738138 + 1.596133694991510372437756115619381943355320422428308993e-8 i
Result:
is: 1.74845553146951715462e-07 0x1.777a5c00000000000000p-23
should be: 1.59613371408795501338e-08 0x1.1236b400000000000000p-26
difference: 1.58884219558785844129e-07 0x1.55338600000000000000p-23
ulp(x) : 1.77635683940025046468e-15 0x1.00000000000000000000p-49
ulp : 89443864.0000
max.ulp : 0.0000
Maximal error of real part of: cpow
is : 5 ulp
accepted: 5 ulp
Maximal error of imaginary part of: cpow
is : 89443864 ulp
accepted: 2 ulp
Test suite completed:
6599 test cases plus 5995 tests for exception flags executed.
2 errors occurred.
These ulp are much better than before and probably the *truth*
when it comes to glibc's cpow. That is to say that the imaginary
result is inaccurate because cexp and clog are inaccurate.
Cheers,
Carlos.