What wrong with this implement of this arcsine approximate in C#

Question

This is a formula to approximate arcsine(x) using Taylor series from this blog

This is my implementation in C#, I don't know where is the wrong place, the code give wrong result when running: When i = 0, the division will be 1/x. So I assign temp = 1/x at startup. For each iteration, I change "temp" after "i". I use a continual loop until the two next value is very "near" together. When the delta of two next number is very small, I will return the value.

My test case: Input is x =1, so excected arcsin(X) will be arcsin (1) = PI/2 = 1.57079633 rad.

class Arc{
            static double abs(double x)
            {
                return x >= 0 ? x : -x;
            }

            static double pow(double mu, long n)
            {
                double kq = mu;
                for(long i = 2; i<= n; i++)
                {
                    kq *= mu;
                }
                return kq;
            }

            static long  fact(long n)
            {
                long gt = 1;
                for (long i = 2; i <= n; i++) {
                    gt *= i;
                }
                return gt;
            }


            #region arcsin
            static double arcsinX(double x) {
              int i = 0;
              double temp = 0;
              while (true)
               {
               //i++;
               var iFactSquare = fact(i) * fact(i);
               var tempNew = (double)fact(2 * i) / (pow(4, i) * iFactSquare * (2*i+1)) * pow(x, 2 * i + 1) ;
            if (abs(tempNew - temp) < 0.00000001)
            {
                return tempNew;
            }
            temp = tempNew;
            i++;
        }
    }

            public static void Main(){
                Console.WriteLine(arcsin());
                Console.ReadLine();
            }
        }

So what input did you give, what was the result, and what did you expect? Note that `fact(2 * i)` is going to overflow pretty quickly... — Jon Skeet, Jan 21 '16 at 14:49
First suspicious place is `fact` which is declared as `long` (as that's why *integer overflow* is quite possible); make `fact` to return `double` — Dmitry Bychenko, Jan 21 '16 at 14:52
My input is x =1, so excected arcsin(X) will be arcsin (1) = PI/2 = 1.57079633 rad. But it return -1.203949002882805E-11 — Andiana, Jan 21 '16 at 14:52
is that ok that you start not from `0` but from `1`? i mean `int i = 0;` and then `i++;` before calculation? Shouldn't `i++;` be in the end? — teo van kot, Jan 21 '16 at 14:58
another problem is in `var tempNew = fact(2 * i)...` it should be `var tempNew = temp + fact(2 * i) ...` — Dmitry Bychenko, Jan 21 '16 at 15:02
@teovankot, I changed it to start from 0 and i++ at the end. Before, I assign temp = 1/x (equal to temp (0)). But this way is more naturally. — Andiana, Jan 21 '16 at 15:08
@DmitryBychenko , the temp is calculate with every index, The temps is not relate with the previous temp. — Andiana, Jan 21 '16 at 15:08
Today would be a good day to learn how to use your debugger. Step through **every line of the code**, and watch **every time a variable changes**. Before each variable changes, make a mental prediction of what the new value will be. Eventually your prediction will be wrong. Either the code is wrong or your mental model is wrong; either way, you've found a problem that needs fixing. — Eric Lippert, Jan 21 '16 at 15:12

Lutz Lehmann · Accepted Answer · 2016-01-22T08:57:04.693

4

In many series evaluations, it is often convenient to use the quotient between terms to update the term. The quotient here is

                (2n)!*x^(2n+1)       4^(n-1)*((n-1)!)^2*(2n-1)
a[n]/a[n-1] = ------------------- * --------------------- -------
              (4^n*(n!)^2*(2n+1))       (2n-2)!*x^(2n-1)  

  =(2n(2n-1)²x²)/(4n²(2n+1)) 
  = ((2n-1)²x²)/(2n(2n+1))

Thus a loop to compute the series value is

sum = 1;
term = 1;
n=1;
while(1 != 1+term) {
    term *= (n-0.5)*(n-0.5)*x*x/(n*(n+0.5));
    sum += term;
    n += 1;
}
return x*sum;

The convergence is only guaranteed for abs(x)<1, for the evaluation at x=1 you have to employ angle halving, which in general is a good idea to speed up convergence.

edited Jan 22 '16 at 08:57

answered Jan 21 '16 at 15:13

Lutz Lehmann

25,219
2
22
51

2

It turns out this does converge at x=1, but because it is on the boundary of the radius of convergence, the convergence will be painfully slow. – Teepeemm Jan 21 '16 at 17:26
Hi @Lutzl, could you explain to me how do you convert x ^ (2i + 1), logarithmically, is converted to pure multiplication and division? – Andiana Jan 22 '16 at 01:52
1

I use that in the quotient, x^(2n+1)/x^(2(n-1)+1)=x^(2n+1)/x^(2n-1)=x^2. There is no logarithm involved. In `term` the factors x^2 accumulate to give x^(2n) in the n-th step, the return value gets multiplied by x to give the remaining factor to x^(2n+1). – Lutz Lehmann Jan 22 '16 at 08:59

score 3 · Answer 2 · answered Jan 21 '16 at 15:27

You are saving two different temp values (temp and tempNew) to check whether or not continuing computation is irrelevant. This is good, except that you are not saving the sum of these two values.

This is a summation. You need to add every new calculated value to the total. You are only keeping track of the most recently calculated value. You can only ever return the last calculated value of the series. So you will always get an extremely small number as your result. Turn this into a summation and the problem should go away.

score 2 · Answer 3 · edited May 23 '17 at 10:27

NOTE: I've made this a community wiki answer because I was hardly the first person to think of this (just the first to put it down in a comment). If you feel that more needs to be added to make the answer complete, just edit it in!

The general suspicion is that this is down to Integer Overflow, namely one of your values (probably the return of fact() or iFactSquare()) is getting too big for the type you have chosen. It's going to negative because you are using signed types — when it gets to too large a positive number, it loops back into the negative.

Try tracking how large n gets during your calculation, and figure out how big a number it would give you if you ran that number through your fact, pow and iFactSquare functions. If it's bigger than the Maximum long value in 64-bit like we think (assuming you're using 64-bit, it'll be a lot smaller for 32-bit), then try using a double instead.

What wrong with this implement of this arcsine approximate in C#

3 Answers3