Acosc1.cin

From TORI
Jump to navigation Jump to search

// Acosc1.cin is set of routines for evaluation acosc1, whih is first branch of function ArcCosc behind the cut line at the negative part of the real axis.

z_type acos(z_type z){
if(Im(z)<0){if(Re(z)>=0){return  I*log( z + sqrt(z*z-1.) );} 
                   else{return  I*log( z - sqrt(z*z-1.) );}}
           if(Re(z)>=0){return -I*log( z + sqrt(z*z-1.) );} 
                  else {return -I*log( z - sqrt(z*z-1.) );} }
z_type cosc(z_type z) {return cos(z)/z;}
z_type cosp(z_type z) {return (-sin(z) - cos(z)/z)/z ;}
z_type cohc(z_type z) {return cosh(z)/z ;}
z_type cohp(z_type z) {return (sinh(z)-cosh(z)/z)/z ;}
z_type acoscL(z_type z){ int n; z_type s,q; z*=-I; q=I*sqrt(1.50887956153832-z);
  s=q*1.1512978931181814 + 1.199678640257734; DO(n,6) s+= (z-cohc(s))/cohp(s);
  return -I*s; }
z_type acoscR(z_type z) {int n; z_type s= (1.-0.5/(z*z))/z;
       DO(n,5) s+=(z-cosc(s))/cosp(s); return s;}
z_type acoscB(z_type z){ z_type t=0.33650841691839534+z, u=sqrt(t), s; int n;
s=    2.798386045783887
+u*(-2.437906425896532
+u*( 0.7079542331649882 
+u*(-0.5009330133042798
+u*( 0.5714459932734446    ))));
DO(n,6) s+=(z-cosc(s))/cosp(s); return s; }
DB Sazae0= 2.798386045783887;
DB Tarao0=-0.33650841691839534;
DB Sazae1= 6.1212504668980685;
DB Tarao1= 0.161228034325064;
z_type acosc(z_type z){DB x1=-0.33650841691839534,x=Re(z),y=Im(z),yy=y*y,r; 
r=x-x1;r*=r;r+=yy; if(r < 1.8 )     return acoscB(z); 
r=x+2.;r*=r;r+=yy; if(r>8. && x>=0) return acoscR(z);
if(y >= 0) return acoscL(z);
      return conj(acoscL(conj(z))); }
z_type acosc1R(z_type z){int n;  z_type s,t,u; t=Tarao1-z; u=sqrt(t);
s=Sazae1- 3.522043522040332*u;
 s+=(z-cosc(s))/cosp(s);
DO(n,5) s+=(z-cosc(s))/cosp(s); return s; }
z_type acosc1B(z_type z){ z_type t=0.33650841691839534+z, u=sqrt(t), s; int n;
s=    2.798386045783887
+u*( 2.437906425896532
+u*( 0.7079542331649882 
+.16*u*( 0.5009330133042798 // I have no idea about this..
+u*( 0.5714459932734446    ))));
DO(n,6) s+=(z-cosc(s))/cosp(s);  return s; }
z_type acosc1H(z_type z){int n; z_type s=2*M_PI - acos(2*M_PI*z);
 DO(n,6) s+=(z-cosc(s))/cosp(s);  return s; }
z_type acosc1(z_type z){ DB x=Re(z), y=Im(z);
if( abs(z) > 1 ) return acosc1H(z); 
if( x < 0 )       return acosc1B(z);
                  return acosc1R(z);  }
// The acosc1(z) returns of order of 14 correct decimal digits while \(|z|<6 \). 
// Let me know if any problem with the use. 
// Copileft 2012 by Dmitrii Kouznetsov