No idea for limit the sin(x) formula loop

Question

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

Referancing by the formula;

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;
}

You can set a convergence epsilon value, stopping when there is no appreciable difference to the series sum. Also, you do not need to compute either powers or factorials, as [this answer](https://stackoverflow.com/a/38265293/4142924) shows. — Weather Vane, Oct 01 '17 at 13:52
If you want integer result(-1.0, 0.0, 1.0) , it can be simplified more. — BLUEPIXY, Oct 01 '17 at 13:58
I had implemented this function. To find the answer of sin() for given angle. What I did was i ran the loop for 8 times. After that the value of fact() grew extremely large therefore making the factor `pow(x,2*k+1)/fact(2*k+1)` very small. I recommend running the loop 8 times rather than complication the loop with some epsilon values. — Sreeram TP, Oct 01 '17 at 14:01
@frogwine I also assume that the argument(`int x`) isn't radian, right? — BLUEPIXY, Oct 01 '17 at 14:04
NB: if performance is an issue, note that your sign term, factorial terms and x^n terms can be _progressively_ calculated as you iterate. — Alnitak, Oct 01 '17 at 14:13
@BLUEPIXY yeah, is not radion example; sin(90)'s return value must be 1. — , Oct 01 '17 at 14:14
@frogwine The formula for series expression is based on Radian. So if you want to use that formula You need to convert it to Radian first. If you only want to get integer result, you can simply simplify it with remainder. — BLUEPIXY, Oct 01 '17 at 14:19
My initial thinking was to use an epsilon count too, but on trying it out I also now think it's simpler to just use a fixed iteration count. — Alnitak, Oct 01 '17 at 14:42
BLUEPIXY and @Sreeram TP, you are the mans! Thank you for you two (and of course other helpers) I runned my loop for 8 times and send parameter to sin function as radian, it's worked! Calculating sin30 = 0.5, sin90 = 1.0 thank you for all! — , Oct 01 '17 at 14:54

gsamaras · Answer 1 · 2017-10-01T14:40:07.447

2

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 +π.

edited Oct 01 '17 at 14:40

answered Oct 01 '17 at 13:52

gsamaras

71,951
46
188
305

1

Not only that, it should be of type `double` – meowgoesthedog Oct 01 '17 at 13:57
@gsamaras tp var is sum variable. I'm summing on it. – Oct 01 '17 at 13:58
1

tp should be double as the value computed by the formula will be a float or double. Also the return type of the function is double. – Sreeram TP Oct 01 '17 at 14:07
1

NB: not only must the value `x` be in radians, it must also typically be in the range -π to +π – Alnitak Oct 01 '17 at 14:34

score 0 · Answer 2 · answered Oct 01 '17 at 13:50

0

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

answered Oct 01 '17 at 13:50

Mitch

3,342
2
21
31

i did something

double sin(int x)
{
    double eps=8.85*pow(10,-12);
    double tp=0;
    double prev_ans;
    double cur_ans;

    for(int k=1;;k++)
    {
        cur_ans+=pow(-1,k)*pow(x,2*k+1)/fact(2*k+1);
        if(abs(prev_ans - cur_ans) < eps) break;
        prev_ans=cur_ans;
    }
    return cur_ans;
}

– Oct 01 '17 at 14:18

@frogwine try changing it so that sin takes a double, also make sure that you initialize `prev_ans` to zero. Also I did make a mistake, you should be using `fabs()` as pointed out in the other answer as you are comparing doubles – Mitch Oct 01 '17 at 14:49
`abs` is an integer function, it needs `fabs`. – Weather Vane Oct 01 '17 at 17:46

score 0 · Answer 3 · answered Oct 01 '17 at 15:02

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.

score 0 · Answer 4 · answered Oct 01 '17 at 16:16

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:

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)

score 0 · Accepted Answer · answered Jul 24 '18 at 13:24

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;
}

No idea for limit the sin(x) formula loop

5 Answers5