Susin.cin

// susin.cin is the complex double implementation in C++ of function SuSin, that is superfunction of function sin, \(\mathrm{SuSin}(z\!+\!1)=\sin(\mathrm{SuSin}(z))\).

//The routine is tested to return of order of 14 correct signigificant figures, although the primitive implementation of arcsin used could be improved.

//

//DB mans double
//z_type means complex<double>
//DO(x,y) means for(x=0;x<y;x++)

z_type susin0(z_type z){ z_type d,c,L;
DB A[9][9]={
{1.0000000000000000, 0.0000000000000000,0.0000000000000000, 0.0000000000000000,0.0000000000000000, 0.0000000000000000,.0000000000000000, .0000000000000000,.000000000000000},
{0.0000000000000000,-0.3000000000000000,0.0000000000000000, 0.0000000000000000,0.0000000000000000, 0.0000000000000000,.0000000000000000, .0000000000000000,.000000000000000},
{0.1128571428571429,-0.1800000000000000,0.1350000000000000, 0.0000000000000000,0.0000000000000000, 0.0000000000000000,.0000000000000000, .0000000000000000,.000000000000000},
{0.1174285714285714,-0.2772857142857143,0.2160000000000000,-0.0675000000000000,0.0000000000000000, 0.0000000000000000,.0000000000000000, .0000000000000000,.000000000000000},
{0.1490034322820037,-0.4129714285714285,0.4207500000000000,-0.1917000000000000,0.0354375000000000, 0.0000000000000000,.0000000000000000, .0000000000000000,.000000000000000},
{0.2002934593977451,-0.6500921243042672,0.8099614285714286,-0.4936950000000000,0.1506600000000000,-0.0191362500000000,.0000000000000000, .0000000000000000,.000000000000000},
{0.2806962658586312,-1.0510236905951191,1.5586288622448981,-1.1871745714285715,0.4976943750000000,-0.1109173500000000,.0105249375000000, .0000000000000000,.000000000000000},
{0.4048298271172445,-1.7253296512057330,2.9846735140074210,-2.7385222637755100,1.4561118321428572,-0.4547520225000000,.0784112400000000,-.0058638937500000,.000000000000000},
{0.5969517517175638,-2.8569320127510403,5.6727958236173519,-6.1201236292764376,3.9545046460331634,-1.5833518624285714,.3881107608750000,-.0539255619642857,.003298440234375}};
c=1./z;
L=log(z);
d=(((((((A[8][8]*L+A[8][7])*L+A[8][6])*L+A[8][5])*L+A[8][4])*L+A[8][3])*L+A[8][2])*L+A[8][1])*L+A[8][0];d*=c;
d+=((((((A[7][7]*L+A[7][6])*L+A[7][5])*L+A[7][4])*L+A[7][3])*L+A[7][2])*L+A[7][1])*L+A[7][0]; d*=c;
d+= (((((A[6][6]*L+A[6][5])*L+A[6][4])*L+A[6][3])*L+A[6][2])*L+A[6][1])*L+A[6][0]; d*=c;
d+= ((((A[5][5]*L+A[5][4])*L+A[5][3])*L+A[5][2])*L+A[5][1])*L+A[5][0]; d*=c;
d+= (((A[4][4]*L+A[4][3])*L+A[4][2])*L+A[4][1])*L+A[4][0]; d*=c;
d+= ((A[3][3]*L+A[3][2])*L+A[3][1])*L+A[3][0]; d*=c;
d+= (A[2][2]*L+A[2][1])*L+A[2][0]; d*=c;
d+= A[1][1]*L; d*=c;
d+=1.; return sqrt(3./z)*d;
}

z_type susin(z_type z){ z_type d; int m,n;
//z+= 1.4304553465285;
z+= 1.4304553465288;
DB y=abs(Im(z)); if(y>40.) return susin0(z);
DB x=Re(z);
m=40-x; if(m<0)m=0;
d=susin0(z+(0.+m));
DO(n,m) d=arcsin(d);
return d;
}

//

// Category:C++