Arcsin Calculation Algorithm

An Algorithm for Calculating the Inverse Sine of a Given Value
Using Newton’s Method Iteration

Leo T. Meire

Because the Pbasic program does not support an inverse sine function, the calculation of the inverse sine of a given number a, where -1 <= a <= 1, requires development of an algorithm that can be translated into a Pbasic function.

Consider that the problem can be stated : Find the value of X such that X = arcsin(a).

Restate the problem as sin(X) = a (by taking the sine of both sides), and define a function f(X) = sin(X) - a = 0.

The problem is now to determine the value of X that makes this statement true.

Newton’s Method is an iterative formula for numeric determination of the “zeroes” of a differentiable function f(X), starting from an initial value X0 sufficiently close to the true value of X :

X1 = X0 - ( f(X0) / f ’ (X0) ), or, in general,

Xn+1 = Xn - ( f (Xn) / f ’ (Xn) ).

In this case we have defined f(X) = sin (X) - a, and so we can calculate

f ’ (X) = cos (X).

The iteration formula then becomes

X n+1 = Xn - ( (sin (Xn) - a) / cos(Xn) ).

The remaining part of the problem to be determined is the value of the “initial guess” X0. Because the sine function is “well behaved” - it does not have asymptotes and imaginary roots - a good initial value is X0 = a. In fact, for small values of a, this is very close to the true value of X. Both x and a must be expressed in radians; the arcsin function takes values between -pi/2 and pi/2 only (-90° to 90°).

Here is an example of a function written in Pascal to calculate the arcsin of a given value.

function arcsin(a:real):real;
{MUST have -1 < a < 1}
{Returns x = arcsin(a), using Newton’s Method}
var x:real;
const epsilon = 1.0E-7;
begin {arcsin}
x := a;
repeat x := x - (sin(x)-a)/cos(x) until (abs(sin(x)-a) <= epsilon);
arcsin := x;
end; {arcsin}

Again, note that the value of x returned is in radians,
-pi/2 <= x <= pi/2.

Here I have used a constant, epsilon, with a small value (1E-7) to determine when the value of X is sufficiently close to the true value that the iteration may be stopped. This value can be adjusted to converge to a sufficiently accurate value in a reasonable time.

The last remaining problem I can see in the translation of this algorithm to Pbasic is that language’s use of “two’s complement” arithmetic for negative values. Frankly, I have not given sufficient thought to the ramifications of this to be able to state in advance how it will affect the calculations. The easiest way to settle the matter is to code it up and run it. If it works for values of a near -1, 0, and 1, it is probably OK. I am not a Pbasic programmer, so I can’t say how best to do the translation, but I know that there are plenty of clever FIRST people out there. Let me know how it goes. Admittedly, using this method is somewhat like climbing down the chimney because the doors are locked, but at least it will get you in the house.

After looking up some information on ‘two’s complement’ arithmetic, I think the arcsin algorithm should work just fine in Pbasic. You run into problems doing ‘two’s complement’ calculations with very large numbers, and get strange results like adding two large positive numbers and getting a negative result. However, all of the functions in the arcsin calculation are bounded. The sine and cosine range between -1 and 1, the arcsin result will be between -pi/2 and pi/2, the epsilon constant will be a small positive number, and the ‘first guess’ will also fall between -1 and 1. So, this algorithm should rapidly converge to the angle you are looking for and not head off into deep space.

By the way, does anyone find any of this useful, or even interesting? I love this kind of stuff, but then, I have gotten strange looks all my life and am used to them. Actually, I like them. But, have I gone too deep into the Nerdly Woods to interest even confirmed FIRST-ers?

TOO MUCH MATH TO TAKE IN AT ONCE!!! OVERLOAD!!!

J/k, only us “FIRST Doirks” can read that. (J/K too!)

Sorry, guys - I guess I go off the deep end every so often.

(But it really is quite simple in concept…)

It really ain’t to hard to do, but I’m not sure if there’s a sine and cosine function in Pbasic. Anyone knows for sure?

There is definitely a sin and cos in Pbasic.

As to the question of going off the deep end, it is not so much that as the fact that TIME is too short to actually implement such algorithms.

By TIME, I mean 2 things:

  1. the time in the loop to get everything done before you miss too many data packets and get shut down by the Master CPU

  2. the time it takes to program a robot is such that it is barely running in time as it is without the complexities of series expansions to approximate arc anything.

I have always been able to make an approximation of a much lower order of magnitude work for me.

Besides, giving the operators time to practice is more important than anything I can program or design into the robot.

So speaks Dr. Joe :wink:

Joe J.

NOW we’re getting a discussion going!

I agree completely with Dr. Joe - TIME is of the essence in these competitions. Any engineering project depends on control of all four of these aspects - methods, materials, money, and people. With time = money (or resources - CPU time, cash, space, etc.), it is impossible to succeed if too much of a limited resource is expended in non-productive ways.

That being said, maybe it would be an interesting intellectual exercise (it being October, and the new game being two months away) to look at whether this algorithm CAN be implemented, and HOW it could be done. Then each team can decide on whether it is worth their while to do it?

First, let me point out that Newton’s Method does not employ any series expansions - that would be out of the question for arcsin calculations!

I did a quick test to see how many loops it would take for the algorithm to converge at different values of the sine from -1 to +1.

At +/- 1, four loops got to within 0.001 of the exact value of the angle, and ten loops to within 1E-6.
At +/- 0.7, two loops converged within 2E-5.
Near zero, the algorithm converges to very nearly the exact value after the first loop (as one would expect).

