CODEWITHSUNDEEP
cmathfloating-point
scifie
scifie

Reputation: 65

Calculate maclaurin series for sin using C

I wrote a code for calculating sin using its maclaurin series and it works but when I try to calculate it for large x values and try to offset it by giving a large order N (the length of the sum) - eventually it overflows and doesn't give me correct results. This is the code and I would like to know is there an additional way to optimize it so it works for large x values too (it already works great for small x values and really big N values). Here is the code:

long double calcMaclaurinPolynom(double x, int N){
    long double result = 0;
    long double atzeretCounter = 2;
    int sign = 1;
    long double fraction = x;

    for (int i = 0; i <= N; i++)
    {
        result += sign*fraction;
        sign = sign*(-1);
        fraction = fraction*((x*x) / ((atzeretCounter)*(atzeretCounter + 1)));
        atzeretCounter += 2;

    }
    return result;
}

Upvotes: 2

Views: 2544

Answers (2)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

You can avoid the sign variable by incorporating it into the fraction update as in (-x*x).

With your algorithm you do not have problems with integer overflow in the factorials.

As soon as x*x < (2*k)*(2*k+1) the error - assuming exact evaluation - is bounded by abs(fraction), i.e., the size of the next term in the series.

For large x the biggest source for errors is truncation resp. floating point errors that are magnified via cancellation of the terms of the alternating series. For k about x/2 the terms around the k-th term have the biggest size and have to be offset by other big terms.

Halving-and-Squaring

One easy method to deal with large x without using the value of pi is to employ the trigonometric theorems where

sin(2*x)=2*sin(x)*cos(x)
cos(2*x)=2*cos(x)^2-1=cos(x)^2-sin(x)^2

and first reduce x by halving, simultaneously evaluating the Maclaurin series for sin(x/2^n) and cos(x/2^n) and then employ trigonometric squaring (literal squaring as complex numbers cos(x)+i*sin(x)) to recover the values for the original argument.

cos(x/2^(n-1)) = cos(x/2^n)^2-sin(x/2^n)^2
sin(x/2^(n-1)) = 2*sin(x/2^n)*cos(x/2^n)

then

cos(x/2^(n-2)) = cos(x/2^(n-1))^2-sin(x/2^(n-1))^2
sin(x/2^(n-2)) = 2*sin(x/2^(n-1))*cos(x/2^(n-1))

etc.

See https://stackoverflow.com/a/22791396/3088138 for the simultaneous computation of sin and cos values, then encapsulate it with

def CosSinForLargerX(x,n):
    k=0
    while abs(x)>1:
        k+=1; x/=2
    c,s = getCosSin(x,n)
    r2=0
    for i in range(k):
        s2=s*s; c2=c*c; r2=s2+c2
        s = 2*c*s
        c = c2-s2
    return c/r2,s/r2

Upvotes: 1

chux
chux

Reputation: 153507

The major issue is using the series outside its range where it well converges.

As OP said "converted x to radX = (x*PI)/180" indicates the OP is starting with degrees rather than radians, the OP is in luck. The first step in finding my_sin(x) is range reduction. When starting with degrees, the reduction is exact. So reduce the range before converting to radians.

long double calcMaclaurinPolynom(double x /* degrees */, int N){
  // Reduce to range -360 to 360
  // This reduction is exact, no round-off error
  x = fmod(x, 360);  

  // Reduce to range -180 to 180
  if (x >= 180) {
    x -= 180;
    x = -x;
  } else if (x <= -180) {
    x += 180;
    x = -x; 
  }

  // Reduce to range -90 to 90
  if (x >= 90) {
    x = 180 - x;
  } else if (x <= -90) {
    x = -180 - x;
  }

  //now convert to radians.
  x = x*PI/180;
  // continue with regular code

Alternative, if using C11, use remquo(). Search SO for sample code.

As @user3386109 commented above, no need to "convert back to degrees".

[Edit]

With typical summation series, summing the least significant terms first improves the precision of the answer. With OP's code this can be done with 

for (int i = N; i >= 0; i--)

Alternatively, rather than iterating a fixed number of times, loop until the term has no significance to the sum. The following uses recursion to sum the least significant terms first. With range reduction in the -90 to 90 range, the number of iterations is not excessive.

static double sin_d_helper(double term, double xx, unsigned i) {
  if (1.0 + term == 1.0)
    return term;
  return term - sin_d_helper(term * xx / ((i + 1) * (i + 2)), xx, i + 2);
}


#include <math.h>
double sin_d(double x_degrees) {

  // range reduction and d --> r conversion from above
  double x_radians = ...

  return x_radians * sin_d_helper(1.0, x_radians * x_radians, 1);
}

Upvotes: 2

Related Questions