[Scipy-svn] r4468 - in trunk/scipy/special: . tests

scipy-svn@scip... scipy-svn@scip...
Mon Jun 23 19:54:24 CDT 2008


Author: wnbell
Date: 2008-06-23 19:54:19 -0500 (Mon, 23 Jun 2008)
New Revision: 4468

Modified:
   trunk/scipy/special/_cephesmodule.c
   trunk/scipy/special/specfun.pyf
   trunk/scipy/special/specfun_wrappers.c
   trunk/scipy/special/specfun_wrappers.h
   trunk/scipy/special/tests/test_basic.py
Log:
addresses ticket #659


Modified: trunk/scipy/special/_cephesmodule.c
===================================================================
--- trunk/scipy/special/_cephesmodule.c	2008-06-23 21:14:11 UTC (rev 4467)
+++ trunk/scipy/special/_cephesmodule.c	2008-06-24 00:54:19 UTC (rev 4468)
@@ -108,7 +108,7 @@
 static void * gdtri_data[] = { (void *)gdtri, (void *)gdtri, };
 */
 static void * hyp2f1_data[] = { (void *)hyp2f1, (void *)hyp2f1, (void *)chyp2f1_wrap, (void *)chyp2f1_wrap};
-static void * hyperg_data[] = { (void *)hyperg, (void *)hyperg, (void *)chyp1f1_wrap, (void *)chyp1f1_wrap};
+static void * hyp1f1_data[] = { (void *)hyp1f1_wrap, (void *)hyp1f1_wrap, (void *)chyp1f1_wrap, (void *)chyp1f1_wrap};
 static void * hypU_data[] = { (void *)hypU_wrap, (void *)hypU_wrap, };
 static void * hyp2f0_data[] = { (void *)hyp2f0, (void *)hyp2f0, };
 static void * threef0_data[] = { (void *)threef0, (void *)threef0, };
@@ -441,7 +441,7 @@
 	f = PyUFunc_FromFuncAndData(cephes4_functions, hyp2f1_data, cephes_5c2_types, 4, 4, 1, PyUFunc_None, "hyp2f1", hyp2f1_doc, 0);
 	PyDict_SetItemString(dictionary, "hyp2f1", f);
 	Py_DECREF(f);
-	f = PyUFunc_FromFuncAndData(cephes3_functions, hyperg_data, cephes_4c_types, 4, 3, 1, PyUFunc_None, "hyp1f1", hyp1f1_doc, 0);
+	f = PyUFunc_FromFuncAndData(cephes3_functions, hyp1f1_data, cephes_4c_types, 4, 3, 1, PyUFunc_None, "hyp1f1", hyp1f1_doc, 0);
 	PyDict_SetItemString(dictionary, "hyp1f1", f);
 	Py_DECREF(f);
 

Modified: trunk/scipy/special/specfun.pyf
===================================================================
--- trunk/scipy/special/specfun.pyf	2008-06-23 21:14:11 UTC (rev 4467)
+++ trunk/scipy/special/specfun.pyf	2008-06-24 00:54:19 UTC (rev 4468)
@@ -242,7 +242,12 @@
         ! eix
         ! e1xb
 
-        ! chgm
+        subroutine chgm(a,b,x,hg) ! in :specfun:specfun.f
+             double precision intent(in) :: a
+             double precision intent(in) :: b
+             double precision intent(in) :: x
+             double precision intent(out) :: hg
+        end subroutine chgm
 
         ! stvh0
 

Modified: trunk/scipy/special/specfun_wrappers.c
===================================================================
--- trunk/scipy/special/specfun_wrappers.c	2008-06-23 21:14:11 UTC (rev 4467)
+++ trunk/scipy/special/specfun_wrappers.c	2008-06-24 00:54:19 UTC (rev 4468)
@@ -29,6 +29,7 @@
 extern void F_FUNC(cpsi,CPSI)(double*,double*,double*,double*);
 extern void F_FUNC(hygfz,HYGFZ)(double*,double*,double*,Py_complex*,Py_complex*);
 extern void F_FUNC(cchg,CCHG)(double*,double*,Py_complex*,Py_complex*);
+extern void F_FUNC(chgm,CHGM)(double*,double*,double*,double*);
 extern void F_FUNC(chgu,CHGU)(double*,double*,double*,double*,int*);
 extern void F_FUNC(itairy,ITAIRY)(double*,double*,double*,double*,double*);
 extern void F_FUNC(e1xb,E1XB)(double*,double*);
@@ -147,6 +148,15 @@
   
 }
 
