specialfunction.h

数学特殊函数算法库

while((ep>0.0000001||FloatEqual(ep,0.0000001))&&(Abs(h)>s)) 
{  
	g=0.0; 
	for(i=1;i<=m;i++) 
	{ 
		aa=(i-1.0)*h; 
		bb=i*h; 
		w=0.0; 
		for(j=0;j<5;j++) 
		{  
			xx=((bb-aa)*t[j]+(bb+aa))/2.0; 
			q=sqrt(1.0-k*k*sin(xx)*sin(xx)); 
			w=w+q*c[j]; 
		} 
		g=g+w; 
	} 
	g=g*h/2.0; 
	ep=Abs(g-p)/(1.0+Abs(g)); 
	p=g;  
	m=m+m;  
	h=0.5*h; 
} 
return(g); 
} 
_MITC_MMATH_END
------------------------------------------------
/*  This header file contains several member functions which returns
    values of some special functions, at certain argument values.  These
    functions are used in Test.h, particularly for the classes 
    SvarMeanTest, DvarMeanTest, and VarTest.  These functions are also used
    in Generate.h, particularly for the member function poisson.
*/    


#include <math.h>
#include <iostream.h>

class SpecFunc {
public:

     SpecFunc();
     ~SpecFunc();
     float Gammln(float xx);			// Natural Log of gamma function.
     float Betai(float a,float b,float x);	// Incomplete Beta Function.
     float Sqr(float x);			// Square
     
     
private:
     int j;
     double x,y,tmp,ser;
     double* cof; 
     float bt;
     int MAXIT;
     float EPS;
     float FPMIN;
     void error(char ErrorText[]);
     int m,m2;
     float aa,c,d,h,qab,qam,qap;
     float Betacf(float a,float b,float x);   	// Function needed only for Betai.

};

SpecFunc::SpecFunc() {

     cof=new double [6];
     cof[0]=76.1800917294714;
     cof[1]=-86.50532032941677;
     cof[2]=24.01409824083091;
     cof[3]=-1.231739572450155;
     cof[4]=0.1208650973866179e-2;
     cof[5]=-0.5395239384953e-5;
     MAXIT=100;
     EPS=3.0e-7;
     FPMIN=1.0e-30;     
}     

SpecFunc::~SpecFunc() {delete [] cof;}

// Note this code for calculating these functions are in "Numerical Recipes in C
// The Art of Scientific Computing Second Edition"

float SpecFunc::Gammln(float xx) {    	// In "Numerical Recipes" Ch. 6-1, p. 214
					
     y=x=xx;
     tmp=x + 5.5;
     tmp -= (x+.5)*log(tmp);
     ser=1.000000000190015;
     for (j=0;j<=5;j++) ser += cof[j]/++y;
          	
     return -tmp + log(2.5066282746210005*ser/x);
}     

float SpecFunc::Betai(float a,float b,float x) {    // In "Numerical Recipes" Ch. 6-4, p. 227

     if(x < 0.0 || x > 1.0) error("Bad x in routine betai");

     if(x == 0.0 || x == 1.0) {bt=0.0;}

     else {bt=exp(Gammln(a+b) - Gammln(a) - Gammln(b)
           + a*log(x) + b*log(1.0 - x));}
     
     if(x < (a+1.0)/(a+b+2.0)) {return bt*Betacf(a,b,x)/a;}
     
     else {return 1.0 - bt*Betacf(b,a,1.0 - x)/b;}
}

float SpecFunc::Sqr(float x) {return x*x;}

