Fw: Error in documentation for tan(), and poor implementation of sin() and cos()

"Carlos O'Donell" <[email protected]> wrote:

> On 7/7/26 10:15 AM, Stefan Kanthak wrote:
>> Hi maintainers,
>>
>> <https://sourceware.org/glibc/manual/latest/html_node/Trig-Functions.html>
>> states:
>
> Stefan,
>
> Thank you for raising this issue.
>
> May you please report this to [email protected] or report a bug
> in bugzilla with component glibc?
>
> Please see:
> https://sourceware.org/glibc/wiki/FilingBug

This page only shows "Dein Browser wird geprüft!" and stalls.
So my two part bug report goes here, per mail.

PART 1
~~~~~~

<https://sourceware.org/glibc/manual/latest/html_node/Trig-Functions.html>
states:

| Mathematically, the tangent function has singularities at odd multiples
| of pi/2. If the argument x is too close to one of these singularities,
| tan will signal overflow.

There exists but no (double-precision) floating-point number for which tan()
overflows or cos() yields 0, so signalling overflow respectively underflow
is (or would be) wrong!

Shown by William Kahan MANY years ago, 6381956970095103 * 2**797 is THE
double-precision floating-point number closest to a multiple of pi/2 -- it
is about 4.687165924254627611e-19 or 0x1.14AE72E6BA22Fp-61 smaller than pi/2.

From the identities
    tan(x) = sin(x) / cos(x),
    sin(x + pi) = cos(x + pi/2),
    cos(x + pi) = -sin(x + pi/2)
and the expansion of the Madhava-Newton series for sin(x) and cos(x),
    cos(x) = x**0 / 0! - x**2 / 2! + ...
           = 1 - x**2 / 2 + ...
    sin(x) = x**1 / 1! - x**3 / 3! + ...
           = x - x**3 / 6 + ...
follows for double-precision floating-point numbers r,
r**2 / 2 + 1 = 1 or |r| < 2**-25.5 ~ 2.107342425544701589e-8e-8,
    cos(r) = -1
    sin(r) = r
    tan(r) = -r
as well as
    cos(pi/2 ± r) = -sin ±r
    sin(pi/2 ± r) = 1
    tan(pi/2 ± r) = -1 / ±r

The maximum absolute value for the double-precision tan() is therefore about
2.133485385753703844e+18, i.e. 290 orders of magnitude beyond overflow!

Evaluation of the maximum/minimum values for long double is left as an
exercise.


PART 2
~~~~~~

>> <https://sourceware.org/glibc/manual/latest/html_node/Errors-in-Math-Functions.html>
>> states:
>>
>> | . Each function with a floating-point result behaves as if it computes an
>> |   infinite-precision result that is within a few ulp of the mathematically
>> |   correct value of the function [...]
>>
>> How much are "a few ulp"? Does 179 count as "few"?
>
> Generally <10 ULP. No, 179 does not count as a few. Though there are known outliers
> that are more than 10 ULP.
>
> Please report them as bugs.

<https://godbolt.org/noscript/z/M9Ks87Wx9> is a slightly bigger demonstration -- it
feeds integers which are near integral multiples of pi, pi/2 or pi/4 and exactly
representable as double-precision floating-point numbers to cos(), sin() and tan()
and prints a line if the value computed during runtime, i.e. by GLIBC, differs by
more than 2 ULP from the value computed during compile time, i.e. by GCC.

It shows errors of 269 ULP for 72 sin/tan, 262 ULP for 66 tan, and 179 ULP for 66 cos!

  61 tan
      65398140378926      -68524021915772.6328  -0x1.f293efe103e51p+45  -0x1.f293efe103e4fp+45      -68524021915772.6172
  62 cos
          74357078147863  -7.42638965257211217e-15  -0x1.0b905acac1b53p-47  -0x1.0b905acac1b51p-47  -7.42638965257210901e-15
  62 tan            74357078147863       134654932852015.438   0x1.e9df2dc274bdcp+46    0x1.e9df2dc274bep+46
134654932852015.5
  63 sin
         139755218526789  -7.16703280049355271e-15  -0x1.023835bd45532p-47  -0x1.023835bd45536p-47  -7.16703280049355902e-15
  63 tan           139755218526789   7.16703280049355271e-15   0x1.023835bd45532p-47   0x1.023835bd45536p-47
7.16703280049355902e-15
  65 sin           148714156295726   1.48527793051442243e-14   0x1.0b905acac1b53p-46   0x1.0b905acac1b51p-46