A value within 0.001 is going to be close enough in any competition - maybe too close. What mechanism is anyone going to devise that needs more accuracy than that? You could implement a version that just goes through twice, calls it done, and still gets very close to the mark each time. I don’t know what is involved with the lookup function in terms of CPU usage, and how it would compare to a two-loop calculation.

Now for the real fun part - how Pbasic implements the sine and cosine functions. According to the BASIC Stamp Programming Manual 2.0c, the SIN and COS operators use a method that is a little unusual - at least to me. These operators take an argument not in degrees or radians, but in “binary radians”, sometimes known as “brads”. That is, they break the circle into 0 to 255 units (brads), each equivalent to 1.406 degrees. So, the argument passed to the SIN function is an angle in brads, expressed as an integer from 0 to 255.

The value returned by the function is based not on a unit circle like the ordinary sine value in our calculators, but on a circle of radius 127. We see that number a lot in this language, don’t we? The results are given in two’s complement form to accommodate negative values. The functions return integer values from -127 to +127.

So, at zero degrees we get sin(0 brads) = 0; at 45 degrees we have sin(32 brads) = 90. At 90 degrees, sin(64 brads) = 127; at 180 degrees, sin(128 brads) = 0; and at 270 degrees, sin (192 brads) = -127.

Brads X (180/128) = degrees & degrees X (128/180) = brads.

Now, whether it is worth implementing this method is something else again. It looks like it will work, but that is only 1/4 of the overall engineering problem. Does your team have the people necessary to do it, the time to get it done, the design that can benefit from it? I don’t know - you tell me.

But, don’t you think the mental exercise is beneficial? There’s a lot to this problem, and a successful implementation would be, if nothing else, satisfying to the person who actually accomplished it. That’s part of the FIRST experience too, is it not?

This is only my very humble opinion. If you actually want to build a robot that gets out on the field in time, follow Dr. Joe’s advice. But, if you are looking for some worthy intellectual opponent during the off-season, give this one a try…

Again, let me know how it goes, and good luck.

Is there an ArcTan function in PBasic?

If there is, the problem of finding ArcSin and ArcCos can be calculated using ArcTan by using:

ArcSin(x) = ArcTan (x / sqrt(1 - sqr (X)))
ArcCos(x) = ArcTan (sqrt(1 - sqr (X)) / x)

Of course, this assumes that there is a sqrt(x) function in PBasic as well.

Anthony.

The real test would be to truncate all your calculations to 16 bits and then to count the number of interations needed to get to within 1 bit of byte result (1 out of 256).

We only have 8 bits of resolution in the inputs. We have LESS than 8 bits of resolution in the outputs (I don’t recall the actual numbers but I THINK it is something like 64 PWM steps in each direction – not more). Finally, the sin and cos functions we need in our newton’s method are only 8 bit functions.

There may be some case when more resolution can be utilized but I are skeptical.

If you got the answer to 8 significant bits, I suppose that would be plenty.

Joe J.

I did a quick search on google and found that arcsin(x) = pi/2 - sqrt(1 - x)(a0 + a1x + a2x^2 + a3*x^3),

where
a0 = 1.5707288
a1 = -0.2121144
a2 = 0.0742610
a3 = -0.0187293

============
I tested this for x = 1/2:
Expected result: pi/6 = 0.523599 rads
Result using above equation: 0.523645 rads
a difference of 0.000046 rads (0.0026 degrees)

and again for x = sqrt(3)/2:
Expected result: pi/3 = 1.047198 rads
Result using above equation: 1.047174 rads
a difference of 0.000024 rads (0.0014 degrees)

then I rounded all of the a values to the nearest thousanth
a0 = 1.571
a1 = -0.212
a2 = 0.074
a3 = -0.019

and did it again

============
x = 1/2:
Expected result: pi/6 = 0.523599 rads
Result using above equation: 0.523483 rads
a difference of 0.000116 rads (0.0066 degrees)

x = sqrt(3)/2:
Expected result: pi/3 = 1.047198 rads
Result using above equation: 1.047174 rads
a difference of 0.000023 rads (0.0013 degrees)

Even the one rounded to the thousandths seems accurate enough for me. However, the equation still seems fairly proccessor intensive. I’m not sure how long it would take to calculate it using this method. Maybe someone could test both this method and the other method to see which one is faster.

BTW, you might even be able to round to the hundredths. I’m too lazy to do that right now :smiley:

The Parallax Manual only mentions SIN, COS, and SQR (page 236 of Manual 1.9). Just imagine the fun of finding TAN (= SIN/COS), let alone computing ArcTAN.

It might be more practical to use a lookup table, and leave large spaces in tabular places you don’t need to “calculate”.

Using an ArcTan calculation, with the additional problem of noticing and debugging an incorrect sign or mistyped number used in a formula, could lead to many late nights.

In going through the Basic Stamp Pbasic manual, the greatest obstacle I can see to implementing the Newton’s Method algorithm is in the division step. The manual recommends that division NOT be performed in cases where either numerator or denominator are negative; results can be unpredictable. For the demoninator, cos(x), this is not a problem over the restricted range -pi/2 to +pi/2. However, the numerator can certainly go negative. The additional steps necessary to handle that make the algorithm needlessly complicated and slow - getting back to Dr. Joe’s wise advice about time of execution.

The second difficulty comes from the different sizes of the sin function argument and result. SIN takes a byte argument (-127 to +127), and outputs a word result (-32,767 to +32,767). Scaling the word result back down to byte-size to send back through the function at the next iteration could be tricky (i.e., an infinite sink of time trying to work out the programming).

It could be that Jay Lundy’s method (which looks to me like a truncated Maclaurin series, or Taylor series about x=0) may actually be more efficient.

Math on a Basic Stamp is definitely not for the faint of heart.