This is the mail archive of the glibc-bugs@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[Bug math/15359] Inaccurate cos for IBM long double


http://sourceware.org/bugzilla/show_bug.cgi?id=15359

Carlos O'Donell <carlos at redhat dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |WAITING

--- Comment #3 from Carlos O'Donell <carlos at redhat dot com> 2013-04-12 06:53:33 UTC ---
OK I have a PPC system to test this with setup now.

pi/2 rounded to an IBM long double is:
1.570796326794896619231321691639740613531845480547870673025773690344164545962257761857472360134124755859375L

An infinite precision cos(x) of this value yields:
cos(1.570796326794896619231321691639740613531845480547870673025773690344164545962257761857472360134124755859375)=
1.0828566739219139682237461698605809743657180846737456545052536933566509276210257011310354810728258977667687179060601794311587878461976409993637130916602872593217973444675992784631595571310528753415771366618685685565006009288560082552130...
à 10^-32

The highest precision answer in an IBM long double is:
1.082856673921913968223746169860583539424442830597581817048602989306514317306816102445404953931906744377790888342946371404131442949282377978537095231104103731922805309295654296875e-32

Notice that the infinite precision answer differs from rounded IBM long double
value around ~1e-64? That's no mistake, that's around 1ulp. So it should be
possible to attain this value.

However the ppc cos(x) implemetnation yields:
1.0828566739219139682237461698548710058685937356711019144992345158825349842613103435536910382894240130262349699496720470547777605274664125545314163900911808013916015625e-32

They differ at around the 29th digit, or a difference of
5.71253355584909492648e-62.

Comparing hex digits:
Is:        0x1.c1cd129024e088a67cc7402000p-107
Should be: 0x1.c1cd129024e088a67cc74020bcp-107

It shows the answer is not as precise as the IBM long double type would allow.
A rough guess is that we have 8-12 bits lost.

Therefore there is nothing wrong with the test, it is now correctly showing
that there are 376 ulps of error in the IBM long double implementation of
cos(x). Previously we compared to zero and used a huge value for ulp, ignoring
a lot of inaccuracy around zero.

Note: At the expected value for cos(~pi/2) the value of 1 ulp is 1.5193e-64

Do you see any flaws in that analysis?

The patch I'll probably apply is to convert all the results to hex literals
with the maximum precision for the type.

e.g.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 0049fcd..f0eb81a 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5337,23 +5337,23 @@ cos_test (void)
      to each type.  */
 #ifdef TEST_FLOAT
   /* 32-bit float.  */
-  TEST_f_f (cos, M_PI_2l, -4.371139000186241438857289400265215e-8L);
+  TEST_f_f (cos, M_PI_2l, -0x1.777a5cp-25L);
 #endif
 #if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MANT_DIG == 53)
   /* 64-bit double or 64-bit long double.  */
-  TEST_f_f (cos, M_PI_2l, 6.123233995736765886130329661375001e-17L);
+  TEST_f_f (cos, M_PI_2l, 0x1.1a62633145c07p-54L);
 #endif
 #if defined TEST_LDOUBLE && LDBL_MANT_DIG == 64
   /* 96-bit long double.  */
-  TEST_f_f (cos, M_PI_2l, -2.50827880633416601177866354016537e-20L);
+  TEST_f_f (cos, M_PI_2l, -0xe.ce675d1fc8f8cbbp-69L);
 #endif
 #if defined TEST_LDOUBLE && LDBL_MANT_DIG == 106
   /* 128-bit IBM long double.  */
-  TEST_f_f (cos, M_PI_2l, 1.082856673921913968223746169860580e-32L);
+  TEST_f_f (cos, M_PI_2l, 0x1.c1cd129024e088a67cc74020bcp-107L);
 #endif
 #if defined TEST_LDOUBLE && LDBL_MANT_DIG == 113
   /* 128-bit long double.  */
-  TEST_f_f (cos, M_PI_2l, 4.335905065061890512398522013021675e-35L);
+  TEST_f_f (cos, M_PI_2l, 0x1.cd129024e088a67cc74020bbea64p-115L);
 #endif

   TEST_f_f (cos, 0.75L, 0.731688868873820886311838753000084544L);
...

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]