Trigonometric functions of large values lose precision

[10110111]

The following program compares gcem::sin with GCC's __builtin_sin for input of the form 2**n with n from 0 to 63. It appears that already for input of 512.0 relative error is more than 5e-13, and for larger input relative error tends to become larger. At input of 2**52 relative error reaches 0.3, and starting from 2**55 the result is simply zero.

I do admit that range reduction is hard to do, but I suggest to at least warn the users (in the README?) that the library is not intended to handle most of the range of the data types it works with.

BTW, why not implement the functions in terms of long double internally (at least when T is a builtin floating-point type)? It would already improve precision, and since it's constexpr, performance is not a problem.

#include <array>
#include <limits>
#include <iomanip>
#include <gcem.hpp>
#include <iostream>

template<unsigned N>
constexpr std::array<double, N> sinPow2_gcem()
{
    static_assert(N<=64);
    std::array<double, N> arr{};
    for(unsigned n=0;n<N;++n)
    {
        arr[n]=gcem::sin(static_cast<double>(1ull<<n));
    }
    return arr;
}

template<unsigned N>
constexpr std::array<double, N> sinPow2_gcc()
{
    static_assert(N<=64);
    std::array<double, N> arr{};
    for(unsigned n=0;n<N;++n)
    {
        arr[n]=__builtin_sin(static_cast<double>(1ull<<n));
    }
    return arr;
}

int main()
{
    constexpr auto table_gcem=sinPow2_gcem<64>();
    constexpr auto table_gcc =sinPow2_gcc <64>();
    for(unsigned n=0;n<table_gcem.size();++n)
    {
        std::cout << std::defaultfloat << "sin(2^" << std::setw(2) << n << ") = "
                  << std::setprecision(std::numeric_limits<double>::max_digits10)
                  << std::setw(21) << table_gcem[n] << " or "
                  << std::setw(21) << table_gcc[n] << ", rel_diff: "
                  << std::scientific << std::setprecision(2)
                  << std::abs((table_gcem[n]-table_gcc[n])/table_gcc[n]) << "\n";
    }
}

Output on my system with GCC version 8.3.0-6ubuntu1~18.04:

