Threaded Mode | Linear Mode

Valentin Albillo · 04-05-2019, 01:58 AM

Hi, all:

At long last, today it's time for my final original solution, namely:

Tier 5 - The Challenge:

Consider the function cin(x) which has the defining property that   cin(cin(cin(x))) = sin(x).

Write a program or function which accepts an argument x in the range [-Pi, Pi] and outputs the corresponding value of cin(x) correct to at least 8-10 digits in the whole range. Use it to tabulate cin(x) for x = 0.0, 0.2, 0.4, ..., 1.0 and also to compute cin(Pi/2), cin(-0.71), cin(2.019).

My original solution:

My original solution for the HP-71B is the following user-defined function (plus initialization code):

       1 DESTROY ALL @ OPTION BASE 1 @ DIM C(7) @ READ C
       2 DATA 1,-1/18,-7/1080,-643/408240,-13583/29393280,-29957/215550720,-24277937/648499737600

       3 DEF FNC(X) @ L=0 @ M=1/3 @ REPEAT @ X=SIN(X) @ L=L+1 @ UNTIL ABS(X)<M
       4 S=0 @ FOR Z=1 TO 7 @ S=S+C(Z)*X^(2*Z-1) @ NEXT Z
       5 FOR Z=1 TO L @ S=ASIN(S) @ NEXT Z @ FNC=S @ END DEF

Instead of tabulating it for 0.0, 0.2, ..., 1.0 as I originally asked, let's better tabulate it for x from 0 to Pi/2 in steps of Pi/10:

       6 FOR X=0 TO PI/2 STEP PI/10
       7 Y=FNC(FNC(FNC(X))) @ DISP X;FNC(X);Y;SIN(X);Y-SIN(X) @ NEXT X

       >FIX 10
       >RUN
                   x           cin(x)        cin(cin(cin(x)))    sin(x)       Error
              ----------------------------------------------------------------------
              0.0000000000   0.0000000000      0.0000000000   0.0000000000     0
              0.3141592654   0.3124163699      0.3090169944   0.3090169944 -1.0E-12
              0.6283185307   0.6138343796      0.5877852523   0.5877852523   2.2E-11
              0.9424777961   0.8897456012      0.8090169944   0.8090169944   4.1E-11
              1.2566370614   1.1122980783      0.9510565164   0.9510565163   1.0E-10
              1.5707963268   1.2103683445      1.0000000000   1.0000000000   1.0E-11

So we've got 10 correct decimals or better, as the error in cin(x) is even smaller than the error in cin(cin(cin(x)))-sin(x) which doesn't exceed 10^-10. As for the discrete values asked in the challenge:

       >FIX 10 @ FNC(PI/2); FNC(-0.71); FNC(2.019)

              1.2103683445   -0.6887785253   1.0269233188

Notes:

Line 4 evaluates the formal series in a simple loop but that's not optimal. I know of several better ways to evaluate the series but I don't want to digress from the main subject, which is the computation of cin(x).
Albert Chan found the correct conjugation (sin/arcsin) procedure to increase accuracy and almost duplicated my original solution but there's an important difference which affects both accuracy and running time. He used up to the x⁷ term in his formal series expansion:

       x - (1/18) x³ - (7/1080) x⁵ - 0.00158 x⁷

and then iterated the sine of the argument till it got < 0.1, while my original solution uses up to the x¹³ term:

       x - 1/18 x³ - 7/1080 x⁵ - 643/408240 x⁷ - 13583/29393280 x⁹ - 29957/215550720 x¹¹ - 24277937/648499737600 x¹³

and iterates until the sin gets < 1/3. This way significantly fewer sin/arcsin iterations are needed and the computation is both more accurate and faster. For instance, to see how many iterations my function performs when computing cin(Pi/2) just execute this:

       >FNC(PI/2);L
                            1.2103683445     24

thus 24 sin/arcsin were necessary for this argument while in J-F Garnier's HP-71B version of Albert Chan's code several hundred sines/arcsines are necessary to bring this argument below 0.1, which explains why it takes much longer and worse, several decimal places are lost in the process.
My solution will also work for tin(x) , defined as tin(tin(x)) = sin(x) , by simply replacing the coefficients in the DATA statement at line 2 by those of its own formal series, namely:

x - x³/12 - x⁵/160 - 53/40320 x⁷ - 23/71680 x⁹ - 92713/1277337600 x¹¹ + ...

and of course it will also work for any other such functions as well.
The coefficients of the formal series for cin(x) and tin(x) can be obtained in a number of ways (even manually for the first 4 or so !), most easily by using some CAS which can deal with formal series (even a version of Newton's method for solving f(x) = 0 can be put to the task), but it's important to be aware that both formal series do not converge.

In fact, their radius of convergence is 0 and thus they behave like asymptotic series, so you can't get arbitrarily accurate results by taking more and more terms, you must instead truncate the series after a certain number of terms to get the most accurate results. Using further terms only worsens the accuracy.
Although at first sight the coefficients of the formal series for cin(x) and tin(x) seem to (slowly) get smaller and smaller, matter of fact they tend to grow ever bigger after a while, tending to infinity. For instance, for tin(x) we find that the smallest coefficient in absolute value is:

       Coeff₃₇ = -0.000000000594338574503

but afterwards we have, e.g.:

       Coeff₁₀₁ = 0.0833756228055

       Coeff₁₅₁ = 388536047335.239

       Coeff₂₀₁ = 6555423874651256623811186991.51

       Coeff₂₅₁ = -35365220492708296140377087748804440170254492009.57

That's all for Tier 5, I could say a whole lot more about this topic and post additional code and results but this post is long enough as it is so I'll stop right now.

Thank you very much to Albert Chan, J-F Garnier, Oulan and Gerson W. Barbosa for your valuable contributions and to Werner for your interest, I hope you enjoyed it all ! Smile

V.