Fac.cin

// Fac.cin is the C++ code of function "fac" that evaluates Factorial of the complex argument.

// z_type is supposed to be defined as complex double (although, perhaps, it might have another meaning, for example, just complex or just "double" for the evaluation for the real values of the argument).

z_type fracti(z_type z){ z_type s; int n;  DB a[17]=
{0.0833333333333333333, 0.0333333333333333333,  .252380952380952381, .525606469002695418, 
1.01152306812684171,   1.51747364915328740,   2.26948897420495996, 3.00991738325939817, 
4.02688719234390123,   5.00276808075403005,   6.28391137081578218, 7.49591912238403393, 
9.04066023436772670,  10.4893036545094823,   12.2971936103862059, 13.9828769539924302, 
16.0535514167049355 };
s=a[16]/(z+19./(z+25./(z)));  for(n=15;n>=0;n--) s=a[n]/(z+s);
return  s + log(2.*M_PI)/2. - z + (z+.5)*log(z);
} //  logfactorial for large values of argument except vicinity of negative part of real axis)

z_type infac0(z_type z){ z_type s; int n;  DB c[28]={ 1.,
 0.57721566490153286061,        -0.65587807152025388108,
-0.042002635034095235529,        0.16653861138229148950,
-0.042197734555544336748,       -0.0096219715278769735621,
0.0072189432466630995424,      -0.0011651675918590651121,
-0.00021524167411495097282,      0.00012805028238811618615,
-0.000020134854780788238656,    -0.0000012504934821426706573,
0.0000011330272319816958824,   -2.0563384169776071035e-7,
6.1160951044814158179e-9,       5.0020076444692229301e-9,
-1.1812745704870201446e-9,       1.0434267116911005105e-10,
7.7822634399050712540e-12,     -3.6968056186422057082e-12,
5.1003702874544759790e-13,     -2.0583260535665067832e-14,
-5.3481225394230179824e-15,      1.2267786282382607902e-15,
-1.1812593016974587695e-16,      1.1866922547516003326e-18,
1.4123806553180317816e-18};
s=c[27]*z; for(n=26;n>0;n--) {s+=c[n]; s*=z;}
s+=c[0];  return s;}

z_type fac0(z_type z){ return 1./infac0(z);}

z_type expaun(z_type z) {int n,m; DB x,y;
               x=Re(z);if(x<-.5) return expaun(z+1.)-log(z+1.);
                       if(x>.6) return expaun(z-1.)+log(z);
               y=Im(z); if(fabs(y)>1.4)return expaun(z/2.)+expaun(z/2.-.5)+z*log(2.)-log(sqrt(M_PI));
                       return -log(infac0(z)); }

z_type lof(z_type z){DB x,y; x=Re(z); y=Im(z);
       if(fabs(y)>5. ) return fracti(z);
       if(x>0 && x*x+y*y>25.) return fracti(z);
       return expaun(z);  } // lof(z) returns 16 digits of complex logfactorial.

 z_type infac1(z_type z){return infac0(z/2.)*infac0((z-1.)/2.)*sqrt(M_PI)/exp(log(2.)*z);}
 z_type infac2(z_type z){return infac1(z/2.)*infac1((z-1.)/2.)*sqrt(M_PI)/exp(log(2.)*z);}
 z_type infac3(z_type z){return infac2(z/2.)*infac2((z-1.)/2.)*sqrt(M_PI)/exp(log(2.)*z);}
 z_type inhalf(z_type z){DB x=Re(z); DB y=Im(z); DB r=x*x+y*y;
                       if(r<2.) return infac0(z); 
                       if(r<5.) return infac1(z);
                                return infac2(z); }

z_type infacmi(z_type z){ if(Re(z)> 1.) return infacmi(z-1.)/z;     return inhalf(z);}
z_type infaclu(z_type z){ if(Re(z)<-.5) return infaclu(z+1.)*(z+1.);return inhalf(z);}

z_type infac(z_type z){DB x=Re(z),y=Im(z),t=x*x+y*y; if(t<1.)return infac0(z);
       if( fabs(y)> 5. || (x>0 && t>25) )              return exp(-fracti(z));
       if( x>0 ) return infacmi(z);
                 return infaclu(z);}

z_type fac(z_type z){ DB x=Re(z),y=Im(z),t=x*x+y*y; if(t<2.)return 1./infac0(z);
               if( (x>0. && t>25.) || fabs(y)>5.)         return exp(fracti(z));
               if(x>0) return 1./infacmi(z);
                       return 1./infaclu(z);}

// Copyleft 2009 by Dmitrii Kouznetsov. The free ue is allowed but attribute the source at the reproduction.

// Such attribution helps to trace the erors if any. // // //