Arcsin.cin

From TORI
Jump to: navigation, search

arcsin.cin is the complex double routine in C++ for evaluation of function ArcSin of complex argument.

For some reason, the compiler recognises exp, log, sin, cos, of complex argument, but fails with asin and acos. For me, it is easier to write the efficient implementation, than to search for the appropriate compiler with these routines. T (talk) 14:02, 4 December 2013 (JST)

The routine and the test are copypasted below. The generators of some figures still use the primitive implementation, that provides of order of 14 decimal digits instead of 15 digits. This does not affect the explicit plots nor the complex maps; but the lost of precision is not so good for the numerical testing of routines for iteration of sin, see articles SuSin and AuSin. However, the future figures are supposed to use the precise routine below.

Contents

Routine

//

z_type asin0(z_type z){ int m,M; z_type q,s; DB c[41]={
1.000000000000000, 0.1666666666666667, 0.07500000000000000, 0.04464285714285714, 0.03038194444444444,
0.02237215909090909, 0.01735276442307692, 0.01396484375000000, 0.01155180089613971, 0.009761609529194079,
0.008390335809616815, 0.007312525873598845, 0.006447210311889648, 0.005740037670841923, 0.005153309682319904,
0.004660143486915096, 0.004240907093679363, 0.003880964558837669, 0.003569205393825935, 0.003297059503473485,
0.003057821649258031, 0.002846178401108942, 0.002657870638207290, 0.002489448678246883, 0.002338091892111975,
0.002201473973710138, 0.002077661032518167, 0.001965033616277284, 0.001862226406403127, 0.001768081120515418,
0.001681609393583107, 0.001601963275351444, 0.001528411596122568, 0.001460320894079115, 0.001397139917630253,
0.001338386951275178, 0.001283639387629029, 0.001232525098500017, 0.001184715256162439, 0.001139918330702224,
0.001097875046591447};
M=40; q=z*z; s=c[M]*q; for(m=M-1;m>0;m--) {s+=c[m]; s*=q;}
return z*(1.+s); }

z_type asin1(z_type z){ int m,M; z_type q,qq,s; DB c[41]=
{2.000000000000000, 0.3333333333333333, 0.1500000000000000, 0.08928571428571429, 0.06076388888888889,
0.04474431818181818, 0.03470552884615385, 0.02792968750000000, 0.02310360179227941, 0.01952321905838816,
0.01678067161923363, 0.01462505174719769, 0.01289442062377930, 0.01148007534168385, 0.01030661936463981,
0.009320286973830192, 0.008481814187358726, 0.007761929117675338, 0.007138410787651869, 0.006594119006946969,
0.006115643298516061, 0.005692356802217884, 0.005315741276414580, 0.004978897356493767, 0.004676183784223950,
0.004402947947420276, 0.004155322065036335, 0.003930067232554567, 0.003724452812806255, 0.003536162241030836,
0.003363218787166214, 0.003203926550702888, 0.003056823192245135, 0.002920641788158231, 0.002794279835260507,
0.002676773902550357, 0.002567278775258057, 0.002465050197000034, 0.002369430512324879, 0.002279836661404447,
0.002195750093182894};
M=29; qq=.5-.5*z; q=sqrt(qq); s=c[M]*q; for(m=M-1;m>0;m--) {s+=c[m]; s*=qq;}
return .5*M_PI - q*(2.+s) ; }

z_type asin2(z_type z){ return I*log(-I*z+sqrt(1.-z*z)); }
z_type asin3(z_type z){ return -I*log(I*z+sqrt(1.-z*z)); }

//z_type asin0123(z_type z){
z_type arcsin(z_type z){
if(abs(z)<0.7) return asin0(z);
if(abs(z-1.)<.7) return asin1(z);
if(abs(z+1.)<.7) return -conj(asin1(-conj(z)));
if(Im(z)>=0) return asin2(z);
if(Im(z)<0) return asin3(z);
printf("asin0123(%24.16e +i* %24.16e) problem. do any!\n",Re(z),Im(z));
getchar();return 0; }
//

Comparison with asin

// In order to avoid confusion of names at the comparison of the TORI implementation of ArcSin above with the built-in function asin, the name of the function here is arcsin, not asin. Comparison of these two functions for 201 real values of the argument can be performed with the routine below:


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)
using namespace std;
#include <complex>
typedef complex<double> z_type;
#define Re(x) x.real()
#define Im(x) x.imag()
#define I z_type(0.,1.)
#include"arcsin.cin"

int main(){ int j,k,m,n; DB x,y, p,q, t; z_type z,c,d;
for(m=-101;m<102;m++) { x=.01*m; z=x; y=asin(x); c=arcsin(z); d=arcsin(z+I*1.e-8);
printf("%5.2f %20.16f %20.16f %20.16f %9.1e %9.1e\n",x,y,Re(c),Re(d), Re(c)-y, Re(d)-Re(c)); }
}

Keywords

ArcSin, C++, Inverse function

References