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); }

//