Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation

Doxygen Source Code Documentation


Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

cdf_00.c

Go to the documentation of this file.
00001 #include "cdflib.h"
00002 double algdiv(double *a,double *b)
00003 /*
00004 -----------------------------------------------------------------------
00005  
00006      COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8
00007  
00008                          --------
00009  
00010      IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
00011      LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
00012  
00013 -----------------------------------------------------------------------
00014 */
00015 {
00016 static double c0 = .833333333333333e-01;
00017 static double c1 = -.277777777760991e-02;
00018 static double c2 = .793650666825390e-03;
00019 static double c3 = -.595202931351870e-03;
00020 static double c4 = .837308034031215e-03;
00021 static double c5 = -.165322962780713e-02;
00022 static double algdiv,c,d,h,s11,s3,s5,s7,s9,t,u,v,w,x,x2,T1;
00023 /*
00024      ..
00025      .. Executable Statements ..
00026 */
00027     if(*a <= *b) goto S10;
00028     h = *b/ *a;
00029     c = 1.0e0/(1.0e0+h);
00030     x = h/(1.0e0+h);
00031     d = *a+(*b-0.5e0);
00032     goto S20;
00033 S10:
00034     h = *a/ *b;
00035     c = h/(1.0e0+h);
00036     x = 1.0e0/(1.0e0+h);
00037     d = *b+(*a-0.5e0);
00038 S20:
00039 /*
00040                 SET SN = (1 - X**N)/(1 - X)
00041 */
00042     x2 = x*x;
00043     s3 = 1.0e0+(x+x2);
00044     s5 = 1.0e0+(x+x2*s3);
00045     s7 = 1.0e0+(x+x2*s5);
00046     s9 = 1.0e0+(x+x2*s7);
00047     s11 = 1.0e0+(x+x2*s9);
00048 /*
00049                 SET W = DEL(B) - DEL(A + B)
00050 */
00051     t = pow(1.0e0/ *b,2.0);
00052     w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;
00053     w *= (c/ *b);
00054 /*
00055                     COMBINE THE RESULTS
00056 */
00057     T1 = *a/ *b;
00058     u = d*alnrel(&T1);
00059     v = *a*(log(*b)-1.0e0);
00060     if(u <= v) goto S30;
00061     algdiv = w-v-u;
00062     return algdiv;
00063 S30:
00064     algdiv = w-u-v;
00065     return algdiv;
00066 } /* END */
 

Powered by Plone

This site conforms to the following standards: