Started adding cdflib routines to cephesmodule for stats computations.

git-svn-id: http://svn.scipy.org/svn/scipy/trunk@364 d6536bca-fef9-0310-8506-e4c0a848fbcf

Lib/special/amos_wrappers.c

@@ -7,7 +7,21 @@
 
 #include "amos_wrappers.h"
 
-/* This must be linked with g77
+#if defined(NO_APPEND_FORTRAN)
+#if defined(UPPERCASE_FORTRAN)
+#define F_FUNC(f,F) F
+#else
+#define F_FUNC(f,F) f
+#endif
+#else
+#if defined(UPPERCASE_FORTRAN)
+#define F_FUNC(f,F) F##_
+#else
+#define F_FUNC(f,F) f##_
+#endif
+#endif
+
+/* This must be linked with fortran
  */
 
 int ierr_to_mtherr( int nz, int ierr) {
@@ -36,15 +50,15 @@ int cairy_wrap(Py_complex z, Py_complex *ai, Py_complex *aip, Py_complex *bi, Py
   int kode = 1;
   int nz;
 
-  zairy_(CADDR(z), &id, &kode, F2C_CST(ai), &nz, &ierr);
+  F_FUNC(zairy,ZAIRY)(CADDR(z), &id, &kode, F2C_CST(ai), &nz, &ierr);
   DO_MTHERR("airy:");
-  zbiry_(CADDR(z), &id, &kode, F2C_CST(bi), &nz, &ierr);
+  F_FUNC(zbiry,ZBIRY)(CADDR(z), &id, &kode, F2C_CST(bi), &nz, &ierr);
   DO_MTHERR("airy:");
 
   id = 1;
-  zairy_(CADDR(z), &id, &kode, F2C_CST(aip), &nz, &ierr);
+  F_FUNC(zairy,ZAIRY)(CADDR(z), &id, &kode, F2C_CST(aip), &nz, &ierr);
   DO_MTHERR("airy:");
-  zbiry_(CADDR(z), &id, &kode, F2C_CST(bip), &nz, &ierr);
+  F_FUNC(zbiry,ZBIRY)(CADDR(z), &id, &kode, F2C_CST(bip), &nz, &ierr);
   DO_MTHERR("airy:");
   return 0;
 }

Lib/special/cdf_wrappers.c

@@ -0,0 +1,57 @@
+/* This file is a collection (more can be added) of wrappers around some
+ *  CDF Fortran algorithms, so that they can be called from
+ *  cephesmodule.so
+ */
+
+#include "cdf_wrappers.h"
+#if defined(NO_APPEND_FORTRAN)
+#if defined(UPPERCASE_FORTRAN)
+#define F_FUNC(f,F) F
+#else
+#define F_FUNC(f,F) f
+#endif
+#else
+#if defined(UPPERCASE_FORTRAN)
+#define F_FUNC(f,F) F##_
+#else
+#define F_FUNC(f,F) f##_
+#endif
+#endif
+
+/* This must be linked with fortran
+ */
+
+
+int status_to_mtherr( int status) {
+     /* Return mtherr equivalents for ierr values */
+
+  if (nz != 0) return UNDERFLOW;
+
+  switch (ierr) {
+  case 1:
+    return DOMAIN;
+  case 2:
+    return OVERFLOW;
+  case 3:
+    return PLOSS;
+  case 4:
+    return TLOSS;
+  case 5:   /* Algorithm termination condition not met */
+    return TLOSS;    
+  }
+  return -1;
+}
+
+Py_complex cwofz_wrap( Py_complex z) {
+  int errflag;
+  Py_complex cy;
+
+  F_FUNC(wofz,WOFZ)(CADDR(z), CADDR(cy), &errflag);
+  if (errflag==1) mtherr("wofz:",3); /* wofz returns a single flag both
+                                        for real overflows and for domain
+                                        errors -- internal overflows from too
+                                        large abs(z)*/
+  return cy;
+}
+
+

Lib/special/cdf_wrappers.h

@@ -0,0 +1,54 @@
+/* This file is a collection of wrappers around the
+ *  Amos Fortran library of functions that take complex
+ *  variables (see www.netlib.org) so that they can be called from
+ *  the cephes library of corresponding name but work with complex
+ *  arguments.
+ */
+
+#ifndef _CDF_WRAPPERS_H
+#define _CDF_WRAPPERS_H
+#ifndef _AMOS_WRAPPERS_H
+#include "Python.h"
+#endif _AMOS_WRAPPERS_H
+extern double cdfbet3_wrap(double p, double x, double b);
+extern double cdfbet4_wrap(double p, double x, double a);
+
+extern double cdfbin2_wrap(double p, double xn, double pr);
+extern double cdfbin3_wrap(double p, double s, double pr);
+
+extern double cdfchi3_wrap(double p, double x);
+
+extern double cdfchn1_wrap(double x, double df, double nc);
+extern double cdfchn2_wrap(double p, double df, double nc);
+extern double cdfchn3_wrap(double p, double x, double nc);
+extern double cdfchn4_wrap(double p, double x, double df);
+
+extern double cdff3_wrap(double p, double f, double dfd);
+extern double cdff4_wrap(double p, double f, double dfn);
+
+extern double cdffnc1_wrap(double f, double dfn, double dfd, double nc);
+extern double cdffnc2_wrap(double p, double dfn, double dfd, double nc);
+extern double cdffnc3_wrap(double p, double f, double dfd, double nc);
+extern double cdffnc4_wrap(double p, double f, double dfn, double nc);
+extern double cdffnc5_wrap(double p, double f, double dfn, double dfd);
+
+extern double cdfgam3_wrap(double p, double x, double scl);
+extern double cdfgam4_wrap(double p, double x, double shp);
+
+
+extern double cdfnbn2_wrap(double p, double xn, double pr);
+extern double cdfnbn3_wrap(double p, double s, double pr);
+
+extern double cdfnor3_wrap(double p, double x, double std);
+extern double cdfnor4_wrap(double p, double x, double mn);
+
+extern double cdfpoi2_wrap(double p, double xlam);
+
+extern double cdft3_wrap(double p, double t);
+
+extern double cdftnc1_wrap(double t, double df, double nc);
+extern double cdftnc2_wrap(double p, double df, double nc);
+extern double cdftnc3_wrap(double p, double t, double nc);
+extern double cdftnc4_wrap(double p, double t, double df);
+
+#endif

Lib/special/cdflib/algdiv.f

@@ -0,0 +1,71 @@
+      DOUBLE PRECISION FUNCTION algdiv(a,b)
+C-----------------------------------------------------------------------
+C
+C     COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8
+C
+C                         --------
+C
+C     IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
+C     LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
+C
+C-----------------------------------------------------------------------
+C     .. Scalar Arguments ..
+      DOUBLE PRECISION a,b
+C     ..
+C     .. Local Scalars ..
+      DOUBLE PRECISION c,c0,c1,c2,c3,c4,c5,d,h,s11,s3,s5,s7,s9,t,u,v,w,
+     +                 x,x2
+C     ..
+C     .. External Functions ..
+      DOUBLE PRECISION alnrel
+      EXTERNAL alnrel
+C     ..
+C     .. Intrinsic Functions ..
+      INTRINSIC dlog
+C     ..
+C     .. Data statements ..
+      DATA c0/.833333333333333D-01/,c1/-.277777777760991D-02/,
+     +     c2/.793650666825390D-03/,c3/-.595202931351870D-03/,
+     +     c4/.837308034031215D-03/,c5/-.165322962780713D-02/
+C     ..
+C     .. Executable Statements ..
+C------------------------
+      IF (a.LE.b) GO TO 10
+      h = b/a
+      c = 1.0D0/ (1.0D0+h)
+      x = h/ (1.0D0+h)
+      d = a + (b-0.5D0)
+      GO TO 20
+
+   10 h = a/b
+      c = h/ (1.0D0+h)
+      x = 1.0D0/ (1.0D0+h)
+      d = b + (a-0.5D0)
+C
+C                SET SN = (1 - X**N)/(1 - X)
+C
+   20 x2 = x*x
+      s3 = 1.0D0 + (x+x2)
+      s5 = 1.0D0 + (x+x2*s3)
+      s7 = 1.0D0 + (x+x2*s5)
+      s9 = 1.0D0 + (x+x2*s7)
+      s11 = 1.0D0 + (x+x2*s9)
+C
+C                SET W = DEL(B) - DEL(A + B)
+C
+      t = (1.0D0/b)**2
+      w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t + c0
+      w = w* (c/b)
+C
+C                    COMBINE THE RESULTS
+C
+      u = d*alnrel(a/b)
+      v = a* (dlog(b)-1.0D0)
+      IF (u.LE.v) GO TO 30
+      algdiv = (w-v) - u
+      RETURN
+
+   30 algdiv = (w-u) - v
+      RETURN
+
+      END

Lib/special/cdflib/alngam.f

@@ -0,0 +1,128 @@
+      DOUBLE PRECISION FUNCTION alngam(x)
+C**********************************************************************
+C
+C     DOUBLE PRECISION FUNCTION ALNGAM(X)
+C                 double precision LN of the GAMma function
+C
+C
+C                              Function
+C
+C
+C     Returns the natural logarithm of GAMMA(X).
+C
+C
+C                              Arguments
+C
+C
+C     X --> value at which scaled log gamma is to be returned
+C                    X is DOUBLE PRECISION
+C
+C
+C                              Method
+C
+C
+C     If X .le. 6.0, then use recursion to get X below 3
+C     then apply rational approximation number 5236 of
+C     Hart et al, Computer Approximations, John Wiley and
+C     Sons, NY, 1968.
+C
+C     If X .gt. 6.0, then use recursion to get X to at least 12 and
+C     then use formula 5423 of the same source.
+C
+C**********************************************************************
+C
+C     .. Parameters ..
+      DOUBLE PRECISION hln2pi
+      PARAMETER (hln2pi=0.91893853320467274178D0)
+C     ..
+C     .. Scalar Arguments ..
+      DOUBLE PRECISION x
+C     ..
+C     .. Local Scalars ..
+      DOUBLE PRECISION offset,prod,xx
+      INTEGER i,n
+C     ..
+C     .. Local Arrays ..
+      DOUBLE PRECISION coef(5),scoefd(4),scoefn(9)
+C     ..
+C     .. External Functions ..
+      DOUBLE PRECISION devlpl
+      EXTERNAL devlpl
+C     ..
+C     .. Intrinsic Functions ..
+      INTRINSIC log,dble,int
+C     ..
+C     .. Data statements ..
+      DATA scoefn(1)/0.62003838007127258804D2/,
+     +     scoefn(2)/0.36036772530024836321D2/,
+     +     scoefn(3)/0.20782472531792126786D2/,
+     +     scoefn(4)/0.6338067999387272343D1/,
+     +     scoefn(5)/0.215994312846059073D1/,
+     +     scoefn(6)/0.3980671310203570498D0/,
+     +     scoefn(7)/0.1093115956710439502D0/,
+     +     scoefn(8)/0.92381945590275995D-2/,
+     +     scoefn(9)/0.29737866448101651D-2/
+      DATA scoefd(1)/0.62003838007126989331D2/,
+     +     scoefd(2)/0.9822521104713994894D1/,
+     +     scoefd(3)/-0.8906016659497461257D1/,
+     +     scoefd(4)/0.1000000000000000000D1/
+      DATA coef(1)/0.83333333333333023564D-1/,
+     +     coef(2)/-0.27777777768818808D-2/,
+     +     coef(3)/0.79365006754279D-3/,coef(4)/-0.594997310889D-3/,
+     +     coef(5)/0.8065880899D-3/
+C     ..
+C     .. Executable Statements ..
+      IF (.NOT. (x.LE.6.0D0)) GO TO 70
+      prod = 1.0D0
+      xx = x
+      IF (.NOT. (x.GT.3.0D0)) GO TO 30
+   10 IF (.NOT. (xx.GT.3.0D0)) GO TO 20
+      xx = xx - 1.0D0
+      prod = prod*xx
+      GO TO 10
+
+   20 CONTINUE
+   30 IF (.NOT. (x.LT.2.0D0)) GO TO 60
+   40 IF (.NOT. (xx.LT.2.0D0)) GO TO 50
+      prod = prod/xx
+      xx = xx + 1.0D0
+      GO TO 40
+
+   50 CONTINUE
+   60 alngam = devlpl(scoefn,9,xx-2.0D0)/devlpl(scoefd,4,xx-2.0D0)
+C
+C
+C     COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)
+C
+C
+      alngam = alngam*prod
+      alngam = log(alngam)
+      GO TO 110
+
+   70 offset = hln2pi
+C
+C
+C     IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET
+C
+C
+      n = int(12.0D0-x)
+      IF (.NOT. (n.GT.0)) GO TO 90
+      prod = 1.0D0
+      DO 80,i = 1,n
+          prod = prod* (x+dble(i-1))
+   80 CONTINUE
+      offset = offset - log(prod)
+      xx = x + dble(n)
+      GO TO 100
+
+   90 xx = x
+C
+C
+C     COMPUTE POWER SERIES
+C
+C
+  100 alngam = devlpl(coef,5,1.0D0/ (xx**2))/xx
+      alngam = alngam + offset + (xx-0.5D0)*log(xx) - xx
+  110 RETURN
+
+      END