sin(2^ 0) =    0.8414709848078965 or    0.8414709848078965, rel_diff: 0.00e+00
sin(2^ 1) =   0.90929742682568182 or   0.90929742682568171, rel_diff: 1.22e-16
sin(2^ 2) =   -0.7568024953079282 or   -0.7568024953079282, rel_diff: 0.00e+00
sin(2^ 3) =   0.98935824662338168 or   0.98935824662338179, rel_diff: 1.12e-16
sin(2^ 4) =  -0.28790331666506591 or   -0.2879033166650653, rel_diff: 2.12e-15
sin(2^ 5) =   0.55142668124169159 or   0.55142668124169059, rel_diff: 1.81e-15
sin(2^ 6) =   0.92002603819679174 or   0.92002603819679063, rel_diff: 1.21e-15
sin(2^ 7) =   0.72103771050172805 or    0.7210377105017316, rel_diff: 4.93e-15
sin(2^ 8) =  -0.99920803410706305 or  -0.99920803410706271, rel_diff: 3.33e-16
sin(2^ 9) =  0.079518494012835325 or  0.079518494012876348, rel_diff: 5.16e-13
sin(2^10) =  -0.15853338004391457 or  -0.15853338004399595, rel_diff: 5.13e-13
sin(2^11) =  -0.31305701278994652 or  -0.31305701279012343, rel_diff: 5.65e-13
sin(2^12) =  -0.59464198760808085 or  -0.59464198760821463, rel_diff: 2.25e-13
sin(2^13) =  -0.95617315284297588 or   -0.9561731528431463, rel_diff: 1.78e-13
sin(2^14) =  -0.55993846566976402 or  -0.55993846566934702, rel_diff: 7.45e-13
sin(2^15) =   0.92785633341405394 or   0.92785633341392471, rel_diff: 1.39e-13
sin(2^16) =   0.69206545382222318 or   0.69206545382272322, rel_diff: 7.23e-13
sin(2^17) =  -0.99911378895071246 or  -0.99911378895077085, rel_diff: 5.84e-14
sin(2^18) = -0.084107027822479849 or -0.084107027809500717, rel_diff: 1.54e-10
sin(2^19) =   0.16761802725644473 or    0.1676180272206543, rel_diff: 2.14e-10
sin(2^20) =   0.33049314009026282 or   0.33049314002173469, rel_diff: 2.07e-10
sin(2^21) =   0.62384439947212311 or   0.62384439935862968, rel_diff: 1.82e-10
sin(2^22) =   0.97512939500607743 or   0.97512939494170703, rel_diff: 6.60e-11
sin(2^23) =   0.43224820173300005 or   0.43224820225679778, rel_diff: 1.21e-09
sin(2^24) =  -0.77956367249015968 or  -0.77956367321777775, rel_diff: 9.33e-10
sin(2^25) =  -0.97651729150386235 or  -0.97651729095092843, rel_diff: 5.66e-10
sin(2^26) =   0.42075989626084753 or   0.42075989775848105, rel_diff: 3.56e-09
sin(2^27) =  -0.76340322666929872 or  -0.76340322880198075, rel_diff: 2.79e-09
sin(2^28) =  -0.98619821358735438 or  -0.98619821183697565, rel_diff: 1.77e-09
sin(2^29) =   0.32656763928319343 or   0.32656766301856333, rel_diff: 7.27e-08
sin(2^30) =  -0.61732637553418868 or  -0.61732641504604213, rel_diff: 6.40e-08
sin(2^31) =  -0.97131013678590572 or  -0.97131017579293921, rel_diff: 4.02e-08
sin(2^32) =  -0.46198671538301844 or  -0.46198657951383493, rel_diff: 2.94e-07
sin(2^33) =   0.81946005342466022 or    0.8194597047356359, rel_diff: 4.26e-07
sin(2^34) =   0.93932460955758246 or   0.93932502694657105, rel_diff: 4.44e-07
sin(2^35) =  -0.64443221186658051 or  -0.64443035102329116, rel_diff: 2.89e-06
sin(2^36) =   0.98554449814305511 or    0.9855441071151041, rel_diff: 3.97e-07
sin(2^37) =   0.33393553242896457 or   0.33393988357522025, rel_diff: 1.30e-05
sin(2^38) =  -0.62953253790928321 or  -0.62953971111701379, rel_diff: 1.14e-05
sin(2^39) =  -0.97826467462457967 or  -0.97826480872807564, rel_diff: 1.37e-07
sin(2^40) =  -0.40568990812427469 or  -0.40570501153282873, rel_diff: 3.72e-05
sin(2^41) =   0.74176170084066606 or   0.74163206513677404, rel_diff: 1.75e-04
sin(2^42) =   0.99494505102476649 or   0.99498379417508209, rel_diff: 3.89e-05
sin(2^43) =  -0.19982643887478732 or  -0.19906887541578028, rel_diff: 3.81e-03
sin(2^44) =   0.39157600746803073 or   0.39016922335187676, rel_diff: 3.61e-03
sin(2^45) =   0.72061401054487106 or   0.71849129172091508, rel_diff: 2.95e-03
sin(2^46) =   0.99925593527723167 or   0.99947305248379947, rel_diff: 2.17e-04
sin(2^47) = -0.077080812954327838 or -0.064884736238273691, rel_diff: 1.88e-01
sin(2^48) =   0.15561499277355606 or   0.12949601773888192, rel_diff: 2.02e-01
sin(2^49) =   0.30743851458038091 or     0.256811307519207, rel_diff: 1.97e-01
sin(2^50) =   0.58509727294046221 or   0.49639651520894085, rel_diff: 1.79e-01
sin(2^51) =   0.94898461935558631 or   0.86183956388130056, rel_diff: 1.01e-01
sin(2^52) =   0.59847214410395666 or   0.87421730262363506, rel_diff: 3.15e-01
sin(2^53) =  -0.95892427466313834 or  -0.84892596481465499, rel_diff: 1.30e-01
sin(2^54) =   0.90929742682568182 or   0.89733475299759258, rel_diff: 1.33e-02
sin(2^55) =                     0 or  -0.79207844079082801, rel_diff: 1.00e+00
sin(2^56) =                     0 or   0.96699996306127067, rel_diff: 1.00e+00
sin(2^57) =                     0 or  -0.49273775680001242, rel_diff: 1.00e+00
sin(2^58) =                     0 or   0.85753897066968443, rel_diff: 1.00e+00
sin(2^59) =                     0 or   0.88226868987759099, rel_diff: 1.00e+00
sin(2^60) =                     0 or  -0.83064921763725463, rel_diff: 1.00e+00
sin(2^61) =                     0 or   0.92500446025316185, rel_diff: 1.00e+00
sin(2^62) =                     0 or  -0.70292244361920886, rel_diff: 1.00e+00
sin(2^63) =                     0 or    0.9999303766734422, rel_diff: 1.00e+00

[emsr]

