Bug 15563 - sincos() is incorrect for long double and large inputs on x86_64
Summary: sincos() is incorrect for long double and large inputs on x86_64
Status: RESOLVED FIXED
Alias: None
Product: glibc
Classification: Unclassified
Component: math (show other bugs)
Version: 2.18
: P2 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords:
Depends on: 13658
Blocks: 13851 13852 13854
  Show dependency treegraph
 
Reported: 2013-06-03 16:05 UTC by Carlos O'Donell
Modified: 2014-06-13 17:32 UTC (History)
5 users (show)

See Also:
Host:
Target:
Build:
Last reconfirmed:
fweimer: security-


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Carlos O'Donell 2013-06-03 16:05:43 UTC
Still broken for long double on x86/x86_64.

Tested with current sources.

+++ This bug was initially created as a clone of Bug #13658 +++

sincos() is inaccurate for large inputs on x86_64: with glibc 2.13,

#define _GNU_SOURCE
#include <stdio.h>
#include <math.h>

int main (void)
{
  volatile double x = 1.0e22;
  double s1, s2, c1;

  sincos (x, &s1, &c1);
  s2 = sin (x);
  printf ("s1 = %.17g\n", s1);
  printf ("s2 = %.17g\n", s2);
  return 0;
}

outputs:

s1 = 0.46261304076460175
s2 = -0.85220084976718879

(s2 is the correct value). I suppose that contrary to the other trig functions, glibc uses the hardware sincos instruction, which has never been meant to be used directly by a C library (the hardware elementary functions of the x86 processors were designed for small inputs, and they must not be used by code where inputs can be large, like here). The sincos() function can simply be implemented by a call to sin() and a call to cos() on this target.

Ditto for sincosf() and sincosl().

Note: x86 (32 bits) has the same problem, but it has been claimed that users don't care about correctness on this target.
Comment 1 jsm-csl@polyomino.org.uk 2013-06-03 17:01:12 UTC
I can't tell what you think the bug is here.  You give a testcase for 
double but talk about long double, and with the obvious substitutions

#define _GNU_SOURCE
#include <stdio.h>
#include <math.h>

int main (void)
{
  volatile long double x = 1.0e22;
  long double s1, s2, c1;

  sincosl (x, &s1, &c1);
  s2 = sinl (x);
  printf ("s1 = %.17Lg\n", s1);
  printf ("s2 = %.17Lg\n", s2);
  return 0;
}

I get on x86_64 (with current glibc)

s1 = -0.8522008497671888
s2 = -0.8522008497671888

