-
Notifications
You must be signed in to change notification settings - Fork 3
/
gamma.adb
111 lines (109 loc) · 5.32 KB
/
gamma.adb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
-- --
-- function Copyright (c) Dmitry A. Kazakov --
-- Gamma Luebeck --
-- Instantiation Spring, 2009 --
-- --
-- Last revision : 15:03 28 Mar 2009 --
-- --
-- This library is free software; you can redistribute it and/or --
-- modify it under the terms of the GNU General Public License as --
-- published by the Free Software Foundation; either version 2 of --
-- the License, or (at your option) any later version. This library --
-- is distributed in the hope that it will be useful, but WITHOUT --
-- ANY WARRANTY; without even the implied warranty of --
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU --
-- General Public License for more details. You should have --
-- received a copy of the GNU General Public License along with --
-- this library; if not, write to the Free Software Foundation, --
-- Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. --
-- --
-- As a special exception, if other files instantiate generics from --
-- this unit, or you link this unit with other files to produce an --
-- executable, this unit does not by itself cause the resulting --
-- executable to be covered by the GNU General Public License. This --
-- exception does not however invalidate any other reasons why the --
-- executable file might be covered by the GNU Public License. --
--____________________________________________________________________--
with Ada.Numerics.Long_Elementary_Functions;
with Long_Float_Chebyshev_Polynomials;
--
-- The implementation of Gamma function uses
--
-- ]-oo, 0[ Euler's reflection formula and Gamma(-x)
-- [0, 1[ Shifted first kind Chebyshev series of Gamma(x)/x
-- [1, 2] Shifted first kind Chebyshev series of Gamma(x+1)
-- ]2, oo[ Shifted first kind Chebyshev series of ln Gamma (x)
--
-- Source: Mathematical functions and their approximations by Yudell L.
-- Luke. 1975
--
function Gamma (X : Float) return Float is
use Ada.Numerics;
use Long_Elementary_Functions;
use Long_Float_Chebyshev_Polynomials;
Ln_2Pi_2 : constant := -- ln 2Pi / 2
( ( 1.1447_29885_84940_01741_43427 -- ln Pi
+ 0.6931_47180_55994_53094_37232_1 -- ln 2
)
/ 2.0
);
C1 : constant Coefficients := -- gamma(x+1), T*n(x), 0 <= x <= 1
( 0.94178_55977_95494_66571, 0.00441_53813_24841_00676,
0.05685_04368_15993_63379, -0.00421_98353_96418_56050,
0.00132_68081_81212_46022, -0.00018_93024_52979_88804,
0.00003_60692_53274_41245, -0.00000_60567_61904_46086,
0.00000_10558_29546_30228, -0.00000_01811_96736_55424,
0.00000_00311_77249_64715, -0.00000_00053_54219_63902,
0.00000_00009_19327_55199, -0.00000_00001_57794_12803,
0.00000_00000_27079_80623, -0.00000_00000_04646_81865,
0.00000_00000_00797_33502, -0.00000_00000_00136_80782,
0.00000_00000_00023_47319, -0.00000_00000_00004_02743,
0.00000_00000_00000_69101, -0.00000_00000_00000_11856,
0.00000_00000_00000_02034 -0.00000_00000_00000_00349,
0.00000_00000_00000_00060, -0.00000_00000_00000_00010,
0.00000_00000_00000_00002
);
C2 : constant Coefficients := -- ln gamma(x), T2n, x > 1
( 0.98575_15540_05098, -0.01357_51199_40355,
0.00060_97577_84871, -0.00005_47947_47404,
0.00000_72256_71298, -0.00000_12227_85636,
0.00000_02472_06061, -0.00000_00571_76475,
0.00000_00147_06251, -0.00000_00041_24183,
0.00000_00012_43043, -0.00000_00003_98333,
0.00000_00001_34579, -0.00000_00000_47620,
0.00000_00000_17553, -0.00000_00000_06710,
0.00000_00000_02651, -0.00000_00000_01079,
0.00000_00000_00451, -0.00000_00000_00193,
0.00000_00000_00085, -0.00000_00000_00040,
0.00000_00000_00017, -0.00000_00000_00008,
0.00000_00000_00004, -0.00000_00000_00002,
0.00000_00000_00001
);
function Gamma (X : Long_Float) return Long_Float is
pragma Inline (Gamma);
begin
if X <= 1.0 then
return Sum_Shifted_T (X, C1) / X;
elsif X <= 2.0 then
return Sum_Shifted_T (X - 1.0, C1);
else
return
( Exp
( (X - 0.5) * Log (X)
- X
+ Ln_2Pi_2
+ Sum_Even_T (1.0 / X, C2) / (12.0 * X)
) );
end if;
end Gamma;
begin
if X < 0.0 then
return
Float -- Euler's reflection formula
( Pi
/ (sin (Pi * Long_Float (X)) * Gamma (1.0 - Long_Float (X)))
);
else
return Float (Gamma (Long_Float (X)));
end if;
end Gamma;