Gamma function: Difference between revisions

</pre>

=={{header|Prolog}}==

This version matches Wolfram Alpha to within a few digits at the end, so the last few digits are a bit off. There's an early check to stop evaluating coefficients once the desired accuracy is reached. Removing this check does not improve accuracy vs. Wolfram Alpha.

<lang Prolog>

gamma_coefficients(

[ 1.0000000000000000000000000000000000000000, 0.5772156649015328606065120900824024310422, -0.6558780715202538810770195151453904812798,

-0.0420026350340952355290039348754298187114, 0.1665386113822914895017007951021052357178, -0.0421977345555443367482083012891873913017,

-0.0096219715278769735621149216723481989754, 0.0072189432466630995423950103404465727099, -0.0011651675918590651121139710840183886668,

-0.0002152416741149509728157299630536478065, 0.0001280502823881161861531986263281643234, -0.0000201348547807882386556893914210218184,

-0.0000012504934821426706573453594738330922, 0.0000011330272319816958823741296203307449, -0.0000002056338416977607103450154130020573,

0.0000000061160951044814158178624986828553, 0.0000000050020076444692229300556650480600, -0.0000000011812745704870201445881265654365,

0.0000000001043426711691100510491540332312, 0.0000000000077822634399050712540499373114, -0.0000000000036968056186422057081878158781,

0.0000000000005100370287454475979015481323, -0.0000000000000205832605356650678322242954, -0.0000000000000053481225394230179823700173,

0.0000000000000012267786282382607901588938, -0.0000000000000001181259301697458769513765, 0.0000000000000000011866922547516003325798,

0.0000000000000000014123806553180317815558, -0.0000000000000000002298745684435370206592, 0.0000000000000000000171440632192733743338

]).

tolerance(1e-17).

gamma(X, _) :- X =< 0.0, !, fail.

gamma(X, Y) :-

X < 1.0, small_gamma(X, Y), !.

gamma(1, 1) :- !.

gamma(1.0, 1) :- !.

gamma(X, Y) :-

X1 is X - 1,

gamma(X1, Y1),

Y is X1 * Y1.

small_gamma(X, Y) :-

gamma_coefficients(Cs),

recip_gamma(X, 1.0, Cs, 1.0, 0.0, Y0),

Y is 1 / Y0.

recip_gamma(_, _, [], _, Y, Y) :- !.

recip_gamma(_, _, [], X0, X1, Y) :- tolerance(Tol), abs(X1 - X0) < Tol, Y = X1, !. % early exit

recip_gamma(X, PrevPow, [C|Cs], _, X1, Y) :-

Power is PrevPow * X,

X2 is X1 + C*Power,

recip_gamma(X, Power, Cs, X1, X2, Y).

</lang>

{{Out}}

<pre>

% see how close gamma(0.5) is to the square root of pi.

?- gamma(0.5,X), Y is sqrt(pi), Err is abs(X - Y).

X = 1.772453850905516,

Y = 1.7724538509055159,

Err = 2.220446049250313e-16.

?- gamma(1.5,X).

X = 0.886226925452758.

?- gamma(4.9,X).

X = 20.667385961857857.

?- gamma(5,X).

X = 24.

?- gamma(5.01,X).

X = 24.364473447872836.

?- gamma(6.9,X).

X = 597.4941281573107.

?- gamma(6.95,X).

X = 655.7662628554252.

?- gamma(7,X).

X = 720.

% 100!

?- gamma(101.0,X).

X = 9.33262154439441e+157.

% Note when passed integer, gamma(101) returns full big int precision

?- gamma(101,X).

X = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000.

?- gamma(100.98,X).

X = 8.510619261391532e+157.