Приблизительное вычисление факториалов

{$N+}
 
const
  precision=16;
 
var
  n,result,i:extended;
  c:array[0..precision]of extended;
  d:array[0..precision]of extended;
  j:longint;
 
function power(a:extended;f:longint):extended;
var
  g:longint;
  h:extended;
begin
  if f = 0 then power:=1 else
  if f = 1 then power:=a else
  begin
    h:=a;
    for g:=2 to f do h:=h*a;
    power:=h;
  end;
end;
 
begin
  writeln('Precision: 6.2 <= n <= 1754.5, Ctrl-C - exit');
  c[0]:=1;
  d[0]:=1;
  c[1]:=1;
  d[1]:=12;
  c[2]:=1;
  d[2]:=288;
  c[3]:= -1.39e+2;
  d[3]:=  5.1840e+4;
  c[4]:= -5.71e+2;
  d[4]:=  2.488320e+6;
  c[5]:=  1.63879e+5;
  d[5]:=  2.09018880e+8;
  c[6]:=  5.246819e+6;
  d[6]:=  7.5246796800e+10;
  c[7]:= -5.34703531e+8;
  d[7]:=  9.02961561600e+11;
  c[8]:= -4.483131259e+9;
  d[8]:=  8.6684309913600e+13;
  c[9]:=  4.32261921612371e+14;
  d[9]:=  5.14904800886784000e+17;
  c[10]:= 6.232523202521089e+15;
  d[10]:= 8.6504006548979712000e+19;
  c[11]:=-2.5834629665134204969e+19;
  d[11]:= 1.3494625021640835072000e+22;
  c[12]:=-1.579029138854919086429e+21;
  d[12]:= 9.716130015581401251840000e+24;
  c[13]:= 7.46590869962651602203151e+23;
  d[13]:= 1.16593560186976815022080000e+26;
  c[14]:= 1.511513601028097903631961e+24;
  d[14]:= 2.798245444487443560529920000e+27;
  c[15]:=-8.849272268392873147705987190261e+30;
  d[15]:= 2.99692087104605205332754432000000e+32;
  c[16]:=-1.42801712490607530608130701097701e+32;
  d[16]:= 5.7540880724084199423888850944000000e+34;
  repeat
    readln(n);
    if n > 1754.5 then writeln('Please, no more 1754.5!') else
    if n = 0 then writeln('= 1') else
    if n < 0 then writeln('Please, positive only!') else
    begin
      i:=0;
      for j := precision downto 0 do i:=i+c[j]/d[j]/power(n,j);
      result:=exp(ln(n/exp(1))*n)*sqrt(2*n*pi)*i;
      writeln('=',result);
    end;
  until false;
end.

Эта тема была разделена из темы "Факториал числа - N!"

New version!

{$N+}
 
var
  n:extended;
  d:array[0..101]of extended;
 
function power(a:extended;f:longint):extended;
var
  g:longint;
  h:extended;
begin
  if f = 0 then power:=1 else
  begin
    h:=a;
    for g:=2 to f do h:=h*a;
    power:=h;
  end;
end;
 
function lnfact(z:extended):extended;
var
  r,i,old:extended;
  j:byte;
begin
  i:=0;
  old:=0;
  for j := 0 to 50 do if (ln(z)*j)>11352 then break else
  begin
    i:=i+d[j]/power(z,j);
    if old = i then break else old:=i;
  end;
  r:=ln(z/exp(1))*z+ln(sqrt(2*z*pi))+ln(i);
  lnfact:=r;
end;
 
function smallfact(x:extended):extended;
var
  r,z:extended;
  a,b:byte;
begin
  a:=trunc(x);
  z:=frac(x);
  r:=exp(lnfact(7+z));
  for b:=7 downto a+1 do r:=r/(b+z);
  smallfact:=r;
end;
 
function negfact(x:extended):extended;
var
  z:extended;
begin
  if (-x) < 7 then z:=smallfact(-x) else z:=exp(lnfact(-x));
  negfact:=pi/sin(n*pi)*n/z;
end;
 
function double(f:longint):extended;
var
  k:longint;
  r:extended;
begin
  r:=1;
  for k:=1 to (f-1) div 2 do r:=r*(k*2+1);
  double:=r;
end;
 
procedure filldata;
var
  t,k:byte;
  s:extended;
begin
  d[0]:=1;
  d[1]:=1;
  for t:=2 to 101 do
  begin
    s:=0;
    for k:=2 to t-1 do s:=s+(k*d[k]*d[t+1-k]);
    d[t]:=(d[t-1]-s)/(t+1);
  end;
  for t:=1 to 50 do d[t]:=double(2*t+1)*d[t*2+1];
end;
 
procedure bigfact;
var
  r1,r2:extended;
begin
   r1:=lnfact(n)/ln(10);
   r2:=exp(frac(r1)*ln(10));
   r1:=int(r1);
   writeln('n! = ',r2:1:17,'E+',r1:0:0);
end;
 
function fact(x:extended):extended;
begin
  if n < 0
  then fact:=negfact(x)
  else if n < 7
  then fact:=smallfact(x)
  else fact:=exp(lnfact(x));
end;
 
begin
  filldata;
  writeln('n!, Ctrl-C - exit');
  repeat
    write('n = ');
    readln(n);
    if (n < -1754.5) or (n > 1E+4928)
    then writeln('Sorry for n exceeds the limits of -1754.5 to 10^4928.')
    else if (n < 0) and (frac(n) = 0)
    then writeln('Sorry for n is negative integer; division by zero.')
    else if n <= 1754.5
    then writeln('n! =',fact(n):26)
    else bigfact;
  until false;
end.

Сообщение отредактировано: Scatman - 06.11.08, 06:26