+double hyp1f1_wrap(double a, double b, double x) {
+   double outy;
+ 
+   F_FUNC(chgm,CHGM)(&a, &b, &x, &outy);
+   if (outy == 1e300) {
+     outy = INFINITY;
+   }
+   return outy;
+}
 
 int itairy_wrap(double x, double *apt, double *bpt, double *ant, double *bnt) {
   double tmp; 

Modified: trunk/scipy/special/specfun_wrappers.h
===================================================================
--- trunk/scipy/special/specfun_wrappers.h	2008-06-23 21:14:11 UTC (rev 4467)
+++ trunk/scipy/special/specfun_wrappers.h	2008-06-24 00:54:19 UTC (rev 4468)
@@ -31,6 +31,7 @@
 Py_complex crgamma_wrap( Py_complex z);
 Py_complex chyp2f1_wrap( double a, double b, double c, Py_complex z);
 Py_complex chyp1f1_wrap( double a, double b, Py_complex z);
+double hyp1f1_wrap( double a, double b, double x);
 double hypU_wrap(double a, double b, double x);
 double exp1_wrap(double x);
 double expi_wrap(double x);

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2008-06-23 21:14:11 UTC (rev 4467)
+++ trunk/scipy/special/tests/test_basic.py	2008-06-24 00:54:19 UTC (rev 4468)
@@ -32,7 +32,7 @@
 #8   test_sh_jacobi
 #8   test_sh_legendre
 
-from numpy import dot
+from numpy import dot, array
 
 from scipy.testing import *
 
@@ -1177,6 +1177,116 @@
         hyp1 = hyp1f1(.1,.1,.3)
         assert_almost_equal(hyp1, 1.3498588075760032,7)
 
+        # test contributed by Moritz Deger (2008-05-29)
+        # http://projects.scipy.org/scipy/scipy/ticket/659
+        
+        # reference data obtained from mathematica [ a, b, x, m(a,b,x)]:
+        # produced with test_hyp1f1.nb
+        ref_data = array([[ -8.38132975e+00,  -1.28436461e+01,  -2.91081397e+01,          1.04178330e+04],
+                          [  2.91076882e+00,  -6.35234333e+00,  -1.27083993e+01,          6.68132725e+00],
+                          [ -1.42938258e+01,   1.80869131e-01,   1.90038728e+01,          1.01385897e+05],
+                          [  5.84069088e+00,   1.33187908e+01,   2.91290106e+01,          1.59469411e+08],
+                          [ -2.70433202e+01,  -1.16274873e+01,  -2.89582384e+01,          1.39900152e+24],
+                          [  4.26344966e+00,  -2.32701773e+01,   1.91635759e+01,          6.13816915e+21],
+                          [  1.20514340e+01,  -3.40260240e+00,   7.26832235e+00,          1.17696112e+13],
+                          [  2.77372955e+01,  -1.99424687e+00,   3.61332246e+00,          3.07419615e+13],
+                          [  1.50310939e+01,  -2.91198675e+01,  -1.53581080e+01,         -3.79166033e+02],
+                          [  1.43995827e+01,   9.84311196e+00,   1.93204553e+01,          2.55836264e+10],
+                          [ -4.08759686e+00,   1.34437025e+01,  -1.42072843e+01,          1.70778449e+01],
+                          [  8.05595738e+00,  -1.31019838e+01,   1.52180721e+01,          3.06233294e+21],
+                          [  1.81815804e+01,  -1.42908793e+01,   9.57868793e+00,         -2.84771348e+20],
+                          [ -2.49671396e+01,   1.25082843e+01,  -1.71562286e+01,          2.36290426e+07],
+                          [  2.67277673e+01,   1.70315414e+01,   6.12701450e+00,          7.77917232e+03],
+                          [  2.49565476e+01,   2.91694684e+01,   6.29622660e+00,          2.35300027e+02],
+                          [  6.11924542e+00,  -1.59943768e+00,   9.57009289e+00,          1.32906326e+11],
+                          [ -1.47863653e+01,   2.41691301e+01,  -1.89981821e+01,          2.73064953e+03],
+                          [  2.24070483e+01,  -2.93647433e+00,   8.19281432e+00,         -6.42000372e+17],
+                          [  8.04042600e-01,   1.82710085e+01,  -1.97814534e+01,          5.48372441e-01],
+                          [  1.39590390e+01,   1.97318686e+01,   2.37606635e+00,          5.51923681e+00],
+                          [ -4.66640483e+00,  -2.00237930e+01,   7.40365095e+00,          4.50310752e+00],
+                          [  2.76821999e+01,  -6.36563968e+00,   1.11533984e+01,         -9.28725179e+23],
+                          [ -2.56764457e+01,   1.24544906e+00,   1.06407572e+01,          1.25922076e+01],
+                          [  3.20447808e+00,   1.30874383e+01,   2.26098014e+01,          2.03202059e+04],
+                          [ -1.24809647e+01,   4.15137113e+00,  -2.92265700e+01,          2.39621411e+08],
+                          [  2.14778108e+01,  -2.35162960e+00,  -1.13758664e+01,          4.46882152e-01],
+                          [ -9.85469168e+00,  -3.28157680e+00,   1.67447548e+01,         -1.07342390e+07],
+                          [  1.08122310e+01,  -2.47353236e+01,  -1.15622349e+01,         -2.91733796e+03],
+                          [ -2.67933347e+01,  -3.39100709e+00,   2.56006986e+01,         -5.29275382e+09],
+                          [ -8.60066776e+00,  -8.02200924e+00,   1.07231926e+01,          1.33548320e+06],
+                          [ -1.01724238e-01,  -1.18479709e+01,  -2.55407104e+01,          1.55436570e+00],
+                          [ -3.93356771e+00,   2.11106818e+01,  -2.57598485e+01,          2.13467840e+01],
+                          [  3.74750503e+00,   1.55687633e+01,  -2.92841720e+01,          1.43873509e-02],
+                          [  6.99726781e+00,   2.69855571e+01,  -1.63707771e+01,          3.08098673e-02],
+                          [ -2.31996011e+01,   3.47631054e+00,   9.75119815e-01,          1.79971073e-02],
+                          [  2.38951044e+01,  -2.91460190e+01,  -2.50774708e+00,          9.56934814e+00],
+                          [  1.52730825e+01,   5.77062507e+00,   1.21922003e+01,          1.32345307e+09],
+                          [  1.74673917e+01,   1.89723426e+01,   4.94903250e+00,          9.90859484e+01],
+                          [  1.88971241e+01,   2.86255413e+01,   5.52360109e-01,          1.44165360e+00],
+                          [  1.02002319e+01,  -1.66855152e+01,  -2.55426235e+01,          6.56481554e+02],
+                          [ -1.79474153e+01,   1.22210200e+01,  -1.84058212e+01,          8.24041812e+05],
+                          [ -1.36147103e+01,   1.32365492e+00,  -7.22375200e+00,          9.92446491e+05],
+                          [  7.57407832e+00,   2.59738234e+01,  -1.34139168e+01,          3.64037761e-02],
+                          [  2.21110169e+00,   1.28012666e+01,   1.62529102e+01,          1.33433085e+02],
+                          [ -2.64297569e+01,  -1.63176658e+01,  -1.11642006e+01,         -2.44797251e+13],
+                          [ -2.46622944e+01,  -3.02147372e+00,   8.29159315e+00,         -3.21799070e+05],
+                          [ -1.37215095e+01,  -1.96680183e+01,   2.91940118e+01,          3.21457520e+12],
+                          [ -5.45566105e+00,   2.81292086e+01,   1.72548215e-01,          9.66973000e-01],
+                          [ -1.55751298e+00,  -8.65703373e+00,   2.68622026e+01,         -3.17190834e+16],
+                          [  2.45393609e+01,  -2.70571903e+01,   1.96815505e+01,          1.80708004e+37],
+                          [  5.77482829e+00,   1.53203143e+01,   2.50534322e+01,          1.14304242e+06],
+                          [ -1.02626819e+01,   2.36887658e+01,  -2.32152102e+01,          7.28965646e+02],
+                          [ -1.30833446e+00,  -1.28310210e+01,   1.87275544e+01,         -9.33487904e+12],
+                          [  5.83024676e+00,  -1.49279672e+01,   2.44957538e+01,         -7.61083070e+27],
+                          [ -2.03130747e+01,   2.59641715e+01,  -2.06174328e+01,          4.54744859e+04],
+                          [  1.97684551e+01,  -2.21410519e+01,  -2.26728740e+01,          3.53113026e+06],
+                          [  2.73673444e+01,   2.64491725e+01,   1.57599882e+01,          1.07385118e+07],
+                          [  5.73287971e+00,   1.21111904e+01,   1.33080171e+01,          2.63220467e+03],
+                          [ -2.82751072e+01,   2.08605881e+01,   9.09838900e+00,         -6.60957033e-07],
+                          [  1.87270691e+01,  -1.74437016e+01,   1.52413599e+01,          6.59572851e+27],
+                          [  6.60681457e+00,  -2.69449855e+00,   9.78972047e+00,         -2.38587870e+12],
+                          [  1.20895561e+01,  -2.51355765e+01,   2.30096101e+01,          7.58739886e+32],
+                          [ -2.44682278e+01,   2.10673441e+01,  -1.36705538e+01,          4.54213550e+04],
+                          [ -4.50665152e+00,   3.72292059e+00,  -4.83403707e+00,          2.68938214e+01],
+                          [ -7.46540049e+00,  -1.08422222e+01,  -1.72203805e+01,         -2.09402162e+02],
+                          [ -2.00307551e+01,  -7.50604431e+00,  -2.78640020e+01,          4.15985444e+19],
+                          [  1.99890876e+01,   2.20677419e+01,  -2.51301778e+01,          1.23840297e-09],
+                          [  2.03183823e+01,  -7.66942559e+00,   2.10340070e+01,          1.46285095e+31],
+                          [ -2.90315825e+00,  -2.55785967e+01,  -9.58779316e+00,          2.65714264e-01],
+                          [  2.73960829e+01,  -1.80097203e+01,  -2.03070131e+00,          2.52908999e+02],
+                          [ -2.11708058e+01,  -2.70304032e+01,   2.48257944e+01,          3.09027527e+08],
+                          [  2.21959758e+01,   4.00258675e+00,  -1.62853977e+01,         -9.16280090e-09],
+                          [  1.61661840e+01,  -2.26845150e+01,   2.17226940e+01,         -8.24774394e+33],
+                          [ -3.35030306e+00,   1.32670581e+00,   9.39711214e+00,         -1.47303163e+01],
+                          [  7.23720726e+00,  -2.29763909e+01,   2.34709682e+01,         -9.20711735e+29],
+                          [  2.71013568e+01,   1.61951087e+01,  -7.11388906e-01,          2.98750911e-01],
+                          [  8.40057933e+00,  -7.49665220e+00,   2.95587388e+01,          6.59465635e+29],
+                          [ -1.51603423e+01,   1.94032322e+01,  -7.60044357e+00,          1.05186941e+02],
+                          [ -8.83788031e+00,  -2.72018313e+01,   1.88269907e+00,          1.81687019e+00],
+                          [ -1.87283712e+01,   5.87479570e+00,  -1.91210203e+01,          2.52235612e+08],
+                          [ -5.61338513e-01,   2.69490237e+01,   1.16660111e-01,          9.97567783e-01],
+                          [ -5.44354025e+00,  -1.26721408e+01,  -4.66831036e+00,          1.06660735e-01],
+                          [ -2.18846497e+00,   2.33299566e+01,   9.62564397e+00,          3.03842061e-01],
+                          [  6.65661299e+00,  -2.39048713e+01,   1.04191807e+01,          4.73700451e+13],
+                          [ -2.57298921e+01,  -2.60811296e+01,   2.74398110e+01,         -5.32566307e+11],
+                          [ -1.11431826e+01,  -1.59420160e+01,  -1.84880553e+01,         -1.01514747e+02],
+                          [  6.50301931e+00,   2.59859051e+01,  -2.33270137e+01,          1.22760500e-02],
+                          [ -1.94987891e+01,  -2.62123262e+01,   3.90323225e+00,          1.71658894e+01],
+                          [  7.26164601e+00,  -1.41469402e+01,   2.81499763e+01,         -2.50068329e+31],
+                          [ -1.52424040e+01,   2.99719005e+01,  -2.85753678e+01,          1.31906693e+04],
+                          [  5.24149291e+00,  -1.72807223e+01,   2.22129493e+01,          2.50748475e+25],
+                          [  3.63207230e-01,  -9.54120862e-02,  -2.83874044e+01,          9.43854939e-01],
+                          [ -2.11326457e+00,  -1.25707023e+01,   1.17172130e+00,          1.20812698e+00],
+                          [  2.48513582e+00,   1.03652647e+01,  -1.84625148e+01,          6.47910997e-02],
+                          [  2.65395942e+01,   2.74794672e+01,   1.29413428e+01,          2.89306132e+05],
+                          [ -9.49445460e+00,   1.59930921e+01,  -1.49596331e+01,          3.27574841e+02],
+                          [ -5.89173945e+00,   9.96742426e+00,   2.60318889e+01,         -3.15842908e-01],
+                          [ -1.15387239e+01,  -2.21433107e+01,  -2.17686413e+01,          1.56724718e-01],
+                          [ -5.30592244e+00,  -2.42752190e+01,   1.29734035e+00,          1.31985534e+00]])
+
+        for a,b,c,expected in ref_data:
+            result = hyp1f1(a,b,c)
+            assert(abs(expected - result)/expected < 1e-4)
+
     def test_hyp1f2(self):
         pass
 



More information about the Scipy-svn mailing list