MCPcopy Create free account
hub / github.com/BoevaLab/FREEC / gamma

Function gamma

src/gammafunc.cpp:52–146  ·  view source on GitHub ↗

Gamma function Input parameters: X - argument Domain: 0 < X < 171.6 -170 < X < 0, X is not an integer. Relative error: arithmetic domain # trials peak rms IEEE -170,-33 20000 2.3e-15 3.3e-16 IEEE -33, 33 20000 9.4e-16 2.2e-16 IEEE 33, 171.6 20000 2.3e-15 3.2e-16 Cephes Math Library Release 2.8

Source from the content-addressed store, hash-verified

50Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
51*************************************************************************/
52double gamma(double x)
53{
54#ifndef ALGLIB_INTERCEPTS_SPECFUNCS
55 double result;
56 double p;
57 double pp;
58 double q;
59 double qq;
60 double z;
61 int i;
62 double sgngam;
63
64 sgngam = 1;
65 q = fabs(x);
66 if( ap::fp_greater(q,33.0) )
67 {
68 if( ap::fp_less(x,0.0) )
69 {
70 p = ap::ifloor(q);
71 i = ap::round_f(p);
72 if( i%2==0 )
73 {
74 sgngam = -1;
75 }
76 z = q-p;
77 if( ap::fp_greater(z,0.5) )
78 {
79 p = p+1;
80 z = q-p;
81 }
82 z = q*sin(ap::pi()*z);
83 z = fabs(z);
84 z = ap::pi()/(z*gammastirf(q));
85 }
86 else
87 {
88 z = gammastirf(x);
89 }
90 result = sgngam*z;
91 return result;
92 }
93 z = 1;
94 while(ap::fp_greater_eq(x,3))
95 {
96 x = x-1;
97 z = z*x;
98 }
99 while(ap::fp_less(x,0))
100 {
101 if( ap::fp_greater(x,-0.000000001) )
102 {
103 result = z/((1+0.5772156649015329*x)*x);
104 return result;
105 }
106 z = z/x;
107 x = x+1;
108 }
109 while(ap::fp_less(x,2))

Callers 2

incompletebetaFunction · 0.85
incompletebetapsFunction · 0.85

Calls 2

round_fFunction · 0.85
gammastirfFunction · 0.85

Tested by

no test coverage detected