Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

BUG: adding Buhring's analytic continuation method to special.hyp2f1; #8054 #8151

Closed
wants to merge 12 commits into from
Closed

BUG: adding Buhring's analytic continuation method to special.hyp2f1; #8054 #8151

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
db3b3bb
BUG: fixed hyp2f1 in small region; closes #8054
FormerPhysicist Nov 12, 2017
d41284e
DOC: updating hyp2f1 docstring
FormerPhysicist Nov 12, 2017
02db8bf
TST: added new test script for hyp2f1
FormerPhysicist Nov 12, 2017
5cfdcbf
DOC: adding myself to THANKS.txt
FormerPhysicist Nov 12, 2017
8214308
TST: changing where random.seed is set
FormerPhysicist Nov 19, 2017
c6dd2cc
Pulling master branch in
FormerPhysicist Jan 12, 2020
38bdba8
Merge remote-tracking branch 'upstream/master' into carry_buhring
h-vetinari Jan 7, 2021
2f76eab
remove trailing whitespace in affected files
h-vetinari Dec 30, 2020
0333b13
fix references in docstring
h-vetinari Dec 30, 2020
902aa2b
black scipy/special/tests/test_hyp2f1.py
h-vetinari Dec 30, 2020
3890845
blacken comments
h-vetinari Dec 30, 2020
60805b0
remove deprecated dtypes / constructor
h-vetinari Dec 30, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
32 changes: 20 additions & 12 deletions scipy/special/add_newdocs.py
View file
Open in desktop
Expand Up @@ -354,7 +354,7 @@ def add_newdoc(name, doc):
Returns
-------
eAi, eAip, eBi, eBip : array_like
Exponentially scaled Airy functions eAi and eBi, and their derivatives
Exponentially scaled Airy functions eAi and eBi, and their derivatives
eAip and eBip

Notes
Expand All @@ -370,11 +370,11 @@ def add_newdoc(name, doc):
.. [1] Donald E. Amos, "AMOS, A Portable Package for Bessel Functions
of a Complex Argument and Nonnegative Order",
http://netlib.org/amos/

Examples
--------
We can compute exponentially scaled Airy functions and their derivatives:

>>> from scipy.special import airye
>>> import matplotlib.pyplot as plt
>>> z = np.linspace(0, 50, 500)
Expand All @@ -386,9 +386,9 @@ def add_newdoc(name, doc):
... ax[ind].legend(data[2])
... ax[ind].grid(True)
>>> plt.show()

We can compute these using usual non-scaled Airy functions by:

>>> from scipy.special import airy
>>> Ai, Aip, Bi, Bip = airy(z)
>>> np.allclose(eAi, Ai * np.exp(2.0 / 3.0 * z * np.sqrt(z)))
Expand All @@ -399,16 +399,16 @@ def add_newdoc(name, doc):
True
>>> np.allclose(eBip, Bip * np.exp(-abs(np.real(2.0 / 3.0 * z * np.sqrt(z)))))
True
Comparing non-scaled and exponentially scaled ones, the usual non-scaled

Comparing non-scaled and exponentially scaled ones, the usual non-scaled
function quickly underflows for large values, whereas the exponentially
scaled function does not.

>>> airy(200)
(0.0, 0.0, nan, nan)
>>> airye(200)
(0.07501041684381093, -1.0609012305109042, 0.15003188417418148, 2.1215836725571093)

""")

add_newdoc("bdtr",
Expand Down Expand Up @@ -5040,15 +5040,24 @@ def add_newdoc(name, doc):
Here :math:`(\cdot)_n` is the Pochhammer symbol; see `poch`. When
:math:`n` is an integer the result is a polynomial of degree :math:`n`.

The implementation for complex values of ``z`` is described in [2]_.
For most complex values of ``z``, the implementation follows the method described in [1]_.
However, in the region given by

.. math::

\Re(z)>0, \ 0.9<|z|<1.1, \ |z-1|>0.75,

where the transformation methods used in [1]_ are not effective, the analytic continuation series given in [2]_ is used to evaluate the hypergeometric function.

References
----------
.. [1] NIST Digital Library of Mathematical Functions
https://dlmf.nist.gov/15.2
.. [2] S. Zhang and J.M. Jin, "Computation of Special Functions", Wiley 1996
.. [2] W. Buhring, "An Analytic Continuation of the Hypergeometric Series",
SIAM J. Math. Anal., Vol. 18, No. 3, May 1987
.. [3] Cephes Mathematical Functions Library,
http://www.netlib.org/cephes/
.. [4] S. Zhang and J.M. Jin, "Computation of Special Functions", Wiley 1996

Examples
--------
Expand Down Expand Up @@ -5094,7 +5103,6 @@ def add_newdoc(name, doc):
0.9272952180016117
>>> np.arctan(z) / z
0.9272952180016123

""")