On 6/3/19 5:41 AM, Ruslan Kabatsayev wrote: The following program compares |gcem::sin| with GCC's |__builtin_sin| for input of the form |2**n| with |n| from 0 to 63. It appears that already for input of |512.0| relative error is more than |5e-13|, and for larger input relative error tends to become larger. At input of |2**52| relative error reaches |0.3|, and starting from |2**55| the result is simply zero. I do admit that range reduction is hard to do, but I suggest to at least warn the users (in the README?) that the library is not intended to handle most of the range of the data types it works with. BTW, why not implement the functions in terms of |long double| internally (at least when |T| is a builtin floating-point type)? It would already improve precision, and since it's |constexpr|, performance is not a problem. #include <array> #include <limits> #include <iomanip> #include <gcem.hpp> #include <iostream> template<unsigned N> std::array<double, N>sinPow2_gcem() { static_assert(N<=64); std::array<double, N> arr; for(unsigned n=0;n<N;++n) { arr[n]=gcem::sin(static_cast<double>(1ull<<n)); } return arr; } template<unsigned N> std::array<double, N>sinPow2_gcc() { static_assert(N<=64); std::array<double, N> arr; for(unsigned n=0;n<N;++n) { arr[n]=__builtin_sin(static_cast<double>(1ull<<n)); } return arr; } int main() { const auto table_gcem=sinPow2_gcem<64>(); const auto table_gcc =sinPow2_gcc <64>(); for(unsigned n=0;n<table_gcem.size();++n) { std::cout << std::defaultfloat <<"sin(2^" <<std::setw(2) << n <<") = " <<std::setprecision(std::numeric_limits<double>::max_digits10) <<std::setw(21) << table_gcem[n] <<" or " <<std::setw(21) << table_gcc[n] <<", rel_diff: " << std::scientific <<std::setprecision(2) <<std::abs((table_gcem[n]-table_gcc[n])/table_gcc[n]) <<"\n"; } } Output on my system with GCC version 8.3.0-6ubuntu1~18.04: |sin(2^ 0) = 0.8414709848078965 or 0.8414709848078965, rel_diff: 0.00e+00 sin(2^ 1) = 0.90929742682568182 or 0.90929742682568171, rel_diff: 1.22e-16 sin(2^ 2) = -0.7568024953079282 or -0.7568024953079282, rel_diff: 0.00e+00 sin(2^ 3) = 0.98935824662338168 or 0.98935824662338179, rel_diff: 1.12e-16 sin(2^ 4) = -0.28790331666506591 or -0.2879033166650653, rel_diff: 2.12e-15 sin(2^ 5) = 0.55142668124169159 or 0.55142668124169059, rel_diff: 1.81e-15 sin(2^ 6) = 0.92002603819679174 or 0.92002603819679063, rel_diff: 1.21e-15 sin(2^ 7) = 0.72103771050172805 or 0.7210377105017316, rel_diff: 4.93e-15 sin(2^ 8) = -0.99920803410706305 or -0.99920803410706271, rel_diff: 3.33e-16 sin(2^ 9) = 0.079518494012835325 or 0.079518494012876348, rel_diff: 5.16e-13 sin(2^10) = -0.15853338004391457 or -0.15853338004399595, rel_diff: 5.13e-13 sin(2^11) = -0.31305701278994652 or -0.31305701279012343, rel_diff: 5.65e-13 sin(2^12) = -0.59464198760808085 or -0.59464198760821463, rel_diff: 2.25e-13 sin(2^13) = -0.95617315284297588 or -0.9561731528431463, rel_diff: 1.78e-13 sin(2^14) = -0.55993846566976402 or -0.55993846566934702, rel_diff: 7.45e-13 sin(2^15) = 0.92785633341405394 or 0.92785633341392471, rel_diff: 1.39e-13 sin(2^16) = 0.69206545382222318 or 0.69206545382272322, rel_diff: 7.23e-13 sin(2^17) = -0.99911378895071246 or -0.99911378895077085, rel_diff: 5.84e-14 sin(2^18) = -0.084107027822479849 or -0.084107027809500717, rel_diff: 1.54e-10 sin(2^19) = 0.16761802725644473 or 0.1676180272206543, rel_diff: 2.14e-10 sin(2^20) = 0.33049314009026282 or 0.33049314002173469, rel_diff: 2.07e-10 sin(2^21) = 0.62384439947212311 or 0.62384439935862968, rel_diff: 1.82e-10 sin(2^22) = 0.97512939500607743 or 0.97512939494170703, rel_diff: 6.60e-11 sin(2^23) = 0.43224820173300005 or 0.43224820225679778, rel_diff: 1.21e-09 sin(2^24) = -0.77956367249015968 or -0.77956367321777775, rel_diff: 9.33e-10 sin(2^25) = -0.97651729150386235 or -0.97651729095092843, rel_diff: 5.66e-10 sin(2^26) = 0.42075989626084753 or 0.42075989775848105, rel_diff: 3.56e-09 sin(2^27) = -0.76340322666929872 or -0.76340322880198075, rel_diff: 2.79e-09 sin(2^28) = -0.98619821358735438 or -0.98619821183697565, rel_diff: 1.77e-09 sin(2^29) = 0.32656763928319343 or 0.32656766301856333, rel_diff: 7.27e-08 sin(2^30) = -0.61732637553418868 or -0.61732641504604213, rel_diff: 6.40e-08 sin(2^31) = -0.97131013678590572 or -0.97131017579293921, rel_diff: 4.02e-08 sin(2^32) = -0.46198671538301844 or -0.46198657951383493, rel_diff: 2.94e-07 sin(2^33) = 0.81946005342466022 or 0.8194597047356359, rel_diff: 4.26e-07 sin(2^34) = 0.93932460955758246 or 0.93932502694657105, rel_diff: 4.44e-07 sin(2^35) = -0.64443221186658051 or -0.64443035102329116, rel_diff: 2.89e-06 sin(2^36) = 0.98554449814305511 or 0.9855441071151041, rel_diff: 3.97e-07 sin(2^37) = 0.33393553242896457 or 0.33393988357522025, rel_diff: 1.30e-05 sin(2^38) = -0.62953253790928321 or -0.62953971111701379, rel_diff: 1.14e-05 sin(2^39) = -0.97826467462457967 or -0.97826480872807564, rel_diff: 1.37e-07 sin(2^40) = -0.40568990812427469 or -0.40570501153282873, rel_diff: 3.72e-05 sin(2^41) = 0.74176170084066606 or 0.74163206513677404, rel_diff: 1.75e-04 sin(2^42) = 0.99494505102476649 or 0.99498379417508209, rel_diff: 3.89e-05 sin(2^43) = -0.19982643887478732 or -0.19906887541578028, rel_diff: 3.81e-03 sin(2^44) = 0.39157600746803073 or 0.39016922335187676, rel_diff: 3.61e-03 sin(2^45) = 0.72061401054487106 or 0.71849129172091508, rel_diff: 2.95e-03 sin(2^46) = 0.99925593527723167 or 0.99947305248379947, rel_diff: 2.17e-04 sin(2^47) = -0.077080812954327838 or -0.064884736238273691, rel_diff: 1.88e-01 sin(2^48) = 0.15561499277355606 or 0.12949601773888192, rel_diff: 2.02e-01 sin(2^49) = 0.30743851458038091 or 0.256811307519207, rel_diff: 1.97e-01 sin(2^50) = 0.58509727294046221 or 0.49639651520894085, rel_diff: 1.79e-01 sin(2^51) = 0.94898461935558631 or 0.86183956388130056, rel_diff: 1.01e-01 sin(2^52) = 0.59847214410395666 or 0.87421730262363506, rel_diff: 3.15e-01 sin(2^53) = -0.95892427466313834 or -0.84892596481465499, rel_diff: 1.30e-01 sin(2^54) = 0.90929742682568182 or 0.89733475299759258, rel_diff: 1.33e-02 sin(2^55) = 0 or -0.79207844079082801, rel_diff: 1.00e+00 sin(2^56) = 0 or 0.96699996306127067, rel_diff: 1.00e+00 sin(2^57) = 0 or -0.49273775680001242, rel_diff: 1.00e+00 sin(2^58) = 0 or 0.85753897066968443, rel_diff: 1.00e+00 sin(2^59) = 0 or 0.88226868987759099, rel_diff: 1.00e+00 sin(2^60) = 0 or -0.83064921763725463, rel_diff: 1.00e+00 sin(2^61) = 0 or 0.92500446025316185, rel_diff: 1.00e+00 sin(2^62) = 0 or -0.70292244361920886, rel_diff: 1.00e+00 sin(2^63) = 0 or 0.9999303766734422, rel_diff: 1.00e+00 | ??? You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#15?email_source=notifications&email_token=AAOYYXYGT5BSQWIPLMSH4H3PYTRLXA5CNFSM4HSGKB22YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4GXHTPGQ>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAOYYXYAGONHMWGF4APURQDPYTRLXANCNFSM4HSGKB2Q>.

What is probably happening is the finite precision of pi means the period of sin is off.?? This makes a larger and larger difference the greater the argument.?? In fact you can see above that for double with a 53-bit mantissa you've lost your last significant bit and from then on the two functions are completely uncorrelated.?? You either need argument reduction or to divide by a pi of greater than working precision. BTW this is why reperiodized trig functions such as sin_pi() are nice to have: sin_pi(x) := sin(pi*x) except the implementation splits x into integer and fractional parts x = n + y and then sin_pi(x) = +-sin(pi*y).