File:Facsma.jpg

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)
#include <complex>
//#define std::z_type complex<double>
typedef std::complex<double> z_type;
#define Re(x) x.real()
#define Im(x) x.imag()
#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");}

#define fac(z) exp(lofac(z))
int 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);
//FILE *o;o=fopen("facplo.eps","w");ado(o,0,0,98,144);
FILE *o;o=fopen("facsmall.eps","w");ado(o,0,0,460,1100);
#define M(x,y) fprintf(o,"%8.5f %8.5f M\n",0.+x, 0.+y);
#define L(x,y) fprintf(o,"%8.5f %8.5f L\n",0.+x, 0.+y);
fprintf(o,"110 10 translate\n 100 100 scale\n");
for(m=-1;m<4;m++) { if(m==0) {M(m,0)L(m,7)} else {M(m,0)L(m,6)}}
for(n=0;n<7;n++) { M( -1,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);}
        DO(m,45){x=-0.914+.1*m; z=x; y=Re(fac(z)); if(m==0)M(x,y) else L(x,y);}
fprintf(o,"1 setlinejoin 1 setlinecap .04 W 1 0 0 RGB S\n");
fprintf(o,"showpage\n%c%cTrailer",'%','%'); fclose(o);
system("epstopdf facsmall.eps");
system("open facsmall.pdf"); //for macintosh
  getchar(); system("killall Preview"); //for macintosh
}

\documentclass[12pt]{article}
\paperwidth 470pt
\paperheight 1140pt
\textwidth 1800pt
\textheight 1800pt
\topmargin -108pt
\oddsidemargin -78pt
\parindent 0pt
\pagestyle{empty}
\usepackage {graphics}
\usepackage{rotating}
\newcommand \rot {\begin{rotate}}
\newcommand \ero {\end{rotate}}
\newcommand \ing {\includegraphics}
\newcommand \sx {\scalebox}
\begin{document}
\begin{picture}(686,1106)
%\put(0,0) {\ing{facit}}
\put(0,0) {\ing{facsmall}}
\put(84,686){\sx{4.1}{$y$}}
\put(84,596){\sx{4}{$6$}}
\put(84,496){\sx{4}{$5$}}
\put(84,396){\sx{4}{$4$}}
\put(84,296){\sx{4}{$3$}}
\put(84,196){\sx{4}{$2$}}
\put(84, 96){\sx{4}{$1$}}
\put(98,-30){\sx{4}{$0$}}
\put(198,-30){\sx{4}{$1$}}
\put(298,-30){\sx{4}{$2$}}
\put(398,-30){\sx{4}{$3$}}
\put(448,-30){\sx{4.1}{$x$}}
\put(400,640){\sx{4.2}{\rot{84}$y\!=\!\mathrm{Factorial}(x)$\ero}}
\end{picture}
\end{document}