CODEWITHSUNDEEP
algorithmmathapproximation
user8704018
user8704018

Reputation:

No idea for limit the sin(x) formula loop

Making c/++ function to compute sin(x)

Referancing by the formula;
http://wiki.ubc.ca/images/math/4/4/7/447a79826774707026bbefcd76962d3a.png

But i don't know where to stop infinity sum, and i don't want libraries for quick-answer. I want to make it. I made factorial function well. (fact() function)

double sin(int x) {
    int tp=0;
    for(int k=1; ;k++)
    {
        tp+=pow(-1,k)*pow(x,2*k+1)/fact(2*k+1);
    }
    return tp;
}

Upvotes: 2

Views: 304

Answers (5)

user8704018
user8704018

Reputation:

I ran the loop by 8 times, it's approaching to true value if you improve loop timing. Thanks @Sreeram TP

double sin(int x) {
    int tp=0;
    for(int k=1;k<8 ;k++)
    {
        tp+=pow(-1,k)*pow(x,2*k+1)/fact(2*k+1);
    }
    return tp;
}

Upvotes: 0

meowgoesthedog
meowgoesthedog

Reputation: 15045

Before we go any further, you must use double for return types and accumulation, otherwise the result will be (heavily) truncated.

You should stop when adding the change term does not change the current accumulated summation value, i.e. there is underflow. This occurs when:

The magnitude of the change term is less than the value of the lowest binary digit of the accumulator.

To check this, it is useful to extract the exponent of the floating point number. This SO post provides a useful function: frexp.

We can then check that the exponent of the change term is smaller than the exponent of the accumulator minus the number of bits in the mantissa of the double type (52 bits). If you need to know more, search for IEEE floating point format.

Another point to improve numerical accuracy and speed: let's examine the relationship between successive terms:

enter image description here

So instead of calculating each term from scratch (using the power and factorial functions), we can calculate each term iteratively as we go.

Sample code:

// best to use an existing C-library define for this, if exists
#define DBL_MANTISSA_BIT 52

double mysin_new(double x)
{
   double term = x;
   double sum = term;

   for (int n = 1; ; n++) {
      // next term in the series
      term *= - (x / (2*n)) * (x / (2*n+1));

      // check if the exponent is small enough for us to stop
      int e_term, e_sum;
      frexp(term, &e_term);
      frexp(sum, &e_sum);
      if (e_term <= e_sum - DBL_MANTISSA_BIT) break;

      sum += term;
   }

   return sum;
}

Some numerical tests:

angle = 0.157080
  library  : 0.15643446504023086896f, 0011111111000100000001100000101101100111101010000101001101110101b
  mysin_old: 0.15643446504023086896f, 0011111111000100000001100000101101100111101010000101001101110101b
  mysin_new: 0.15643446504023089672f, 0011111111000100000001100000101101100111101010000101001101110110b

angle = 0.314159
  library  : 0.30901699437494739575f, 0011111111010011110001101110111100110111001011111110100101001111b
  mysin_old: 0.30901699437494745126f, 0011111111010011110001101110111100110111001011111110100101010000b
  mysin_new: 0.30901699437494739575f, 0011111111010011110001101110111100110111001011111110100101001111b

angle = 0.471239
  library  : 0.45399049973954674897f, 0011111111011101000011100010111000101011010001001101111000000000b
  mysin_old: 0.45399049973954674897f, 0011111111011101000011100010111000101011010001001101111000000000b
  mysin_new: 0.45399049973954674897f, 0011111111011101000011100010111000101011010001001101111000000000b

angle = 0.628319
  library  : 0.58778525229247313710f, 0011111111100010110011110010001100000100011101010101101001011110b
  mysin_old: 0.58778525229247313710f, 0011111111100010110011110010001100000100011101010101101001011110b
  mysin_new: 0.58778525229247313710f, 0011111111100010110011110010001100000100011101010101101001011110b

angle = 0.785398
  library  : 0.70710678118654746172f, 0011111111100110101000001001111001100110011111110011101111001100b
  mysin_old: 0.70710678118654746172f, 0011111111100110101000001001111001100110011111110011101111001100b
  mysin_new: 0.70710678118654746172f, 0011111111100110101000001001111001100110011111110011101111001100b

angle = 0.942478
  library  : 0.80901699437494745126f, 0011111111101001111000110111011110011011100101111111010010101000b
  mysin_old: 0.80901699437494734024f, 0011111111101001111000110111011110011011100101111111010010100111b
  mysin_new: 0.80901699437494734024f, 0011111111101001111000110111011110011011100101111111010010100111b

angle = 1.099557
  library  : 0.89100652418836778779f, 0011111111101100100000110010000000011101001111010010110001101100b
  mysin_old: 0.89100652418836767676f, 0011111111101100100000110010000000011101001111010010110001101011b
  mysin_new: 0.89100652418836767676f, 0011111111101100100000110010000000011101001111010010110001101011b

angle = 1.256637
  library  : 0.95105651629515353118f, 0011111111101110011011110000111000010011010001000101010011111111b
  mysin_old: 0.95105651629515353118f, 0011111111101110011011110000111000010011010001000101010011111111b
  mysin_new: 0.95105651629515353118f, 0011111111101110011011110000111000010011010001000101010011111111b

angle = 1.413717
  library  : 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
  mysin_old: 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
  mysin_new: 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b

