+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Arcsin in Valentin Albillo's HP-12C Tried & Tricky Trigonometrics (/thread-14048.html)
Arcsin in Valentin Albillo's HP-12C Tried & Tricky Trigonometrics - jklsadf - 11-25-2019 06:40 AM
I've been going through Valentin Albillo's "HP-12C Tried & Tricky Trigonometrics"
article trying to understand the algorithms used, and am so far completely stumped as to how the arcsin calculation is done. http://www.hpcc.org/datafile/hp12/12c_TrigonometryFunctions.pdf
The sin calculation uses a fairly straightforward application of the taylor series expansion of sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ... = sum(-1^n * x^(2n+1) / (2n+1)! from 0 to infinity) as stated in the article, with register 5 holding the running sum, and register 6 holding n.
For arcsin, the article states that "To guarantee fast convergence for all x, a suitable scaling of the input value is previously performed, followed by a change of variable, and finally a Taylor Series Expansion is used to generate an intermediate result, which is then scaled back to give the final result." The motivation of guaranteeing fast convergence makes sense, as does calculating arctan in terms of arcsin (since all of the terms in the arcsin series have the same sign) instead of the other way around, but I don't fully understand what the scaling or change of variable is.
As far as I can tell, register 5 (which holds the running sum), first has the scaling (1 - sqrt(1 - x^2)) / 2 done repeatedly (3 times), which seems like it might be trying to convert to arccos or something (at least the sqrt(1 - x^2) part)? Or maybe it's related to the fact that arcsin(x) = integral(dt/sqrt(1 - t^2) from 0 to x)? There's some additional operations before getting the initial values for register 5 (which hold the running sum) and register 6 (which holds n), which I understand even less. Line 59 seems to always compute 0!, which is 1.
Then after the change of variable, it looks like the series expansion used is simply: x - x^3/3 * x^5/5 - x^7/7 + ... which is surprisingly simple compared to the original series expansion for arcsin(x).
The final scaling back to give the final result doesn't make much sense either, and seems to consist of just multiplying by 8 and correcting the sign. The only thing I can think of is 2^3 is 8, and 3 was the number of scaling iterations done initially, maybe from a half angle formula. Miraculously this all seems to work and give the correct answer, but I have no idea how.
Maybe this will make more sense if I look at it more, but I was hoping Valentin (if he reads this) or someone else could shed some more light on how exactly arcsin(x) is being calculated, and why. The last time I've had to look at anything similar was in school, and although I'm probably a bit younger than most people on this forum, that was still a long time ago.
RE: Arcsin in Valentin Albillo's HP-12C Tried & Tricky Trigonometrics - jklsadf - 11-25-2019 06:56 AM
After a little more googling, the initial scaling does seem to be using the half angle formula somehow sin(x/2) = +/- sqrt((1 - cos(x)) / 2). There was an additional square root on line 52 that gets applied to register 5 on all but the last iteration. The final multiplication by 8 is to undo this scaling.
The taylor series used is apparently that for arctan(x) = x - x^3/3 * x^5/5 - x^7/7 + ... which alternates signs for each term in the series, but I guess convergence is fine after the scaling. There must have been a change of variable to allow using the series for arctan, but I'm not sure how.
RE: Arcsin in Valentin Albillo's HP-12C Tried & Tricky Trigonometrics - Albert Chan - 11-25-2019 04:32 PM
(11-25-2019 06:56 AM)jklsadf Wrote: The taylor series used is apparently that for arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + ... which alternates signs for each term in the series, but I guess convergence is fine after the scaling. There must have been a change of variable to allow using the series for arctan, but I'm not sure how.
θ = asin(x)
After applying half-angle formula 3 times: θ/8 = asin(x') = atan(x' / √(1-x'²))
RE: Arcsin in Valentin Albillo's HP-12C Tried & Tricky Trigonometrics - jklsadf - 11-26-2019 01:02 AM
Yes that's it. Now that I look at it again, lines 61 to 66 are just using the identity asin(x) = atan(x / sqrt(1 - x^2)), where the scaled x^2 is stored in register 5, since the square root on line 52 isn't applied to register 5 in the last iteration, but its value is on the stack for the divide on line 64.
I got confused keeping track of stuff with the instructions in between. Line 59 has to be factorial (or could've been e^x or 10^x) to keep from having an extra 0 on the stack.
As far as I can tell, the only reason for doing the change of variable and using the arctan series instead of arcsin is the series is much simpler and thus fewer program steps. More advanced RPN programmables could've even combined lines 05 through 21 with lines 66 through 81 since they're very similar series, but the 12C doesn't have flags or subroutines.
RE: Arcsin in Valentin Albillo's HP-12C Tried & Tricky Trigonometrics - Don Shepherd - 11-26-2019 03:33 AM
(11-26-2019 01:02 AM)jklsadf Wrote: but the 12C doesn't have flags or subroutines.
Officially, that's true. But see this for an example of a 12c program that has not only flags and subroutines, but indirect addressing as well!
RE: Arcsin in Valentin Albillo's HP-12C Tried & Tricky Trigonometrics - Albert Chan - 11-27-2019 06:31 PM
(11-25-2019 06:56 AM)jklsadf Wrote: initial scaling does seem to be using the half angle formula somehow sin(x/2) = +/- sqrt((1 - cos(x)) / 2).
Above half-angle formula, 1 - cos(x) might introduce big subtraction cancellation errors:
Example, using HP-12C
asin(x) = 2 sign(x) asin(√((1 - √(1-x²)) / 2))
→ asin(0.123) = 2 asin(0.06161708083) = 4 asin(0.03082318608) = 8 asin(0.01541342434)
If there is enough room, this is more accurate, and avoided code to fix sign.
sin(θ/2) = sin(θ) / √(2 + 2 cos(θ)), where -pi/2 ≤ θ ≤ pi/2
asin(x) = 2 asin(x / √(2 + 2 √(1-x²)))
→ asin(0.123) = 2 asin(0.06161708093) = 4 asin(0.03082318602) = 8 asin(0.01541342403)
If there is still room for 2 steps, we can skip 1 asin reduction loop.
asin(x) = atan(x / √(1-x²)) = 2 atan(x / (1 + √(1-x²)))
→ 4 asin(0.03082318602) = 8 atan(0.01541525527)