- To: libc-gnats at gnu dot org, gnats-admin at gnu dot org
- Subject: libc/2269: triginometric argument reduction error in libm
- From: ballen at uwm dot edu
- Date: Thu, 17 May 2001 00:10:45 -0400
>Number: 2269
>Category: libc
>Synopsis: triginometric argument reduction error in libm
>Confidential: no
>Severity: serious
>Priority: high
>Responsible: libc-gnats
>State: open
>Quarter:
>Keywords:
>Class: sw-bug
>Submitter-Id: gnatsweb
>Arrival-Date: Thu May 17 00:10:44 -0400 2001
>Cases:
>Originator: Bruce Allen
>Release: current
>Organization:
U. of Wisconsin - Milwaukee
>Environment:
Intel Linux.
>Description:
Argument reduction is not being properly handled in the sin
and cos functions. Unacceptable errors (gross violation
of IEEE754 standard) occur when x>~10^18.
Further detail is in attached file.
>How-To-Repeat:
Run the following block of code on several platforms
(solaris, alpha linux with ccc & compaq math lib)
and compare output with the bc script below
/* compile with gcc demo.c -o demo -lm */
#include <stdio.h>
#include <math.h>
int main() {
float x;
for (x=1.0;x<1.e38;x*=2.0)
printf("%e %f\n",x,sin(x));
return 0;
}
/* to run: bc -l script */
scale=60
x=1.0
while (x < 10^38) {
print s(x), "\n"
x=2.0*x
}
halt
>Fix:
Do some systematic testing of argument reduction functions,
isolate problem, and fix. I'd start with
/usr/lib/libc/sysdeps/ieee754/flt-32/e_rem_pio2f.c
>Unformatted:
Column 1 shows values of x. It is always EXACTLY a power of 2 (what
is printed as "x" is a user-friendly approximation). Note that these
powers of 2 are exactly represented in IEEE754 floating-point format.
Columns 2-6 shows values of sin(x) computed in different ways (code below).
Column 2: Intel linux
Column 3: Solaris 2.8
Column 4: Alpha linux with Compaq C compiler and Compaq Extended Math Library
Column 5: Alpha linux with gcc math library (I think!)
Column 6: Arbitrary precision calculator bc
A few comments:
(A) The code that produced this output is below. It's only a few lines: the results can
be easily reproduced. Compare results when x is greater than of order 10^18.
(B) Please DON'T tell me "ahh, of course the sin/cos of a large floating point number
is intrinsically inaccurate." Go read the IEEE754 standard more
carefully and see (E) below. Every normalized IEEE754 float has an
EXACT value. The IEEE standard clearly states that math functions
should return the closest IEEE float to the true correct value.
Bottom line: columns 2-6 should agree to at least 6 decimal places.
(C) The values printed in column 1 are increasing integer powers of 2. Printing them
exactly takes up a lot of space so I've just printed the approximate value.
(D) The function /usr/lib/libc/sysdeps/ieee754/flt-32/e_rem_pio2f.c is I think the
generic argument reduction code. I can't tell if a bug slipped into it or not, or
if instead the problem is some architecture-specific "fast" routine.
(E) The original paper that described how to do trig argument reduction by storing only
about 400 binary digits of 2/pi is by Mary H. Payne and Robert N. Hanek (Digital
Equipment Corporation) ACM Signum Newsletter #16 (1983) page 19. I was showing
my students in my graduate course on mathematical methods how nice this was, and how
every standard library implements this technique. Imagine my surprise when I found
that my favorite library got it wrong.
(F) For testing purposes there is a nice paper by Roger Alan Smith, IEEE Transactions
on computers, Vol 44 No 11, November 1995 page 1348. It gives a table with
the worst-case values of x (the ones that require the most bits of 2/pi to
calculate correctly). Testing an argument reduction routine with these values
is a good idea.
(G) My code is after the table (go to the bottom).
(1) (2) (3) (4) (5) (6)
x x86-gcc Solaris AXP-CCC AXP-GCC BC
1.000000e+00 0.841471 0.841471 0.841471 0.841471 .84147098480
2.000000e+00 0.909297 0.909297 0.909297 0.909297 .90929742682
4.000000e+00 -0.756802 -0.756802 -0.756802 -0.756802 -.7568024953
8.000000e+00 0.989358 0.989358 0.989358 0.989358 .98935824662
1.600000e+01 -0.287903 -0.287903 -0.287903 -0.287903 -.2879033166
3.200000e+01 0.551427 0.551427 0.551427 0.551427 .55142668124
6.400000e+01 0.920026 0.920026 0.920026 0.920026 .92002603819
1.280000e+02 0.721038 0.721038 0.721038 0.721038 .72103771050
2.560000e+02 -0.999208 -0.999208 -0.999208 -0.999208 -.9992080341
5.120000e+02 0.079518 0.079518 0.079518 0.079518 .07951849401
1.024000e+03 -0.158533 -0.158533 -0.158533 -0.158533 -.1585333800
2.048000e+03 -0.313057 -0.313057 -0.313057 -0.313057 -.3130570127
4.096000e+03 -0.594642 -0.594642 -0.594642 -0.594642 -.5946419876
8.192000e+03 -0.956173 -0.956173 -0.956173 -0.956173 -.9561731528
1.638400e+04 -0.559938 -0.559938 -0.559938 -0.559938 -.5599384656
3.276800e+04 0.927856 0.927856 0.927856 0.927856 .92785633341
6.553600e+04 0.692065 0.692065 0.692065 0.692065 .69206545382
1.310720e+05 -0.999114 -0.999114 -0.999114 -0.999114 -.9991137889
2.621440e+05 -0.084107 -0.084107 -0.084107 -0.084107 -.0841070278
5.242880e+05 0.167618 0.167618 0.167618 0.167618 .16761802722
1.048576e+06 0.330493 0.330493 0.330493 0.330493 .33049314002
2.097152e+06 0.623844 0.623844 0.623844 0.623844 .62384439935
4.194304e+06 0.975129 0.975129 0.975129 0.975129 .97512939494
8.388608e+06 0.432248 0.432248 0.432248 0.432248 .43224820225
1.677722e+07 -0.779564 -0.779564 -0.779564 -0.779564 -.7795636732
3.355443e+07 -0.976517 -0.976517 -0.976517 -0.976517 -.9765172909
6.710886e+07 0.420760 0.420760 0.420760 0.420760 .42075989775
1.342177e+08 -0.763403 -0.763403 -0.763403 -0.763403 -.7634032288
2.684355e+08 -0.986198 -0.986198 -0.986198 -0.986198 -.9861982118
5.368709e+08 0.326568 0.326568 0.326568 0.326568 .32656766301
1.073742e+09 -0.617326 -0.617326 -0.617326 -0.617326 -.6173264150
2.147484e+09 -0.971310 -0.971310 -0.971310 -0.971310 -.9713101757
4.294967e+09 -0.461987 -0.461987 -0.461987 -0.461987 -.4619865795
8.589935e+09 0.819460 0.819460 0.819460 0.819460 .81945970473
1.717987e+10 0.939325 0.939325 0.939325 0.939325 .93932502694
3.435974e+10 -0.644430 -0.644430 -0.644430 -0.644430 -.6444303510
6.871948e+10 0.985544 0.985544 0.985544 0.985544 .98554410711
1.374390e+11 0.333940 0.333940 0.333940 0.333940 .33393988357
2.748779e+11 -0.629540 -0.629540 -0.629540 -0.629540 -.6295397111
5.497558e+11 -0.978265 -0.978265 -0.978265 -0.978265 -.9782648087
1.099512e+12 -0.405705 -0.405705 -0.405705 -0.405705 -.4057050115
2.199023e+12 0.741632 0.741632 0.741632 0.741632 .74163206513
4.398047e+12 0.994984 0.994984 0.994984 0.994984 .99498379417
8.796093e+12 -0.199069 -0.199069 -0.199069 -0.199069 -.1990688754
1.759219e+13 0.390169 0.390169 0.390169 0.390169 .39016922335
3.518437e+13 0.718491 0.718491 0.718491 0.718491 .71849129172
7.036874e+13 0.999473 0.999473 0.999473 0.999473 .99947305248
1.407375e+14 -0.064885 -0.064885 -0.064885 -0.064885 -.0648847362
2.814750e+14 0.129496 0.129496 0.129496 0.129496 .12949601773
5.629500e+14 0.256812 0.256811 0.256811 0.256811 .25681130751
1.125900e+15 0.496398 0.496397 0.496397 0.496397 .49639651520
2.251800e+15 0.861841 0.861840 0.861840 0.861840 .86183956388
4.503600e+15 0.874214 0.874217 0.874217 0.874217 .87421730262
9.007199e+15 -0.848932 -0.848926 -0.848926 -0.848926 -.8489259648
1.801440e+16 0.897325 0.897335 0.897335 0.897335 .89733475299
3.602880e+16 -0.792107 -0.792078 -0.792078 -0.792078 -.7920784407
7.205759e+16 0.966976 0.967000 0.967000 0.967000 .96699996306
1.441152e+17 -0.492899 -0.492738 -0.492738 -0.492738 -.4927377568
2.882304e+17 0.857730 0.857539 0.857539 0.857539 .85753897066
5.764608e+17 0.881919 0.882269 0.882269 0.882269 .88226868987
1.152922e+18 -0.831475 -0.830649 -0.830649 -0.830649 -.8306492176
2.305843e+18 0.923872 0.925004 0.925004 0.925004 .92500446025
4.611686e+18 -0.707133 -0.702922 -0.702922 -0.702922 -.7029224436
9.223372e+18 0.987373 0.999930 0.999930 0.999930 .99993037667
1.844674e+19 0.312821 0.023599 0.023599 0.023599 .02359850990
3.689349e+19 -0.594243 -0.047184 -0.047184 -0.047184 -.0471838762
7.378698e+19 -0.955882 -0.094263 -0.094263 -0.094263 -.0942626475
1.475740e+20 -0.561582 -0.187686 -0.187686 -0.187686 -.1876858605
2.951479e+20 0.929330 -0.368701 -0.368701 -0.368701 -.3687010302
5.902958e+20 0.686311 -0.685451 -0.685451 -0.685451 -.6854506367
1.180592e+21 -0.998319 -0.998179 -0.998179 -0.998179 -.9981794021
2.361183e+21 -0.115712 -0.120410 -0.120410 -0.120410 -.1204100802
4.722366e+21 0.229869 0.239068 0.239068 0.239068 .23906801037
9.444733e+21 0.447426 0.464271 0.464271 0.464271 .46427142695
1.888947e+22 0.800285 0.822404 0.822404 0.822404 .82240388067
3.777893e+22 0.959733 0.935738 0.935738 0.935738 .93573785320
7.555786e+22 -0.539203 -0.660063 -0.660063 -0.660063 -.6600625307
1.511157e+23 0.908207 0.991692 0.991692 0.991692 .99169201857
3.022315e+23 0.760208 0.255132 0.255132 0.255132 .25513242885
6.044629e+23 -0.987784 -0.493378 -0.493378 -0.493378 -.4933782134
1.208926e+24 0.307855 -0.858295 -0.858295 -0.858295 -.8582954304
2.417852e+24 -0.585806 -0.880879 -0.880879 -0.880879 -.8808786886
4.835703e+24 -0.949535 0.833914 0.833914 0.833914 .83391392223
9.671407e+24 -0.595666 -0.920465 -0.920465 -0.920465 -.9204650614
1.934281e+25 0.956916 0.719481 0.719481 0.719481 .71948125641
3.868563e+25 0.555710 -0.999377 -0.999377 -0.999377 -.9993765291
7.737125e+25 -0.924008 0.070569 0.070569 0.070569 .07056908810
1.547425e+26 -0.706633 -0.140786 -0.140786 -0.140786 -.1407863037
3.094850e+26 0.999999 -0.278768 -0.278768 -0.278768 -.2787681465
6.189700e+26 0.002681 -0.535435 -0.535435 -0.535435 -.5354346809
1.237940e+27 -0.005361 -0.904431 -0.904431 -0.904431 -.9044312486
2.475880e+27 -0.010723 -0.771696 -0.771696 -0.771696 -.7716958419
4.951760e+27 -0.021444 0.981584 0.981584 0.981584 .98158440403
9.903520e+27 -0.042878 -0.375022 -0.375022 -0.375022 -.3750220659
1.980704e+28 -0.085677 0.695303 0.695303 0.695303 .69530282426
3.961408e+28 -0.170723 0.999452 0.999452 0.999452 .99945178104
7.922816e+28 -0.336434 0.066180 0.066180 0.066180 .06617962947
1.584563e+29 -0.633644 -0.132069 -0.132069 -0.132069 -.1320690910
3.169127e+29 -0.980405 -0.261824 -0.261824 -0.261824 -.2618244672
6.338253e+29 -0.386262 -0.505382 -0.505382 -0.505382 -.5053817087
1.267651e+30 0.712568 -0.872184 -0.872184 -0.872184 -.8721836054
2.535301e+30 0.999880 -0.853307 -0.853307 -0.853307 -.8533072094
5.070602e+30 -0.031006 0.889843 0.889843 0.889843 .88984323544
1.014120e+31 0.061982 -0.812011 -0.812011 -0.812011 -.8120111168
2.028241e+31 0.123726 0.947848 0.947848 0.947848 .94784753153
4.056482e+31 0.245552 -0.604204 -0.604204 -0.604204 -.6042037178
8.112964e+31 0.476067 0.962895 0.962895 0.962895 .96289515929
1.622593e+32 0.837316 0.519724 0.519724 0.519724 .51972407717
3.245186e+32 0.915554 -0.888036 -0.888036 -0.888036 -.8880360820
6.490371e+32 -0.736462 -0.816591 -0.816591 -0.816591 -.8165913896
1.298074e+33 0.996402 0.942700 0.942700 0.942700 .94269950227
2.596148e+33 -0.168898 -0.629050 -0.629050 -0.629050 -.6290501714
5.192297e+33 0.332943 0.978003 0.978003 0.978003 .97800279968
1.038459e+34 0.627894 0.408007 0.408007 0.408007 .40800665745
2.076919e+34 0.977379 -0.745003 -0.745003 -0.745003 -.7450029813
4.153837e+34 0.413426 -0.993925 -0.993925 -0.993925 -.9939250685
8.307675e+34 -0.752880 0.218781 0.218781 0.218781 .21878056866
1.661535e+35 -0.991028 -0.426961 -0.426961 -0.426961 -.4269608179
3.323070e+35 0.264914 -0.772176 -0.772176 -0.772176 -.7721758248
6.646140e+35 -0.510899 -0.981295 -0.981295 -0.981295 -.9812948137
1.329228e+36 -0.878379 0.377820 0.377820 0.377820 .37782010936
2.658456e+36 -0.839669 -0.699631 -0.699631 -0.699631 -.6996314273
5.316912e+36 0.912046 -0.999779 -0.999779 -0.999779 -.9997788086
1.063382e+37 -0.748037 -0.042054 -0.042054 -0.042054 -.0420541594
2.126765e+37 0.992880 0.084034 0.084034 0.084034 .08403391098
4.253530e+37 -0.236543 0.167473 0.167473 0.167473 .16747334849
8.507059e+37 0.459661 0.330216 0.330216 0.330216 .33021611202
Here is a routine for printing out the table.
/* Compile as: gcc -o demo demo.c -lm */
#include <stdio.h>
#include <math.h>
int main() {
float x;
for (x=1.0;x<1.e38;x*=2.0)
printf("%e %f\n",x,sin(x));
return 0;
}
Here is a routine that uses the arbitrary precision calculator to
verify the correct values found by solaris, ccc, etc.
/* To run: bc -l script */
scale=60
x=1.0
while (x < 10^38) {
print s(x), "\n"
x=2.0*x
}
halt