Factorial/Code/FactoReal

// FactoReal.cc, source of the figure // http://en.citizendium.org/wiki/Image:FactoReal.jpg // which represents the factorial of the real argument. // Copyleft 2009 by Dmitrii Kouznetsov. // Feel free to use, post, redistribute and even modify the code below, // but please indicate the source and the modifications if any. // Such indications help a lot to correct bugs if any.


 * 1) include 
 * 2) include 
 * 3) include 
 * 4) define DB double
 * 5) define DO(x,y) for(x=0;x
 * 7) define z_type complex
 * 8) define Re(x) x.real
 * 9) define Im(x) x.imag
 * 10) define I z_type(0.,1.)

DB expaunoc[31]= {1., 0.5772156649015329,   -0.6558780715202537,	 -0.04200263503409518,    0.16653861138229112, -0.04219773455554465,   -0.009621971527877027,   0.0072189432466631676, -0.0011651675918590183, -0.0002152416741149077,  0.00012805028238804805,-0.00002013485478102872,-1.2504934818746705e-6, 1.1330272320364543e-6, -2.0563384228733383e-7, 6.1160952968819515e-9,  5.00200766282674e-9, -1.1812748557105124e-9, 1.0434320074637071e-10, 7.782441358017422e-12, -3.696820627396846e-12, 5.10702591327572e-13, -2.0650148258027912e-14,-6.217248937900877e-15,  7.771561172376096e-16, -9.992007221626409e-16, -3.3306690738754696e-16, 5.551115123125783e-16, -1.1102230246251565e-16, 1.3322676295501878e-15, 9.992007221626409e-16 }; z_type expauno(z_type z) {int n,m; DB x,y; z_type s; s=expaunoc[24]; x=Re(z);if(x<-.9) return expauno(z+1.)-log(z+1.); if(x>.5) return expauno(z-1.)+log(z); y=Im(z); if(fabs(y)>.7)return expauno(z/2.)+expauno(z/2.-.5)+z*log(2.)-log(sqrt(M_PI)); for(n=23; n>=0; n--) { s*=z;s+=expaunoc[n]; } 			return -log(s); }

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 }; /* a[0]=1./12.; a[1]=1./30.; a[2]=53./210.; a[3]=195./371.; a[4]=22999./22737.; a[5]=29944523./19773142.; a[6]=109535241009./48264275462.; a[7]=29404527905795295658./9769214287853155785.; a[8]=455377030420113432210116914702./113084128923675014537885725485.; a[9]=26370812569397719001931992945645578779849./5271244267917980801966553649147604697542.; a[10]=152537496709054809881638897472985990866753853122697839./24274291553105128438297398108902195365373879212227726.; a[11]= too long... */ 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); }

z_type lofac(z_type z){DB x,y,r; x=Re(z); y=Im(z); if(fabs(y)>5 ) return fracti(z); if(x>0 && (x-3)*(x-3.)+y*y >25) return fracti(z); return expauno(z); }

void ado(FILE *O, int x, int y, int X, int Y) {	fprintf(O,"%c!PS-Adobe-2.0 EPSF-2.0\n",'%'); fprintf(O,"%c%cBoundingBox: %d %d %d %d\n",'%','%',x,y,X,Y); fprintf(O,"/M {moveto} bind def\n"); fprintf(O,"/L {lineto} bind def\n"); fprintf(O,"/S {stroke} bind def\n"); fprintf(O,"/s {show newpath} bind def\n"); fprintf(O,"/C {closepath} bind def\n"); fprintf(O,"/F {fill} bind def\n"); fprintf(O,"/times-Roman findfont 20 scalefont setfont\n"); fprintf(O,"/W {setlinewidth} bind def\n"); fprintf(O,"/RGB {setrgbcolor} bind def\n");}
 * 1) define fac(z) exp(lofac(z))

main{ int j,k,m,n; DB x,y, p,q, t; z_type z,c,d, cu,cd; FILE *o;o=fopen("FactoReal.eps","w");ado(o,0,0,98,144); fprintf(o,"41 71 translate\n 10 10 scale\n"); for(m=-3;m<4;m++) {	if(m==0){M(m,-6.2)L(m,6.2)} else	{M(m,-6)L(m,6)}			} for(n=-6;n<7;n++) {	M( -3,n)L(3,n)} fprintf(o,".006 W 0 0 0 RGB S\n"); DO(m,24){x=-3.9744+.03933*m; z=x; y=Re(fac(z)); if(m==0)M(x,y) else L(x,y);} DO(m,21){x=-2.921+.039*m; z=x; y=Re(fac(z)); if(m==0)M(x,y) else L(x,y);} DO(m,20){x=-1.832+.0347*m; z=x; y=Re(fac(z)); if(m==0)M(x,y) else L(x,y);} DO(m,41){x=-0.866+.0999*m; z=x; y=Re(fac(z)); if(m==0)M(x,y) else L(x,y);} fprintf(o,"1 setlinejoin 1 setlinecap .02 W 1 0 0 RGB S\n"); DO(m,176){x=-4+.05*m; z=x; y=Re(exp(-lofac(z))); if(m==0)M(x,y) else L(x,y);} fprintf(o,"1 setlinejoin 1 setlinecap .02 W 0 0 1 RGB S\n"); DO(m,61){y=.4616321450 +.04*m; x=Re(fac(y)); if(m==0)M(x,y) else L(x,y);} fprintf(o,"1 setlinejoin 1 setlinecap .02 W 0 .8 0 RGB S\n"); M(-1.504083  ,0) L(-1.504083,    Re(fac(-1.504083  )))  L(0,Re(fac(-1.504083))) M(-1.504083,   Re(1./fac(-1.504083  )))  L(0,1./Re(fac(-1.504083))) M( .4616321450,0) L( .4616321450,1./Re(fac(.4616321450)))  L(0,1/Re(fac(.4616321450))) M( .4616321450,Re(fac(.4616321450)))  L(0,Re(fac(.4616321450))) M(0, .4616321450) L( Re(fac(.4616321450)),.4616321450) L(  Re(fac(.4616321450)),0) fprintf(o,".01 W 1 0 1 RGB S\n"); fprintf(o,"showpage\n\%\%\%Trailer"); fclose(o); //system( "ggv FactoReal.eps"); //for linux system("open FactoReal.eps"); //for macintosh system("ps2pdf FactoReal.eps"); getchar; system("killall Preview"); //for macintosh } // end of heneratot of Image FactoReal.jpg
 * 1) define M(x,y) fprintf(o,"%8.5f %8.5f M\n",0.+x, 0.+y);
 * 2) define L(x,y) fprintf(o,"%8.5f %8.5f L\n",0.+x, 0.+y);