Facp.cin

// facp.cin is complex dopble routine for evaluation of derivative of Factorial function. The routines required are copypasted below

z_type infp0(z_type z)
{ int n; z_type s; DB c[30]=
{0.57721566490153286061,        -1.3117561430405077622,
-0.12600790510228570659,         0.66615444552916595801,
-0.21098867277772168374,        -0.057731829167261841373,
 0.050532602726641696797,       -0.0093213407348725208969,
-0.0019371750670345587553,       0.0012805028238811618615,
-0.00022148340258867062521,     -0.000015005921785712047888,
 0.000014729354015762046471,    -2.8788737837686499448e-6,
 9.1741426567221237268e-8,       8.0032122311507566881e-8,
-2.0081667698279342458e-8,       1.8781680810439809189e-9,
 1.4786300535819635383e-10,     -7.3936112372844114164e-11,
 1.0710777603654399556e-11,     -4.5283173178463149231e-13,
-1.2300681840672941359e-13,      2.9442687077718258964e-14,
-2.9531482542436469238e-15,      3.0853998623541608647e-17,
 3.8134277693586858102e-17,     -6.4364879164190365785e-18,
 4.9717783335892785568e-19,      4.0120551914810793446e-21};
s=c[29]; for(n=28;n>=0;n--){ s*=z; s+=c[n];}   return s;}

z_type lofp0(z_type z){ return -infp0(z)/infac(z); }
z_type facp0(z_type z){ z_type f=fac(z); f*=f;  return -infp0(z)*f; }

z_type lofpa(z_type z){ int n; DB q[11];  q[0] = 12.; q[1] = 5./6.;  q[2] = 252./79.;  
q[3] = 6241./14460.; q[4 ]= 7666692./4146631.; q[5 ]= 179081182865./612465549066.;
q[6 ]= 4881681043696812./3754087889491759.;
q[7 ]= 86960333299682003491937./392729697097736725384440.;
q[8 ]= 378191910699307315313565647105916./377413323237205130354503096392253.;
q[9 ]= 696148976661357653747206985359295786942014225./
      3903889440300118372577892204070110729027524454.;
q[10]= 36675782764501469367480729990524142326314524131790623634298644./
      45019657243089322180478800624616560743983830599801241354133773.;
z_type c=1./(z*z),   s=c/q[10];   for(n=9;n>=0;n--) s=c/(q[n]+s);
return -s + .5/z + log(z);}

       z_type lofp2(z_type z){ return log(2.)+(lofp0(z/2.-.5)+lofp0(z/2.))/2.;}
       z_type lofp4(z_type z){ return log(2.)+(lofp2(z/2.-.5)+lofp2(z/2.))/2.;}
       z_type lofp8(z_type z){ return log(2.)+(lofp4(z/2.-.5)+lofp4(z/2.))/2.;}

z_type lofp1(z_type z){DB x=Re(z),y=Im(z), t=x*x+y*y; 
               if(x>1) return lofp1(z-1.)+1./z;
               if(x<-.5) return lofp1(z+1.)-1./(z+1.);
               if(t<2.) return lofp0(z); 
               return lofp4(z); }

z_type lofp(z_type z){DB x=Re(z),y=Im(z), u=y*y;
if(x>=0. && (x+1.)*(x+1.)+u> 30.)return lofpa(z);
if(x<=0. && (x-1.)*(x-1.)+u> 30.)return lofpa(-z)+1./z-M_PI/tan(M_PI*z);
return lofp1(z);} 

z_type infp(z_type z){ DB x=Re(z),y=Im(z);
               if(x*x+y*y<2.)  return infp0(z);
                               return -infac(z)*lofp(z);}

z_type facp(z_type z){ DB x=Re(z),y=Im(z), u=x*x+y*y; z_type c;
if(u<2){c=infac0(z);    return -infp0(z)/(c*c);}
       return fac(z)*lofp(z);
//if(x>0|| fabs(y)>2.)return M_PI*insincp(M_PI*z)*infp(-z)-insinc(M_PI*z)*infp(-z);
}

//