1.4852779305144218e-14
  65 tan
        148714156295726  -1.48527793051442243e-14  -0x1.0b905acac1b53p-46  -0x1.0b905acac1b51p-46   -1.4852779305144218e-14
  66 cos           214112296674652   2.59356852078558913e-16   0x1.2b04a1af8c415p-52   0x1.2b04a1af8c362p-52
2.59356852078550088e-16
  66 tan           214112296674652          3855691461342618   0x1.b65763fd56b34p+51   0x1.b65763fd56c3ap+51
3855691461342749
  67 sin           279510437053578   1.43340656009871054e-14   0x1.023835bd45532p-46   0x1.023835bd45536p-46
1.4334065600987118e-14
  67 tan           279510437053578   1.43340656009871054e-14   0x1.023835bd45532p-46   0x1.023835bd45536p-46
1.4334065600987118e-14
  68 sin
       288469374822515  -7.68574650465067005e-15  -0x1.14e87fd83e173p-47  -0x1.14e87fd83e16cp-47    -7.685746504650659e-15
  68 tan
       288469374822515  -7.68574650465067005e-15  -0x1.14e87fd83e173p-47  -0x1.14e87fd83e16cp-47    -7.685746504650659e-15
  69 sin           428224593349304   5.18713704157117826e-16   0x1.2b04a1af8c415p-51   0x1.2b04a1af8c362p-51
5.18713704157100176e-16
  69 tan
         428224593349304  -5.18713704157117826e-16  -0x1.2b04a1af8c415p-51  -0x1.2b04a1af8c362p-51  -5.18713704157100176e-16
  70 sin           567979811876093   6.64831909633643459e-15   0x1.df0fd744991e1p-48   0x1.df0fd744991ffp-48
6.64831909633645826e-15
  70 tan           567979811876093   6.64831909633643459e-15   0x1.df0fd744991e1p-48   0x1.df0fd744991ffp-48
6.64831909633645826e-15
  71 sin
         856449186698608  -1.03742740831423565e-15  -0x1.2b04a1af8c415p-50  -0x1.2b04a1af8c362p-50  -1.03742740831420035e-15
  71 tan
         856449186698608  -1.03742740831423565e-15  -0x1.2b04a1af8c415p-50  -0x1.2b04a1af8c362p-50  -1.03742740831420035e-15
  72 sin          1284673780047912   1.55614111247135338e-15   0x1.c086f2875261fp-50   0x1.c086f28752513p-50
1.55614111247130053e-15
  72 tan
        1284673780047912  -1.55614111247135338e-15  -0x1.c086f2875261fp-50  -0x1.c086f28752513p-50  -1.55614111247130053e-15

>> With -DLIBM, i.e. sin() and cos() evaluated during runtime,
>> <https://godbolt.org/noscript/z/x33TPr6vE> yields the following results:
>>
>>   5.31937264832654142e+255  -4.68716592425461995e-19  -0x1.14ae72e6ba227p-61
>>            214112296674652   2.59356852078558913e-16   0x1.2b04a1af8c415p-52
>>                                                                      ~~~
>>             74357078147863  -7.42638965257211217e-15  -0x1.0b905acac1b53p-47
>>             65398140378926   1.45934224530656633e-14   0x1.06e4484403842p-46
>>
>>            139755218526789  -7.16703280049355271e-15  -0x1.023835bd45532p-47
>>            428224593349304   5.18713704157117826e-16   0x1.2b04a1af8c415p-51
>>                                                                      ~~~
>>            856449186698608  -1.03742740831423565e-15  -0x1.2b04a1af8c415p-50
>>                                                                      ~~~
>>
>> Without -DLIBM, i.e. when GCC evaluates the functions at compile time, the
>> results are:
>>
>>   5.31937264832654142e+255  -4.68716592425462765e-19  -0x1.14ae72e6ba22fp-61
>>            214112296674652   2.59356852078550088e-16   0x1.2b04a1af8c362p-52
>>                                                                      ~~~
>>             74357078147863  -7.42638965257210901e-15  -0x1.0b905acac1b51p-47
>>             65398140378926   1.45934224530656665e-14   0x1.06e4484403843p-46
>>
>>            139755218526789  -7.16703280049355902e-15  -0x1.023835bd45536p-47
>>            428224593349304   5.18713704157100176e-16   0x1.2b04a1af8c362p-51
>>                                                                      ~~~
>>            856449186698608  -1.03742740831420035e-15  -0x1.2b04a1af8c362p-50
>>                                                                      ~~~
>>
>> The values computed by GCC are correct, the underlined values computed by
>> GLIBC are 179 ULP off!
>
> Thank you.
>
> -- 
> Cheers,
> Carlos.