Conference rusure::math

I suspect it uses a straight-forward algorithm implemented in H-floating.
It doesn't take long; try the following Pascal program:

program x939 ( input, output ) ;
type hfloat = quadruple;

function factorial (n: integer): hfloat;
var i : integer; x,y : hfloat;
begin
y := 1;
for i := 1 to n do begin x := i; y := y * x; end;
factorial := y;
end;

procedure main;
var i : integer;
begin
for i := 1700 to 1800 do writeln (i, ' ', factorial(i));
end;

begin
main;
end.

    Or, possibly, Stirling's approximation:
    
    N! ~= (N/e)^N * SQRT(2*PI*N)

There is a function (gamma?) defined by an integral which evaluates to
(n-1)! or (n+1)! (can't remember which) for integer n.  The integral
is probably harder to do than any straightforward iteration, but it may
be that Stirling's formula (or other approximations to n!) come from applying
simplifying assumptions to evaluate this integral for large n.

If I can find the form of the integrand, I will post it.

The form of the gamma function integral is:

gamma( u ) = integral(  0,  oo, e^(-t) * t^(u-1), t )

with 'integral' having arguments as in:

integral( lower_limit, upper_limit, integrand, variable_of_integration ) 

and 'oo' denoting infinity.

Reference:

p. 362, problem #10
Elementary Partial Differential Equations
Berg, Paul W. & McGregor, James L.
Holden-Day
1966

					     6
Yikes. Those were fast responses. Thanks a 10 !.

Re: .1 Thanks, I'll try that if I can figure out how to program. :-/
Re: .2 Stirlings' formula is what I was referring to in .0 as the not very
	good approximation, I just couldn't remember the guy's name.
Re: .4 Thanks, I'll try that if I can remember how to evaluate an integral. :-\

Mike

    re .2 Actually Stirling's formula is terrible.  It in fact diverges
    from the factorial! (pardon the punctuational pun) The ratio of
    Stirling's formula to the factorial converges to one as n gets large
    but so does (n^2+n)/n^2, the top is only close to the bottom when
    n is very small (like zero).
                                
    There are more accurate forms of Stirling's formula, I believe the
    Monthly had an article a year or so ago about these, but they all
    diverge from the factorial just the same as the simple one does.
    
    re .4 The gamma function is not much help for evaluating factorials
    for various reasons but it allows an extension of the factorial
    to all real numbers except (if I remember correctly) negative integers
    where it becomes undefined.
    
    Chuck

    Here is a routine which I wrote in VAX Pascal for calculating
    factorials.  Since my use freqeuntly involved large factorials
    which could cause floating overflow and which were to be divided
    by other large quantities, I chose to calculate the (natural) log
    of the factorial, and then call the exponential function.  It
    might be possible to write a "clever" exponential function which
    takes advantage of the fact that anything non-integral is noise,
    but it is probably not worth it.  The whole thing could be converted
    quite easily to calculate factorial rather than log factorial
    pretty easily -- just some changed constants and multiplications
    and simple algebraic manipulations (the tricky part is calculating
    Z^Z efficiently, see Knuth).
    
    I didn't care much about the absolute error, only the relative
    error, and have not checked the convergence properties for the
    absolute error for this extended version of Stirlings approximation.
    
    Numerical Recipies by Press, Flannery, Teukolsky and Vetterling
    contains another method of calculating large factorials, by using
    a general, very accurate extension of Stirlings approximation
    to the gamma function (actually the log gamma function).  Like
    the code below, they keep a table of relatively small values of
    factorial which is built as needed.
    
    Under other circumstances I would probably hard-wire the table
    rather than calculating it as needed.  When a large number of
    factorials are needed, the incremental cost of building the table
    is trivial.
    
    I should warn people that I made no effort to examine the effects
    of round-off error in this calculation.  I simply coded it naively
    and tried it out against some values calculated via brute force
    and since it seemed to work I kept it.
    
    					Topher
-----------------------------------------------------------------------------
MODULE log_factorial;
 

CONST
    C0 = 9.189385332046727417803296D-1;  {ln(2*pi)/2}
    C1 = 1D0/12D0;
    
    calc_thresh = 100;
 
VAR
    table : ARRAY [1..calc_thresh] OF Double	    {Holds values already}
			:= (calc_thresh of 0D0);    {  calculated}
    table_limit: Integer := 1;			    {Limit of values already
						       calculated}
 
 
[GLOBAL] FUNCTION Ln_Factorial(N: Integer): Double;

    VAR
	Z: Double;
	i: Integer;
	double_i: Double;
    BEGIN
    IF N > calc_thresh
    THEN
	{ The following approximation, based on Stirlings formula, has
	  a maximum absolute error, with the above threshold > 99 of
	  less than 3*10^-9.  The formula was taken from Abramowitz and
	  Stegun "Handbook of Mathematical Functions" which is available
	  both from the US Government Printing Office and from Dover Books.
	  If this accuracy is not acceptable a substitution can be made
	  resulting in a (theoretical) maximum absolute error of less than
	  8*10^-14.  This substitution consists of replacing the term
	  "C1/Z" with "(C1*Z3 - C2)/Z3" where Z3 = Z**3, and C2 = "1D0/360D0".}
	BEGIN
	Z := N + 1.0D0;
	Ln_Factorial := C0 -
			Z +
			(Z - 0.5D0)*LN(Z) +
			C1/Z
	END
    ELSE IF N < 2
    THEN
	Ln_Factorial := 0D0
    ELSE
	BEGIN
	IF N > table_limit
	THEN
	    BEGIN
	    FOR i := table_limit + 1 TO N DO
		BEGIN
		double_i := i;
		Table[i] := Table[i-1] + LN(double_i);
		END;
	    table_limit := N;
	    END;

	Ln_Factorial := Table[N];
	END
    END;
END.

     Stirling's formula "diverges" from the factorial in the
     sense that
     
          | Stirling(n) - n! | -> oo as n -> oo
     
     However, the percentage error converges to zero:
     
          | Stirling(n) - n! |
          -------------------- -> 0 as n -> oo
          |        n!        |
     
     So more and more decimal places of the approximation
     are correct for larger and larger values of n.
     
     Dan

    Having worked for DEC for 8 years, my english is nearly gone but
    I thought that's what I said.

     Sorry, you're right, that is what you said.  I remembered
     disagreeing with something you said when I first read
     it, but by the time I replied I lost track of what.
     
     The part I disagreed with is calling it a "terrible"
     approximation, since most uses I've had for the
     approximation use the fact that LIM ratio -> 1 and
     it doesn't matter that LIM | difference | -> oo.
     
     So it was your emphasis on a fact that I considered
     to be irrelevant that I meant to dispute.
     
     Dan

    Re. .10
    You are very right.  Normally Stirling's formula is used for computing
    probabilities where ratios matter in which case it is very very
    good.  The one I remember best is computing the pressure of a gas
    from the total energy (kinetic theory).  Here the limit desired
    is "exact."
    
    Chuck
    

>...the largest factorial the calculator can do on our machine is 1754!.
>This results in the answer .197926189010501 *10^4931

If you *really* dig factorials, MAPLE will give you all 4930 digits of
1754! and will give larger factorials to complete precision. Of course,
they take a little time, but hey... 

Lynn 

>I have been playing with Factorials lately and found it neat that the 
>calculator in All-In-1 can spit out a very large factorial to 15 significant
>digits in an instant. I've looked around in the literature for algorithms
>which approximate factorials, but have found only one, and that one was very
>inaccurate. does anyone know what algorithm All-In-1 uses? How accurate are
>the factorials this program? do you know any other algorithms like this? 

It increments from 2 to ARG.
and does a floatH add 1 to I, and a floatH multiply by I.

BTW, a nit: ALL-IN-1 is spelled in all caps. Please help protect 
a valuable trademark. Thanks.

John Frazer - Charlotte OP/ACES.

RE:       < Note 939.14 by CLT::GILBERT "Multiple inheritence happens" >
                    -< re .-1:  I corrected the note title >-

Thanks for the correction. Sorry for the mistake. 

Mike

						   TM
The current Digital advertisement mentions All-In-1  , spelled just like that.

>< Note 939.16 by CLT::GILBERT "Multiple inheritence happens" >
>
>
>						   TM
>The current Digital advertisement mentions All-In-1  , spelled just like that.
>

Yeah, they never get it right either! %^)

Cheers,

John F.

>>						   TM
>>The current Digital advertisement mentions All-In-1  , spelled just like that.
>>
>
>Yeah, they never get it right either! %^)

	Given the amount of effort that people continually expend in reminding
one another how to spell ALL-IN-1, several thoughts occur to me:

	(a) Why not copyright all the reasonable mutations of ALL-IN-1 at
the outset (eg Allinone, All_in_1 etc)?

	(b) Why not pick a name which is not open to such mal-rendering (eg
ENTIRE, perhaps?)

	(c) If copyrights are case-sensitive, can I frivolously gain a lot of
money by forcing large corporations to buy off me copyrights such as:
ibM, UNix, Texa�o, McDon�lds, Su�?

Andrew

>	(a) Why not copyright all the reasonable mutations of ALL-IN-1 at
>the outset (eg Allinone, All_in_1 etc)?

Now if we can list all of those, we may have found a
justification for having a Factorial function in ALL-IN-1 in the
first place. %^) 

Cheers,

John F.