[Fortran] 纯文本查看 复制代码
SUBROUTINE DEXPI(X,VV,ICON) 00010580
IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00010590
DIMENSION EXWK(10),DTA(-40:40) 00010600
PARAMETER(MMM=-38,LLL=-1) 00010610
DATA (DTA(I),I=MMM,LLL)/-0.426609511901546924D+1, 00010620
& -0.424077731103007936D+1,-0.421480182462791870D+1, 00010630
& -0.418813357754575738D+1,-0.416073460335764294D+1, 00010640
& -0.413256372639094662D+1,-0.410357618951769436D+1, 00010650
& -0.407372322636801330D+1,-0.404295156770125985D+1, 00010660
& -0.401120286938668024D+1,-0.397841304656369131D+1, 00010670
& -0.394451149488801543D+1,-0.390942017507676071D+1, 00010680
& -0.387305253090592919D+1,-0.383531220292320443D+1, 00010690
& -0.379609148977026855D+1,-0.375526949525099083D+1, 00010700
& -0.371270988083496528D+1,-0.366825811827199814D+1, 00010710
& -0.362173810265948791D+1,-0.357294793855387915D+1, 00010720
& -0.352165464434875627D+1,-0.346758742360130662D+1, 00010730
& -0.341042901126418097D+1,-0.334980439378180138D+1, 00010740
& -0.328526588519002188D+1,-0.321627305017310262D+1, 00010750
& -0.314216518454969151D+1,-0.306212278979771567D+1, 00010760
& -0.297511233800020768D+1,-0.287980491486450826D+1, 00010770
& -0.277445268959193027D+1,-0.268669487220421274D+1, 00010780
& -0.252324129568845652D+1,-0.236933521658175054D+1, 00010790
& -0.218780187292690857D+1,-0.196728937843127240D+1, 00010800
& -0.168887633466383960D+1/ 00010810
DATA (DTA(I),I=1,38)/0.825800461705577393D+1, 00010820
& 0.176673644440347963D+2,0.379986217784675437D+2, 00010830
& 0.836207870083096158D+2,0.188981617521544546D+3, 00010840
& 0.437723242328256892D+3,0.103510385047485182D+4, 00010850
& 0.248934917548398214D+4,0.606843126316091148D+4, 00010860
& 0.149564705440828390D+5,0.371945463256666719D+5, 00010870
& 0.931892973609708525D+5,0.234952567224902294D+6, 00010880
& 0.595557648866449846D+6,0.151663448361350786D+7, 00010890
& 0.387790086301002053D+7,0.995090372939220018D+7, 00010900
& 0.256156490911086490D+8,0.661271827337468153D+8, 00010910
& 0.171144667632105541D+9,0.443966366117561315D+9, 00010920
& 0.115411538809391339D+10,0.300595090272945703D+10, 00010930
& 0.784294098806287371D+10,0.204964971160077586D+11, 00010940
& 0.536451185884052706D+11,0.140599195754462385D+12, 00010950
& 0.368973209403295761D+12,0.969455575964382700D+12, 00010960
& 0.255004356635374381D+13,0.671464018407242340D+13, 00010970
& 0.176980372441121638D+14,0.466905501446574597D+14, 00010980
& 0.123285207991205599D+15,0.325798899867222188D+15, 00010990
& 0.861638819996574379D+15,0.228044620030189819D+16, 00011000
& 0.603971826361123688D+16/ 00011010
R=0.57721 56649 01532 86061D0 00011020
EPS=1.0D-17 00011030
ICON=0 00011040
X1=ABS(X) 00011050
IF(ABS(X1).LT.EPS) THEN 00011060
ICON=1 00011070
RETURN 00011080
ENDIF 00011090
IF(X1.GE.40.0D0) GO TO 80 00011100
IF(X1.GT.3.0D0) GO TO 30 00011110
C TAYLOR EXPANSION 00011120
VV=R+DLOG(X1) 00011130
DS=1.0D0 00011140
WK=1.0D0 00011150
N=100 00011160
DO 10 I=1,N 00011170
WK=X*WK/DS 00011180
WK1=WK/DS 00011190
VV=VV+WK1 00011200
IF(ABS(WK1).LT.EPS) RETURN 00011210
DS=DS+1.0D0 00011220
10 CONTINUE 00011230
ICON=2 00011240
RETURN 00011250
C 00011260
C ROMBERG METHOD 00011270
30 VV=R+DLOG(X1) 00011280
IX=X1+0.5D0 00011290
X0=FLOAT(IX) 00011300
H=X1-X0 00011310
N=1 00011320
SOLD=0.0D0 00011330
WK=((DEXP(SIGN(X0,X))-1.0D0)/X0+(DEXP(SIGN(X1,X))-1.0D0) 00011340
& /X1)*H*0.5D0 00011350
ITR=12 00011360
DO 60 M=1,ITR 00011370
EXWK(M)=WK 00011380
EXD=4.0D0 00011390
J=M-1 00011400
DO 40 K=2,M 00011410
EXWK(J)=EXWK(J+1)+(EXWK(J+1)-EXWK(J))/(EXD-1.0D0) 00011420
EXD=EXD*4.0D0 00011430
J=J-1 00011440
40 CONTINUE 00011450
S=EXWK(1) 00011460
IF(ABS(S).LT.EPS) GO TO 65 00011470
IF(ABS(S-SOLD).LT.EPS) GO TO 65 00011480
SOLD=S 00011490
CC=0.0D0 00011500
X2=X0-H*0.5D0 00011510
DO 50 MM=1,N 00011520
X2=X2+H 00011530
CC=CC+(DEXP(SIGN(X2,X))-1.0D0)/X2 00011540
50 CONTINUE 00011550
CC=CC*H 00011560
WK=(WK+CC)*0.5D0 00011570
H=H*0.5D0 00011580
N=N*2 00011590
60 CONTINUE 00011600
ICON=2 00011610
65 ONDO=ABS(X)/X*(IX-2)
INDO=INT(ONDO)
VV=VV+S+DTA(INDO) 00011620
RETURN 00011630