/*
* Open Source Physics software is free software as described near the bottom of this code file.
*
* For additional information and documentation on Open Source Physics please see:
* <http://www.opensourcephysics.org/>
*/
/*
* Gamma function by J E Hasbun for the OSP library.
* Ref : http : //www.alglib.net/specialfunctions
* Using the Algorithm of Stephen L. Moshier
*
* @author Javier E Hasbun 2008.
*/
package org.opensourcephysics.numerics.specialfunctions;
import org.opensourcephysics.numerics.Function;
public class Gamma implements Function {
/**
* Implements the Function interface for the gamma function .
*
* @param x
* @return gamma function at x
*/
public double evaluate(double x) {
return gamma(x);
}
public static double gamma(double x) {
/*
* 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
*/
double sgngam, q, z, y, p1, q1;
int ip, p;
double[] pp = {1.60119522476751861407E-4, 1.19135147006586384913E-3, 1.04213797561761569935E-2, 4.76367800457137231464E-2, 2.07448227648435975150E-1, 4.94214826801497100753E-1, 9.99999999999999996796E-1};
double[] qq = {-2.31581873324120129819E-5, 5.39605580493303397842E-4, -4.45641913851797240494E-3, 1.18139785222060435552E-2, 3.58236398605498653373E-2, -2.34591795718243348568E-1, 7.14304917030273074085E-2, 1.00000000000000000320};
sgngam = 1;
q = Math.abs(x);
if(q>33.0) {
if(x<0.0) {
p = (int) Math.floor(q);
ip = Math.round(p);
if(ip%2==0) {
sgngam = -1;
}
z = q-p;
if(z>0.5) {
p = p+1;
z = q-p;
}
z = q*Math.sin(Math.PI*z);
z = Math.abs(z);
z = Math.PI/(z*gammastirf(q));
} else {
z = gammastirf(x);
}
y = sgngam*z;
return y;
}
z = 1;
while(x>=3) {
x = x-1;
z = z*x;
}
while(x<0) {
if(x>-0.000000001) {
y = z/((1+0.5772156649015329*x)*x);
return y;
}
z = z/x;
x = x+1;
}
while(x<2) {
if(x<0.000000001) {
y = z/((1+0.5772156649015329*x)*x);
return y;
}
z = z/x;
x = x+1.0;
}
if(x==2) {
y = z;
return y;
}
x = x-2.0;
p1 = pp[0];
for(int i = 1; i<7; i++) {
p1 = pp[i]+p1*x;
}
q1 = qq[0];
for(int i = 1; i<8; i++) {
q1 = qq[i]+q1*x;
}
return z*p1/q1;
}
private static double gammastirf(double x) {
double p1, w, y, v;
w = 1/x;
double[] pp = {7.87311395793093628397E-4, -2.29549961613378126380E-4, -2.68132617805781232825E-3, 3.47222221605458667310E-3, 8.33333333333482257126E-2};
p1 = pp[0];
for(int i = 1; i<5; i++) {
p1 = pp[i]+p1*x;
}
w = 1+w*p1;
y = Math.exp(x);
if(x>143.01608) {
v = Math.pow(x, 0.5*x-0.25);
y = v*(v/y);
} else {
y = Math.pow(x, x-0.5)/y;
}
return 2.50662827463100050242*y*w;
}
}
/*
* Open Source Physics software is free software; you can redistribute
* it and/or modify it under the terms of the GNU General Public License (GPL) as
* published by the Free Software Foundation; either version 2 of the License,
* or(at your option) any later version.
* Code that uses any portion of the code in the org.opensourcephysics package
* or any subpackage (subdirectory) of this package must also be released
* under the GNU GPL license.
*
* This software 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; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
* or view the license online at http://www.gnu.org/copyleft/gpl.html
*
* Copyright (c) 2007 The Open Source Physics project
* http://www.opensourcephysics.org
*/