float SpecFunc::Betacf(float a,float b,float x) {   // In "Numerical Recipes" Ch. 6-4, p. 227 - 228
					    // It is used for the calculation of the
     qab=a+b;				    // above function Betai only.
     qap=a+1.0;
     qam=a-1.0;
     c=1.0;
     d=1.0-qab*x/qap;
     
     if(fabs(d) < FPMIN) d=FPMIN;
     
     d=1.0/d;
     h=d;
     for (m=1;m<=MAXIT;m++) {
     
          m2=2*m;
          aa=m*(b-m)*x/((qam+m2)*(a+m2));          
          d=1.0 + aa*d;
          if (fabs(d) < FPMIN) d=FPMIN;
          c=1.0 + aa/c;
          if (fabs(c) < FPMIN) c=FPMIN;
          d=1.0/d;
          h *= d*c;
          aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
          d=1.0+aa*d;
          if (fabs(d) < FPMIN) d=FPMIN;
          c=1.0 + aa/c;
          if (fabs(c) < FPMIN) c=FPMIN;
          d=1.0/d;
          h *= d*c;
          if (fabs(d*c - 1.0) < EPS) break;
     }
      
     if (m > MAXIT) error("a or b too big, or MAXIT too small in betacf");
      
     return h;
} 

void SpecFunc::error(char ErrorText[]) {       

     cerr<<endl<<ErrorText<<endl;
     cout.flush();
     abort();
--------------------------------------------------
//IncompleteGammaFunction.cpp           
//不完全伽马函数    
   
#include <IOSTREAM>             //模板类输入输出流标准头文件    
#include <SPECIALFUNCTION.H>    //数学变换头文件    
using namespace std;            //名字空间    
   
void main(void)   
{      
    double a[3] = {0.5,5.0,50.0};   
    double x[3] = {0.1,1.0,10.0};   
       
    for(int i=0; i<3; i++)   
    for(int j=0; j<3; j++)   
    {   
        double s = a[i];   
        double t = x[j];   
        double y = IncompleteGammaFunction(s, t);   
        cout << " GammaFunction(" << a[i] << "," << x[j]    
                                            << ")\t = " << y << endl;   
    }   
    cout << endl;   
}   
-=------------------
 1    /* 
 2     * Licensed to the Apache Software Foundation (ASF) under one or more 
 3     * contributor license agreements.  See the NOTICE file distributed with 
 4     * this work for additional information regarding copyright ownership. 
 5     * The ASF licenses this file to You under the Apache License, Version 2.0 
 6     * (the "License"); you may not use this file except in compliance with 
 7     * the License.  You may obtain a copy of the License at 
 8     * 
 9     *      http://www.apache.org/licenses/LICENSE-2.0 
 10     * 
 11     * Unless required by applicable law or agreed to in writing, software 
 12     * distributed under the License is distributed on an "AS IS" BASIS, 
 13     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 
 14     * See the License for the specific language governing permissions and 
 15     * limitations under the License. 
 16     */ 
 17    package org.apache.commons.math.special; 
 18     
 19    import java.io.Serializable; 
 20     
 21    import org.apache.commons.math.MathException; 
 22    import org.apache.commons.math.MaxIterationsExceededException; 
 23    import org.apache.commons.math.util.ContinuedFraction; 
 24     
 25    /** 
 26     * This is a utility class that provides computation methods related to the 
 27     * Gamma family of functions. 
 28     * 
 29     * @version $Revision: 762121 $ $Date: 2009-04-05 19:10:23 +0200 (dim., 05 avr. 2009) $ 
 30     */ 
 31    public class Gamma implements Serializable { 
 32         
 33        /** Serializable version identifier */ 
 34        private static final long serialVersionUID = -6587513359895466954L; 
 35     
 36        /** Maximum allowed numerical error. */ 
 37        private static final double DEFAULT_EPSILON = 10e-15; 
 38     
 39        /** Lanczos coefficients */ 
 40  1      private static final double[] lanczos = 
 41        { 
 42            0.99999999999999709182, 
 43            57.156235665862923517, 
 44            -59.597960355475491248, 
 45            14.136097974741747174, 
 46            -0.49191381609762019978, 
 47            .33994649984811888699e-4, 
 48            .46523628927048575665e-4, 
 49            -.98374475304879564677e-4, 
 50            .15808870322491248884e-3, 
 51            -.21026444172410488319e-3, 
 52            .21743961811521264320e-3, 
 53            -.16431810653676389022e-3, 
 54            .84418223983852743293e-4, 
 55            -.26190838401581408670e-4, 
 56            .36899182659531622704e-5, 
 57        }; 
 58     
 59        /** Avoid repeated computation of log of 2 PI in logGamma */ 
 60  1      private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI); 
 61     
 62         
 63        /** 
 64         * Default constructor.  Prohibit instantiation. 
 65         */ 
 66        private Gamma() { 
 67  0          super(); 
 68  0      } 
 69     
 70        /** 
 71         * Returns the natural logarithm of the gamma function &#915;(x). 
 72         * 
 73         * The implementation of this method is based on: 
 74         * <ul> 
 75         * <li><a href="http://mathworld.wolfram.com/GammaFunction.html"> 
 76         * Gamma Function</a>, equation (28).</li> 
 77         * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html"> 
 78         * Lanczos Approximation</a>, equations (1) through (5).</li> 
 79         * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on 
 80         * the computation of the convergent Lanczos complex Gamma approximation 
 81         * </a></li> 
 82         * </ul> 
 83         *  
 84         * @param x the value. 
 85         * @return log(&#915;(x)) 
 86         */ 
 87        public static double logGamma(double x) { 
 88            double ret; 
 89     
 90  19968          if (Double.isNaN(x) || (x <= 0.0)) { 
 91  3              ret = Double.NaN; 
 92            } else { 
 93  19965              double g = 607.0 / 128.0; 
 94                 
 95  19965              double sum = 0.0; 
 96  299475              for (int i = lanczos.length - 1; i > 0; --i) { 
 97  279510                  sum = sum + (lanczos[i] / (x + i)); 
 98                } 
 99  19965              sum = sum + lanczos[0]; 
 100     
 101  19965              double tmp = x + g + .5; 
 102  19965              ret = ((x + .5) * Math.log(tmp)) - tmp + 
 103                    HALF_LOG_2_PI + Math.log(sum / x); 
 104            } 
 105     
 106  19968          return ret; 
 107        } 
 108     
 109        /** 
 110         * Returns the regularized gamma function P(a, x). 
 111         *  
 112         * @param a the a parameter. 
 113         * @param x the value. 
 114         * @return the regularized gamma function P(a, x) 
 115         * @throws MathException if the algorithm fails to converge. 
 116         */ 
 117        public static double regularizedGammaP(double a, double x) 
 118            throws MathException 
 119        { 
 120  806          return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); 
 121        } 
 122             
 123             
 124        /** 
 125         * Returns the regularized gamma function P(a, x). 
 126         *  
 127         * The implementation of this method is based on: 
 128         * <ul> 
 129         * <li> 
 130         * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> 
 131         * Regularized Gamma Function</a>, equation (1).</li> 
 132         * <li> 
 133         * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html"> 
 134         * Incomplete Gamma Function</a>, equation (4).</li> 
 135         * <li> 
 136         * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html"> 
 137         * Confluent Hypergeometric Function of the First Kind</a>, equation (1). 
 138         * </li> 
 139         * </ul> 
 140         *  
 141         * @param a the a parameter. 
 142         * @param x the value. 
 143         * @param epsilon When the absolute value of the nth item in the 
 144         *                series is less than epsilon the approximation ceases 
 145         *                to calculate further elements in the series. 
 146         * @param maxIterations Maximum number of "iterations" to complete.  
 147         * @return the regularized gamma function P(a, x) 
 148         * @throws MathException if the algorithm fails to converge. 
 149         */ 
 150        public static double regularizedGammaP(double a,  
 151                                               double x,  
 152                                               double epsilon,  
 153                                               int maxIterations)  
 154            throws MathException 
 155        { 
 156            double ret; 
 157     
 158  3790          if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { 
 159  5              ret = Double.NaN; 
 160  3785          } else if (x == 0.0) { 
 161  62              ret = 0.0; 
 162  3723          } else if (a >= 1.0 && x > a) { 
 163                // use regularizedGammaQ because it should converge faster in this 
 164                // case. 
 165  308              ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations); 
 166            } else { 
 167                // calculate series 
 168  3415              double n = 0.0; // current element index 
 169  3415              double an = 1.0 / a; // n-th element in the series 
 170  3415              double sum = an; // partial sum 
 171  3061399              while (Math.abs(an) > epsilon && n < maxIterations) { 
 172                    // compute next element in the series 
 173  3057984                  n = n + 1.0; 
 174  3057984                  an = an * (x / (a + n)); 
 175     
 176                    // update partial sum 
 177  3057984                  sum = sum + an; 
 178                } 
 179  3415              if (n >= maxIterations) { 
 180  24                  throw new MaxIterationsExceededException(maxIterations); 
 181                } else { 
 182  3391                  ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * sum; 
 183                } 
 184            } 
 185     
 186  3766          return ret; 
 187        } 
 188         
 189        /** 
 190         * Returns the regularized gamma function Q(a, x) = 1 - P(a, x). 
 191         *  
 192         * @param a the a parameter. 
 193         * @param x the value. 
 194         * @return the regularized gamma function Q(a, x) 
 195         * @throws MathException if the algorithm fails to converge. 
 196         */ 
 197        public static double regularizedGammaQ(double a, double x) 
 198            throws MathException 
 199        { 
 200  7          return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); 
 201        } 
 202         
 203        /** 
 204         * Returns the regularized gamma function Q(a, x) = 1 - P(a, x). 
 205         *  
 206         * The implementation of this method is based on: 
 207         * <ul> 
 208         * <li> 
 209         * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> 
 210         * Regularized Gamma Function</a>, equation (1).</li> 
 211         * <li> 
 212         * <a href="    http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/"> 
 213         * Regularized incomplete gamma function: Continued fraction representations  (formula 06.08.10.0003)</a></li> 
 214         * </ul> 
 215         *  
 216         * @param a the a parameter. 
 217         * @param x the value. 
 218         * @param epsilon When the absolute value of the nth item in the 
 219         *                series is less than epsilon the approximation ceases 
 220         *                to calculate further elements in the series. 
 221         * @param maxIterations Maximum number of "iterations" to complete.  
 222         * @return the regularized gamma function P(a, x) 
 223         * @throws MathException if the algorithm fails to converge. 
 224         */ 
 225        public static double regularizedGammaQ(final double a,  
 226                                               double x,  
 227                                               double epsilon,  
 228                                               int maxIterations)  
 229            throws MathException 
 230        { 
 231            double ret; 
 232     
 233  3568          if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { 
 234  5              ret = Double.NaN; 
 235  3563          } else if (x == 0.0) { 
 236  1              ret = 1.0; 
 237  3562          } else if (x < a || a < 1.0) { 
 238                // use regularizedGammaP because it should converge faster in this 
 239                // case. 
 240  2709              ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations); 
 241            } else { 
 242                // create continued fraction 
 243  853              ContinuedFraction cf = new ContinuedFraction() { 
 244     
 245                    private static final long serialVersionUID = 5378525034886164398L; 
 246     
 247                    @Override 
 248                    protected double getA(int n, double x) { 
 249  20690                      return ((2.0 * n) + 1.0) - a + x; 
 250                    } 
 251     
 252                    @Override 
 253                    protected double getB(int n, double x) { 
 254  19837                      return n * (a - n); 
 255                    } 
 256                }; 
 257                 
 258  853              ret = 1.0 / cf.evaluate(x, epsilon, maxIterations); 
 259  853              ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * ret; 
 260            } 
 261     
 262  3568          return ret; 
 263        } 
 264    } 

Report generated by Cobertura 1.