Lib/special/cdflib/alnrel.f

@@ -0,0 +1,33 @@
+      DOUBLE PRECISION FUNCTION alnrel(a)
+C-----------------------------------------------------------------------
+C            EVALUATION OF THE FUNCTION LN(1 + A)
+C-----------------------------------------------------------------------
+C     .. Scalar Arguments ..
+      DOUBLE PRECISION a
+C     ..
+C     .. Local Scalars ..
+      DOUBLE PRECISION p1,p2,p3,q1,q2,q3,t,t2,w,x
+C     ..
+C     .. Intrinsic Functions ..
+      INTRINSIC abs,dble,dlog
+C     ..
+C     .. Data statements ..
+      DATA p1/-.129418923021993D+01/,p2/.405303492862024D+00/,
+     +     p3/-.178874546012214D-01/
+      DATA q1/-.162752256355323D+01/,q2/.747811014037616D+00/,
+     +     q3/-.845104217945565D-01/
+C     ..
+C     .. Executable Statements ..
+C--------------------------
+      IF (abs(a).GT.0.375D0) GO TO 10
+      t = a/ (a+2.0D0)
+      t2 = t*t
+      w = (((p3*t2+p2)*t2+p1)*t2+1.0D0)/ (((q3*t2+q2)*t2+q1)*t2+1.0D0)
+      alnrel = 2.0D0*t*w
+      RETURN
+C
+   10 x = 1.D0 + dble(a)
+      alnrel = dlog(x)
+      RETURN
+
+      END