How would I make Taylor Series cosine function deal with larger angles

Question

I have implemented my own Cosine function based on Taylor Series and it looks like it doesn't work well with angles greater than pi. Here is my implementation:

#include <cmath>
#include <iostream>

static constexpr float PI_F32 = static_cast<float>(M_PI);
constexpr float mod(float x, float y)
{
    return (x - (int(x / y) * y));
}
constexpr float Cosine(float radians) 
{
    const float x = mod(fabsf(radians), float(2.f * PI_F32));
    const float x2 = x * x;

    constexpr float c2 = 1.f / 24.f;
    constexpr float c3 = -1.f / 720.f;
    constexpr float c4 = 1.f / 40320.f;
    constexpr float c5 = -1.f / 3628800.f;

    return 1.f - (x2 * 0.5f) + (x2 * (x2 * c2 + x2 * (x2 * c3 + x2 * (x2 * c4 + x2 * x2 * c5))));
}
int main() 
{
    static constexpr float x = 3.5f * PI_F32;
    printf("standard: %0.9f\n", cosf(x));
    printf("my version: %0.9f\n", Cosine(x));
    return 0;
}

Output

standard: 0.000000664
my version: -0.222440109

I would like to only have 5 terms and doing a mod doesn't help for some reason. I have read somewhere that this implementation is only good for angles in range [-pi, pi]. I would like to know how would I fit a greater angle in that range. I could not find any solutions and if there is already an answer please do link it in comments.

Answer 1

This is mostly about the symmetry discussed in the comments. Here's the basic idea of what the cosine is:

The radius of the circle is 1 (we don't care what unit--just 1 of something). I've marked radians around the outside, starting from 0 at the right, and proceeding counter-clockwise through Pi/2 radians at the top, Pi radians on the left, and 3/2 Pi radians at the bottom.

The blue line represents an angle starting from the 0 at the right side, that's greater than Pi radians. The red lines below depict what the cosine represents--the horizontal distance from the center of the circle to wherever the line at that angle intersects the circle. Since the circle has a radius of 1, the cosine will always be a value between 0 and 1 inclusive.

Since the cosine is the measure of the horizontal distance, if we reflect that angle across the horizontal axis so it's above the center of the circle instead of below the center, the horizontal distance will obviously remain the same. So if that angle is Pi + N radians, then Pi - N radians will have exactly the same cosine.

We can also go a step further: after we've reflected the angle so it's above the center, we can also reflect it horizontally so it's to the right of the center, though we have to negate the cosine. So, once we've moved it above the center, we can represent the angle as Pi/2 + M, and we can immediately see that cosine(Pi/2+M) = - cosine(Pi/2-M).

So we can pretty easily restrict the real computation to the range 0..Pi/2, and then if the original angle was between Pi/2 and 3/2 Pi, we know the cosine will be negated.

For what it's worth, the sine is essentially the same, except it deals with the vertical distance instead of the horizontal. So clearly it'll remain the same when we reflect an angle across the vertical axis, and change sign (but otherwise remain the same) if we reflect across the horizontal axis.

Answer 2

This is more of a calculus problem You are stopping too early on your power series expansion. You need to ensure that sum of remaining terms is negligible, compared to the current sum. In your snippet, next term would be pow(1.5*π,12)/fact(12)≈0.25, which is well comparable to the maximum value of cosine. Use geometric identities to get the angle in range [0,π/4] and do the approximation. Since max(pow(x,n)) < 1, now you can stop at the 6th term with acceptable error magnitude.

Answer 3

Another approach is to use recursion to get the base angle small enough, then switch to a simple approximant (e.g. Taylor series or Padé approximant) for relatively small x.

Any of the multi-angle trig formulae give suitable recursion formulae; e.g.

cos(x)=8 cos⁴(x/4) - 8 cos²(x/4) + 1

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

const double TWOPI = 8 * atan( 1.0 );

double cosine( double x )
{
   x = remainder( x, TWOPI );                         // remove multiples of 2.PI

   if ( abs( x ) < 0.01 )                             // if small, use polynomial approximant
   {
      double x2 = x * x;
//    return ( 12 - 5 * x2 ) / ( 12 + x2 );           // Pade approximant
      return 1 - 0.5 * x2 * ( 1 - x2 / 12 );          // Taylor series
   }
   else                                               // otherwise use recursion
   {
      double C = cosine( x / 4 ), C2 = C * C;
      return 8 * C2 * ( C2 - 1 ) + 1;
   }
}

int main()
{
   #define FMT << setw( 15 ) <<
   cout FMT "x" FMT "solution" FMT "std::cos" FMT "Error" << '\n';
   for ( int i = 0; i < 16; i++ )
   {
       double x = i * 1.0;
       cout FMT x FMT cosine( x ) FMT cos( x ) FMT abs( cos(x) - cosine(x) ) << '\n';
   }
}

Output:

          x       solution       std::cos          Error
          0              1              1              0
          1       0.540302       0.540302    8.23785e-13
          2      -0.416147      -0.416147    8.93358e-12
          3      -0.989992      -0.989992     8.1779e-13
          4      -0.653644      -0.653644    1.51822e-11
          5       0.283662       0.283662    2.08528e-12
          6        0.96017        0.96017    1.57652e-14
          7       0.753902       0.753902    1.86762e-12
          8        -0.1455        -0.1455     5.2357e-12
          9       -0.91113       -0.91113    9.83325e-13
         10      -0.839072      -0.839072    1.24134e-12
         11      0.0044257      0.0044257    3.96636e-14
         12       0.843854       0.843854    2.73004e-12
         13       0.907447       0.907447    4.76175e-13
         14       0.136737       0.136737    8.99475e-13
         15      -0.759688      -0.759688    1.73251e-11