which appears to be what you'd expect.  So what program are you building, 
with what options, and running with current glibc, and getting what 
results, and what would you like instead?
Comment 2 Carlos O'Donell 2013-06-03 17:47:29 UTC
(In reply to joseph@codesourcery.com from comment #1)
> I can't tell what you think the bug is here.  You give a testcase for 
> double but talk about long double, and with the obvious substitutions

I hadn't seen and could not find any discussion of this fix for long double (remember at the time ldbl-96 did not correctly do argument reduction).

The additional test case I was using was:

cat >> test.c <<EOF
#include <stdio.h>
#include <math.h>
int main() {
        long double x;

        for (x=1.0;x<1.e38;x*=2.0)
                printf("%Le %Lf\n",x,sinl(x));
        return 0;
}
EOF
gcc -Wall -pedantic -Wl,-rpath=/home/carlos/build/glibc:/home/carlos/build/glibc/math -Wl,--dynamic-linker=/home/carlos/build/glibc/elf/ld.so -o test test.c -lm

Producing the following results for long double on x86-64 with master:

1.000000e+00 0.841471
2.000000e+00 0.909297
4.000000e+00 -0.756802
8.000000e+00 0.989358
1.600000e+01 -0.287903
3.200000e+01 0.551427
6.400000e+01 0.920026
1.280000e+02 0.721038
2.560000e+02 -0.999208
5.120000e+02 0.079518
1.024000e+03 -0.158533
2.048000e+03 -0.313057
4.096000e+03 -0.594642
8.192000e+03 -0.956173
1.638400e+04 -0.559938
3.276800e+04 0.927856
6.553600e+04 0.692065
1.310720e+05 -0.999114
2.621440e+05 -0.084107
5.242880e+05 0.167618
1.048576e+06 0.330493
2.097152e+06 0.623844
4.194304e+06 0.975129
8.388608e+06 0.432248
1.677722e+07 -0.779564
3.355443e+07 -0.976517
6.710886e+07 0.420760
1.342177e+08 -0.763403
2.684355e+08 -0.986198
5.368709e+08 0.326568
1.073742e+09 -0.617326
2.147484e+09 -0.971310
4.294967e+09 -0.461987
8.589935e+09 0.819460
1.717987e+10 0.939325
3.435974e+10 -0.644430
6.871948e+10 0.985544
1.374390e+11 0.333940
2.748779e+11 -0.629540
5.497558e+11 -0.978265
1.099512e+12 -0.405705
2.199023e+12 0.741632
4.398047e+12 0.994984
8.796093e+12 -0.199069
1.759219e+13 0.390169
3.518437e+13 0.718491
7.036874e+13 0.999473
1.407375e+14 -0.064885
2.814750e+14 0.129496
5.629500e+14 0.256811
1.125900e+15 0.496397
2.251800e+15 0.861840
4.503600e+15 0.874217
9.007199e+15 -0.848926
1.801440e+16 0.897335
3.602880e+16 -0.792078
7.205759e+16 0.967000
1.441152e+17 -0.492738
2.882304e+17 0.857539
5.764608e+17 0.882269
1.152922e+18 -0.830649
2.305843e+18 0.925004
4.611686e+18 -0.702922
9.223372e+18 0.999930
1.844674e+19 0.023599
3.689349e+19 -0.047184
7.378698e+19 -0.094263
1.475740e+20 -0.187686
2.951479e+20 -0.368701
5.902958e+20 -0.685451
1.180592e+21 -0.998179
2.361183e+21 -0.120410
4.722366e+21 0.239068
9.444733e+21 0.464271
1.888947e+22 0.822404
3.777893e+22 0.935738
7.555786e+22 -0.660063
1.511157e+23 0.991692
3.022315e+23 0.255132
6.044629e+23 -0.493378
1.208926e+24 -0.858295
2.417852e+24 -0.880879
4.835703e+24 0.833914
9.671407e+24 -0.920465
1.934281e+25 0.719481
3.868563e+25 -0.999377
7.737125e+25 0.070569
1.547425e+26 -0.140786
3.094850e+26 -0.278768
6.189700e+26 -0.535435
1.237940e+27 -0.904431
2.475880e+27 -0.771696
4.951760e+27 0.981584
9.903520e+27 -0.375022
1.980704e+28 0.695303
3.961408e+28 0.999452
7.922816e+28 0.066180
1.584563e+29 -0.132069
3.169127e+29 -0.261824
6.338253e+29 -0.505382
1.267651e+30 -0.872184
2.535301e+30 -0.853307
5.070602e+30 0.889843
1.014120e+31 -0.812011
2.028241e+31 0.947848
4.056482e+31 -0.604204
8.112964e+31 0.962895
1.622593e+32 0.519724
3.245186e+32 -0.888036
6.490371e+32 -0.816591
1.298074e+33 0.942700
2.596148e+33 -0.629050
5.192297e+33 0.978003
1.038459e+34 0.408007
2.076919e+34 -0.745003
4.153837e+34 -0.993925
8.307675e+34 0.218781
1.661535e+35 -0.426961
3.323070e+35 -0.772176
6.646140e+35 -0.981295
1.329228e+36 0.377820
2.658456e+36 -0.699631
5.316912e+36 -0.999779
1.063382e+37 -0.042054
2.126765e+37 0.084034
4.253530e+37 0.167473
8.507059e+37 0.330216

This looks perfectly correct, and therefore it seems I must have made a mistake in my initial testing.

I will note that this is not fixed in 2.15, but is fixed on master.

I haven't done a thorough check to find which patch fixed this.
Comment 3 jsm-csl@polyomino.org.uk 2013-06-03 19:01:43 UTC
commit 8848d99dce1e57168a492d146f5e72195c7665a5
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Fri Mar 16 12:28:25 2012 +0000

    Implement ldbl-96 sinl / cosl / sincosl (bug 13851).