angle = 1.570796
  library  : 1.00000000000000000000f, 0011111111110000000000000000000000000000000000000000000000000000b
  mysin_old: 1.00000000000000022204f, 0011111111110000000000000000000000000000000000000000000000000001b
  mysin_new: 1.00000000000000000000f, 0011111111110000000000000000000000000000000000000000000000000000b

angle = 1.727876
  library  : 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
  mysin_old: 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
  mysin_new: 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b

angle = 1.884956
  library  : 0.95105651629515364220f, 0011111111101110011011110000111000010011010001000101010100000000b
  mysin_old: 0.95105651629515364220f, 0011111111101110011011110000111000010011010001000101010100000000b
  mysin_new: 0.95105651629515375323f, 0011111111101110011011110000111000010011010001000101010100000001b

angle = 2.042035
  library  : 0.89100652418836789881f, 0011111111101100100000110010000000011101001111010010110001101101b
  mysin_old: 0.89100652418836800983f, 0011111111101100100000110010000000011101001111010010110001101110b
  mysin_new: 0.89100652418836800983f, 0011111111101100100000110010000000011101001111010010110001101110b

angle = 2.199115
  library  : 0.80901699437494745126f, 0011111111101001111000110111011110011011100101111111010010101000b
  mysin_old: 0.80901699437494756229f, 0011111111101001111000110111011110011011100101111111010010101001b
  mysin_new: 0.80901699437494756229f, 0011111111101001111000110111011110011011100101111111010010101001b

angle = 2.356194
  library  : 0.70710678118654757274f, 0011111111100110101000001001111001100110011111110011101111001101b
  mysin_old: 0.70710678118654768376f, 0011111111100110101000001001111001100110011111110011101111001110b
  mysin_new: 0.70710678118654757274f, 0011111111100110101000001001111001100110011111110011101111001101b

angle = 2.513274
  library  : 0.58778525229247324813f, 0011111111100010110011110010001100000100011101010101101001011111b
  mysin_old: 0.58778525229247324813f, 0011111111100010110011110010001100000100011101010101101001011111b
  mysin_new: 0.58778525229247324813f, 0011111111100010110011110010001100000100011101010101101001011111b

angle = 2.670354
  library  : 0.45399049973954685999f, 0011111111011101000011100010111000101011010001001101111000000010b
  mysin_old: 0.45399049973954685999f, 0011111111011101000011100010111000101011010001001101111000000010b
  mysin_new: 0.45399049973954691550f, 0011111111011101000011100010111000101011010001001101111000000011b

angle = 2.827433
  library  : 0.30901699437494750677f, 0011111111010011110001101110111100110111001011111110100101010001b
  mysin_old: 0.30901699437494745126f, 0011111111010011110001101110111100110111001011111110100101010000b
  mysin_new: 0.30901699437494745126f, 0011111111010011110001101110111100110111001011111110100101010000b

angle = 2.984513
  library  : 0.15643446504023097998f, 0011111111000100000001100000101101100111101010000101001101111001b
  mysin_old: 0.15643446504023095223f, 0011111111000100000001100000101101100111101010000101001101111000b
  mysin_new: 0.15643446504023095223f, 0011111111000100000001100000101101100111101010000101001101111000b

angle = 3.141593
  library  : 0.00000000000000012246f, 0011110010100001101001100000000000000000000000000000000000000000b
  mysin_old: 0.00000000000000023566f, 0011110010110000111110110010110111010110111000110001010101001000b
  mysin_new: 0.00000000000000023566f, 0011110010110000111110110010110111010110111000110001010101001001b

The new adaptive method wins out some of the times compared to the original, draws most of them, and loses a number of others. (Note it is assumed that the library method sin is perfectly accurate upto the precision of double)

Upvotes: 0

Alnitak
Alnitak

Reputation: 340055

Posted for reference, here's my version that uses progressive calculation of the various components of the series expansion, instead of calling out to an external O(n) factorial function:

double mysin(double x)
{
    int sign = 1;
    int iteration = 1;
    double xn = x;
    double factorial = 1.0;
    double result = 0.0;
    double delta;

    while (iteration < 31) {
        delta = xn / factorial;
        result += sign * delta;
        xn *= x * x;
        sign = -sign;
        factorial *= ++iteration;
        factorial *= ++iteration;
    }
    return result;
}

From my tests there doesn't appear to be any significant advantage in accuracy (or performance) by jumping out early if fabs(delta) becomes too small, unless the initial value of x is also very small.

Upvotes: 0

gsamaras
gsamaras

Reputation: 73444

You should stop when the value computed (tp) is not changed much. Check for the difference of the newly computed value with the old one, and if the absolute difference (use fabs() for that) is less than epsilon (you define that), stop iterating.

When you stop, you will end up with an approximation of the result of sin.

Moreover, int tp=0; doesn't seem correct to me. You should initialize it with the formula result, when k is 0. And more importantly, it should be of type double, since int is for integers.

Furthermore, not only must the value x be in radians, it must also typically be in the range -π to +π.

Upvotes: 2

Mitch
Mitch

Reputation: 3416

What you need is an epsilon value, a value that you say when I am more accurate than that number we want to stop. It would look something like if(abs(prev_ans - cur_ans) < epsilon) break;

In math words the n and n+1 partial sum of the series vary by so little that continuing the series wouldn't produce any more desired accuracy

Upvotes: 0

Related Questions