add_newdoc("hyperu",
Expand Down
172 changes: 144 additions & 28 deletions scipy/special/specfun/specfun.f
View file
Open in desktop
Expand Up @@ -5769,7 +5769,10 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Y)
IMPLICIT COMPLEX *16 (Z)
LOGICAL L0,L1,L2,L3,L4,L5,L6
DIMENSION ZEA(3), ZEB(3), ZFA(3), ZFB(3)
LOGICAL L0,L1,L2,L3,L4,L5,L6,L8
PARAMETER (PI=3.141592653589793D0)
PARAMETER (EL=.5772156649015329D0)
X=DBLE(Z)
Y=DIMAG(Z)
EPS=1.0D-15
Expand All @@ -5785,8 +5788,6 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
BB=B
A0=CDABS(Z)
IF (A0.GT.0.95D0) EPS=1.0D-8
PI=3.141592653589793D0
EL=.5772156649015329D0
IF (L0.OR.L1) THEN
ISFER=3
RETURN
Expand Down Expand Up @@ -5823,7 +5824,7 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
ZR=ZR*(C-A+K-1.0D0)*(C-B+K-1.0D0)/(K*(C+K-1.0D0))*Z
15 ZHF=ZHF+ZR
ZHF=(1.0D0-Z)**(C-A-B)*ZHF
ELSE IF (A0.LE.1.0D0) THEN
ELSE IF (A0.LE.1.1D0) THEN
IF (X.LT.0.0D0) THEN
Z1=Z/(Z-1.0D0)
IF (C.GT.A.AND.B.LT.A.AND.B.GT.0.0) THEN
Expand All @@ -5839,11 +5840,26 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
IF (CDABS(ZHF-ZW).LT.CDABS(ZHF)*EPS) GO TO 25
20 ZW=ZHF
25 ZHF=ZC0*ZHF
ELSE IF (A0.GE.0.90D0) THEN
ELSE IF (A0.LT.0.9D0) THEN
Z00=(1.0D0,0.0D0)
IF (C-A.LT.A.AND.C-B.LT.B) THEN
Z00=(1.0D0-Z)**(C-A-B)
A=C-A
B=C-B
ENDIF
ZHF=(1.0D0,0.D0)
ZR=(1.0D0,0.0D0)
DO 100 K=1,1500
ZR=ZR*(A+K-1.0D0)*(B+K-1.0D0)/(K*(C+K-1.0D0))*Z
ZHF=ZHF+ZR
IF (CDABS(ZHF-ZW).LE.CDABS(ZHF)*EPS) GO TO 105
100 ZW=ZHF
105 ZHF=Z00*ZHF
ELSE IF (CDABS(Z-1.0D0).LT.0.75D0) THEN
GM=0.0D0
MCAB=INT(C-A-B+EPS*DSIGN(1.0D0,C-A-B))
IF (DABS(C-A-B-MCAB).LT.EPS) THEN
M=INT(C-A-B)
M=INT(C-A-B+EPS*DSIGN(1.0D0,C-A-B))
L8=(X.GT.1.0D0).AND.(Y.EQ.0.0D0)
IF (DABS(C-A-B-M).LT.EPS) THEN
CALL GAMMA2(A,GA)
CALL GAMMA2(B,GB)
CALL GAMMA2(C,GC)
Expand All @@ -5860,8 +5876,10 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
ZF0=(1.0D0,0.0D0)
ZR0=(1.0D0,0.0D0)
ZR1=(1.0D0,0.0D0)
SP0=0.D0
SP0=0.0D0
SP=0.0D0
ZLG=CDLOG(1.0D0-Z)
IF (L8) ZLG=ZLG-2.0D0*PI*(0.0D0,1.0D0)
IF (M.GE.0) THEN
ZC0=GM*GC/(GAM*GBM)
ZC1=-GC*(Z-1.0D0)**M/(GA*GB*RM)
Expand All @@ -5870,7 +5888,7 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
40 ZF0=ZF0+ZR0
DO 45 K=1,M
45 SP0=SP0+1.0D0/(A+K-1.0D0)+1.0/(B+K-1.0D0)-1.D0/K
ZF1=PA+PB+SP0+2.0D0*EL+CDLOG(1.0D0-Z)
ZF1=PA+PB+SP0+2.0D0*EL+ZLG
DO 55 K=1,500
SP=SP+(1.0D0-A)/(K*(A+K-1.0D0))+(1.0D0-B)/
& (K*(B+K-1.0D0))
Expand All @@ -5879,7 +5897,7 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
SM=SM+(1.0D0-A)/((J+K)*(A+J+K-1.0D0))
& +1.0D0/(B+J+K-1.0D0)
50 CONTINUE
ZP=PA+PB+2.0D0*EL+SP+SM+CDLOG(1.0D0-Z)
ZP=PA+PB+2.0D0*EL+SP+SM+ZLG
ZR1=ZR1*(A+M+K-1.0D0)*(B+M+K-1.0D0)/(K*(M+K))
& *(1.0D0-Z)
ZF1=ZF1+ZR1*ZP
Expand All @@ -5896,14 +5914,14 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
65 ZF0=ZF0+ZR0
DO 70 K=1,M
70 SP0=SP0+1.0D0/K
ZF1=PA+PB-SP0+2.0D0*EL+CDLOG(1.0D0-Z)
ZF1=PA+PB-SP0+2.0D0*EL+ZLG
DO 80 K=1,500
SP=SP+(1.0D0-A)/(K*(A+K-1.0D0))+(1.0D0-B)/(K*
& (B+K-1.0D0))
SM=0.0D0
DO 75 J=1,M
75 SM=SM+1.0D0/(J+K)
ZP=PA+PB+2.0D0*EL+SP-SM+CDLOG(1.0D0-Z)
ZP=PA+PB+2.0D0*EL+SP-SM+ZLG
ZR1=ZR1*(A+K-1.D0)*(B+K-1.D0)/(K*(M+K))*(1.D0-Z)
ZF1=ZF1+ZR1*ZP
IF (CDABS(ZF1-ZW).LT.CDABS(ZF1)*EPS) GO TO 85
Expand All @@ -5920,6 +5938,7 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
CALL GAMMA2(A+B-C,GABC)
ZC0=GC*GCAB/(GCA*GCB)
ZC1=GC*GABC/(GA*GB)*(1.0D0-Z)**(C-A-B)
IF (L8) ZC1=ZC1*CDEXP(-2.0D0*PI*(C-A-B)*(0.0D0,1.0D0))
ZHF=(0.0D0,0.0D0)
ZR0=ZC0
ZR1=ZC1
Expand All @@ -5933,24 +5952,121 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
95 ZHF=ZHF+ZC0+ZC1
ENDIF
ELSE
Z00=(1.0D0,0.0D0)
IF (C-A.LT.A.AND.C-B.LT.B) THEN
Z00=(1.0D0-Z)**(C-A-B)
A=C-A
B=C-B
ZR=0.5D0*(A+B+1.0D0)-C
MAB=INT(A-B+EPS*DSIGN(1.0D0,A-B))
IF (DABS(A-B-MAB)>EPS) THEN
CALL GAMMA2(A,GA)
CALL GAMMA2(B,GB)
CALL GAMMA2(C,GC)
CALL GAMMA2(B-A,GBA)
CALL GAMMA2(A-B,GAB)
CALL GAMMA2(C-B,GCB)
CALL GAMMA2(C-A,GCA)
ZCA=GC*GBA/(GB*GCA*(0.5D0-Z)**A)
ZCB=GC*GAB/(GA*GCB*(0.5D0-Z)**B)
ZEA(1)=(1.0D0,0.0D0)
ZEA(2)=A*ZR/(1.0D0+A-B)
ZEB(1)=(1.0D0,0.0D0)
ZEB(2)=B*ZR/(1.0D0+B-A)
ZINV=1.0D0/(Z-0.5D0)
ZHF=ZCA+ZCB+(ZCA*ZEA(2)+ZCB*ZEB(2))*ZINV
DO 900 K=2,1500
ZEA(3)=(A+K-1.0D0)/(K*(K+A-B))*
& (ZR*ZEA(2)+0.25D0*(K+A-2.0D0)*ZEA(1))
ZEB(3)=(B+K-1.0D0)/(K*(K+B-A))*
& (ZR*ZEB(2)+0.25D0*(K+B-2.0D0)*ZEB(1))
ZINV=ZINV/(Z-0.5D0)
ZTMP=(ZCA*ZEA(3)+ZCB*ZEB(3))*ZINV
ZHF=ZHF+ZTMP
IF (CDABS(ZTMP).LE.CDABS(ZHF)*EPS) GO TO 905
ZEA(1)=ZEA(2)
ZEA(2)=ZEA(3)
ZEB(1)=ZEB(2)
900 ZEB(2)=ZEB(3)
905 IF (K.GT.150) ISFER=5
ELSE
IF (A.LT.B) THEN
A=BB
B=AA
MAB=-MAB
ENDIF
CALL GAMMA2(A,GA)
CALL GAMMA2(B,GB)
CALL GAMMA2(C,GC)
CALL GAMMA2(C-B,GCB)
CALL GAMMA2(C-A,GCA)
ZLG=LOG(0.5D0-Z)
FMAB=1.0D0
DO 910 K=1,MAB
910 FMAB=FMAB*DBLE(K)
CALL PSI_SPEC(A,PA)
CALL PSI_SPEC(C-A,PCA)
PMAB=-EL
DO 915 K=1,MAB
915 PMAB=PMAB+1.0D0/DBLE(K)
ZEA(2)=(1.0D0,0.0D0)
ZEA(1)=(0.0D0,0.0D0)
ZEB(2)=(1.0D0,0.0D0)
ZEB(1)=(0.0D0,0.0D0)
ZFA(2)=(0.0D0,0.0D0)
ZFA(1)=(0.0D0,0.0D0)
ZFB(2)=(0.0D0,0.0D0)
ZFB(1)=(0.0D0,0.0D0)
ZV1=(0.0D0,0.0D0)
IF (MAB.GT.0) THEN
ZV1=ZV1+1.0D0
ZTMP=(1.0D0,0.0D0)
DO 920 K=1,MAB
ZEB(3)=(ZR*ZEB(2)+0.25D0*(K-MAB-1.0D0)*ZEB(1))/K
ZFB(3)=(0.5D0*ZEB(2)-0.25D0*ZEB(1)+
& ZR*ZFB(2)+0.25D0*(K-1.0D0-MAB)*ZFB(1))/K
ZEB(1)=ZEB(2)
ZEB(2)=ZEB(3)
ZFB(1)=ZFB(2)
ZFB(2)=ZFB(3)
IF (K.EQ.MAB) GO TO 925
ZTMP=ZTMP*(B+K-1.0D0)/((K-MAB)*(Z-0.5D0))
920 ZV1=ZV1+ZTMP*ZEB(3)
925 ZV1=GC*FMAB/(MAB*GA*GCB*(0.5D0-Z)**B)*ZV1
ENDIF
ZINV=(1.0D0,0.0D0)
UA=PCA-PMAB
UB=-EL-PA
WA=(-1)**(MAB+1)*GC/(GB*GCA*FMAB)
WB=GC/(GB*GCB)
ZHF=WA*((UA-ZLG)*ZEA(2)+ZFA(2))
& +WB*(UB*ZEB(2)+ZFB(2))
DO 930 K=1,15000
WA=(A+K-1.0D0)*WA/(MAB+K)
WB=(A+K-1.0D0)*WB/K
UA=UA+1.0D0/(A+K-1.0D0)-1.0D0/(K+MAB)
UB=UB+1.0D0/K
ZEA(3)=(ZR*ZEA(2)+0.25D0*(K+MAB-1.0D0)*ZEA(1))/K
ZEB(3)=(ZR*ZEB(2)+0.25D0*(K-1.0D0)*ZEB(1))/(K+MAB)
ZFA(3)=(0.5D0*ZEA(2)+0.25D0*ZEA(1) +
& ZR*ZFA(2)+0.25D0*(K+MAB-1.0D0)*ZFA(1))/K
ZFB(3)=(0.5D0*ZEB(2)-0.25D0*ZEB(1)+
& ZR*ZFB(2)+0.25D0*(K-1.0D0)*ZFB(1))/(K+MAB)
ZTMP=WA*((UA-ZLG)*ZEA(3)+ZFA(3))+
& WB*(UB*ZEB(3)+ZFB(3))
ZINV=ZINV/(Z-0.5D0)
ZTMP=ZTMP*ZINV
ZHF=ZHF+ZTMP
IF (CDABS(ZTMP).LE.CDABS(ZHF)*EPS) GO TO 935
ZEA(1)=ZEA(2)
ZEA(2)=ZEA(3)
ZEB(1)=ZEB(2)
ZEB(2)=ZEB(3)
ZFA(1)=ZFA(2)
ZFA(2)=ZFA(3)
ZFB(1)=ZFB(2)
930 ZFB(2)=ZFB(3)
935 ZHF=ZV1+ZHF/((0.5D0-Z)**A)
IF(K.GT.150) ISFER=5
ENDIF
ZHF=(1.0D0,0.D0)
ZR=(1.0D0,0.0D0)
DO 100 K=1,1500
ZR=ZR*(A+K-1.0D0)*(B+K-1.0D0)/(K*(C+K-1.0D0))*Z
ZHF=ZHF+ZR
IF (CDABS(ZHF-ZW).LE.CDABS(ZHF)*EPS) GO TO 105
100 ZW=ZHF
105 ZHF=Z00*ZHF
ENDIF
ELSE IF (A0.GT.1.0D0) THEN
ELSE IF (A0.GT.1.1D0) THEN
MAB=INT(A-B+EPS*DSIGN(1.0D0,A-B))
IF (DABS(A-B-MAB).LT.EPS.AND.A0.LE.1.1D0) B=B+EPS
IF (DABS(A-B-MAB).GT.EPS) THEN
CALL GAMMA2(A,GA)
CALL GAMMA2(B,GB)
Expand Down