[Back to MATH SWAG index]  [Back to Main SWAG index]  [Original]

(*
                Daniel Doubrovkine - dblock@infomaniak.ch
   part of the Expression Calculator 2.0 Multithread, destribute freely
              http://www.infomaniak.ch/~dblock/express.htm
  (ref. Wanner / Hairer - Analysis by it's History - University of Geneva)

  calculating a factor of any defined value inclusing negative and positive
      non integers! using Euler's (another swiss guy) Gamma function:
       n!=Gamma(n+1) and (n-1)!=Gamma(n)/(n-1) for negative values
        Gamma(Alpha):=Integral zero->infinity of x^(Alpha-1)*E^-x

 in GammaIntegral function TOTAL PRECISION (not just approx!) is reached by
   calculating an integral from 0 to 100 with a step if 0.01 from Gamma(4)
*)
function Gamma(alpha,step:extended):extended;
 function GammaIntegral(alpha:extended):extended;
  (*x^y:=e^(y*ln(x));*)
  function Power(base, exponent: extended): extended;
   begin
      Power:=exp(exponent*ln(base));
      end;
  function IntegralStep(x: extended):extended;
  begin
     (*Gamma Integral Step...just to have less mess*)
     IntegralStep:=power(x,alpha-1)/power(Exp(1),x);
     end;
 (*Gamma Integral*)
 var
   GammaTempIntegral:extended;
   l: extended;
 begin
   l:=0;
   GammaTempIntegral:=0;
   while l<100 do begin
         l:=l+Step;
         GammaTempIntegral:=GammaTempIntegral+IntegralStep(l)*step;
       end;
       GammaIntegral:=GammaTempIntegral;
       end;
(*Gamma*)
var
  NewGamma: extended;
  i: integer;
  begin
  if (alpha<=0) then begin
    if (trunc(alpha)=alpha) then begin
       writeln('factor results to infinite value');
       halt;
       end
       else begin
            NewGamma:=GammaIntegral(abs(1+frac(alpha))+1);
            for i:=0 to trunc(abs(alpha))-1 do begin
                NewGamma:=NewGamma/(frac(alpha)-i);
                end;
            Gamma:=NewGamma;
            exit;
            end
  end  else begin
     (*the following values return trailing decimals, this corrects it*)
     if alpha=1 then Gamma:=1 else
     if alpha=2 then Gamma:=1 else
     if alpha=3 then Gamma:=2 else
     Gamma:=GammaIntegral(alpha);
     end
end;

function Factor(n: extended):extended;
 begin
   Factor:=Gamma(n+1,0.01);
   end;

begin
     writeln(Factor(6):0);
end.

[Back to MATH SWAG index]  [Back to Main SWAG index]  [Original]