The source code was obtained from: http://www.netlib.org/toms/579
The theoretical article documenting the code is: http://dl.acm.org/citation.cfm?id=355979
The source code was obtained from: http://www.netlib.org/toms/579
The theoretical article documenting the code is: http://dl.acm.org/citation.cfm?id=355979
C CPSC - COMPLEX POWER SERIES COEFFICIENTS | |
C | |
C BY BENGT FORNBERG | |
C | |
C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, DECEMBER 1981 | |
C | |
C PROGRAM TEST (OUTPUT) | |
C | |
C THIS PROGRAM WITH THE FUNCTION F BELOW USES CPSC TO EVALUATE | |
C THE LEADING DERIVATIVES OF F(Z) = CEXP(Z)/(CSIN(Z)**3+CCOS(Z)**3) | |
C AT Z = 0. | |
EXTERNAL F | |
DIMENSION A(51),IRANGE(4),ER(51) | |
COMPLEX A | |
DATA IRANGE/6,12,25,51/ | |
DO 10 I=1,4 | |
IR = IRANGE(I) | |
R = 1. | |
CALL CPSC(F,(0.,0.),IR,1,R,A,ER) | |
WRITE(6,20) I | |
WRITE(6,30) (A(J),J=1,IR) | |
10 WRITE(6,40) (ER(J),J=1,IR) | |
20 FORMAT(/20H DERIVATIVES , RANGE,I3) | |
30 FORMAT(4(1X,E18.10,E9.1)) | |
40 FORMAT(/17H ESTIMATED ERRORS/(1X,16E8.1)) | |
STOP | |
END | |
COMPLEX FUNCTION F(Z) | |
C | |
C TEST FUNCTION FOR USE WITH CPSC. | |
C | |
COMPLEX Z | |
F = CEXP(Z)/(CSIN(Z)**3+CCOS(Z)**3) | |
RETURN | |
END | |
SUBROUTINE CPSC(F,Z,N,IC,R,RS,ER) | |
C | |
C EVALUATION OF COMPLEX POWER SERIES COEFFICIENTS OR DERIVATIVES. | |
C | |
C *** INPUT PARAMETERS *** | |
C F COMPLEX FUNCTION, OF WHICH THE COEFFICIENTS OR DERIVATIVES | |
C ARE SOUGHT. THIS FUNCTION MUST BE DECLARED EXTERNAL IN THE | |
C CALLING PROGRAM. | |
C Z COMPLEX POINT AROUND WHICH F SHALL BE EXPANDED OR AT WHICH | |
C DERIVATIVES SHALL BE EVALUATED. | |
C N INTEGER, NUMBER OF COEFFICIENTS OR DERIVATIVES WANTED. | |
C N MUST BE GE 1 AND LE 51. | |
C IC SELECTS BETWEEN POWER SERIES COEFFICIENTS AND DREIVATIVES. | |
C IC .EQ. 0 ROUTINE RETURNS POWER SERIES COEFFICIENTS IN RS | |
C IC .NE. 0 ROUTINE RETURNS DERIVATIVES IN RS . | |
C *** INPUT AND OUTPUT PARAMETER *** | |
C R INITIAL RADIUS USED IN SEARCH FOR OPTIMAL RADIUS. THE RESULTING | |
C RADIUS IS LEFT IN R. THE PROVIDED GUESS MAY BE IN ERROR WITH AT | |
C MOST A FACTOR OF 3.E4 . | |
C *** OUTPUT PARAMETERS *** | |
C RS COMPLEX ARRAY RS(N) CONTAINING THE N FIRST | |
C COEFFICIENTS (CORRESPONDING TO THE POWERS 0 TO N-1) OR DERIVA- | |
C TIVES (ORDERS 0 TO N-1) . | |
C ER REAL ARRAY ER(N) CONTAINING ERROR ESTIMATES FOR THE | |
C NUMBERS IN RS(N). | |
C | |
DIMENSION IP(32),A(64),RS(N),ER(N),RT(51,3),FV(6), | |
* IW(7),SC(4),RV(3),C(4),FC(3) | |
COMPLEX F,A,V,RS,RT,FV,U,W,T,Z,RV,RQ,S,XK,MULT,CO | |
C | |
C LIST OF THE VARIABLES INITIALIZED IN THE DATA STATEMENT BELOW. | |
C EPS MACHINE ACCURACY. THIS CONSTANT HAS TO BE SUPPLIED BY | |
C THE USER. IN THIS LIST, IT IS GIVEN AS 1.E-14 CORRESPONDING | |
C TO THE 48 BIT FLOATING POINT MANTISSAS ON CDC CYBER MACHINES. | |
C IND INTEGER FLAG. | |
C L2 INTEGER FLAG. | |
C IW 2**( 0 , 1 , 2 , 3 , 4 , 5 , 6 ) . | |
C IP PERMUTATION CONSTANTS FOR THE FFT. | |
C RV CONSTANTS FOR THE LAURENT SERIES TEST . | |
C | |
DATA EPS/1.E-14/,IND/0/,L2/1/,IW/1,2,4,8,16,32,64/, | |
* IP/64,32,48,16,56,24,40,8,60,28,44,12,52,20,36,4,62,30,46,14, | |
* 54,22,38,6,58,26,42,10,50,18,34,2/, | |
* RV/(-.4,.3),(.7,.2),(.02,-.06)/ | |
C | |
C STATEMENT FUNCTION FOR MULTIPLICATION OF A COMPLEX NUMBER | |
C BY A REAL NUMBER. | |
C | |
MULT(RE,CO) = CMPLX(RE*REAL(CO),RE*AIMAG(CO)) | |
C | |
C EVALUATE SOME CONSTANTS THE FIRST TIME THE CODE IS EXECUTED. | |
C | |
IF(IND.EQ.1) GOTO 20 | |
IND = 1 | |
SC(1) = .125 | |
C(1) = EPS**(1./28.) | |
EP6 = C(1)**6 | |
PI = 4.*ATAN(1.) | |
FV(1) = (-1.,0.) | |
FV(2) = (0.,-1.) | |
R1 = SQRT(.5) | |
RA = 1./R1 | |
FV(3) = CMPLX(R1,-R1) | |
DO 10 I=2,4 | |
SC(I) = .5*SC(I-1) | |
C(I) = SQRT(C(I-1)) | |
ANG = PI*SC(I-1) | |
10 FV(I+2) = CMPLX(COS(ANG),-SIN(ANG)) | |
20 CONTINUE | |
C | |
C START EXECUTION. | |
C | |
IF(N.GT.51.OR.N.LT.1) GOTO 260 | |
LF = 0 | |
NP = 0 | |
M = 0 | |
NR = -1 | |
C | |
C FIND IF A FFT OVER 8, 16, 32, OR 64 POINTS SHOULD BE USED. | |
C | |
KL = 1 | |
IF(N.GT.6) KL=2 | |
IF(N.GT.12) KL=3 | |
IF(N.GT.25) KL=4 | |
KM = KL+2 | |
KN = 7-KM | |
IX = IW(KM+1) | |
IS = IW(KN) | |
30 V = CMPLX(R,0.) | |
C | |
C FUNCTION VALUES OF F ARE STORED READY PERMUTATED FOR THE FFT. | |
C | |
DO 40 I=IS,32,IS | |
IQ = IP(I) | |
V = V*FV(KM) | |
A(IQ) = F(Z+V) | |
40 A(IQ-1) = F(Z-V) | |
LN = 2 | |
JN = 1 | |
C | |
C THE LOOP DO 70 ... CONSTITUTES THE FFT. | |
C | |
DO 70 L=1,KM | |
U = (1.,0.) | |
W = FV(L) | |
DO 60 J=1,JN | |
DO 50 I=J,IX,LN | |
IT = I+JN | |
T = A(IT)*U | |
A(IT) = A(I)-T | |
50 A(I) = A(I)+T | |
60 U = U*W | |
LN = LN+LN | |
70 JN = JN+JN | |
CX = 0. | |
B = 1. | |
C | |
C TEST ON HOW FAST THE COEFFICIENTS OBTAINED DECREASE. | |
C | |
DO 80 I=1,IX | |
CT = CABS(A(I))/B | |
IF(CT.LT.CX) GOTO 80 | |
CX = CT | |
INR = I | |
80 B = B*C(KL) | |
IF(M.LE.1) GOTO 100 | |
C | |
C ESTIMATE OF THE ROUNDING ERROR LEVEL FOR THE LAST RADIUS. | |
C | |
ER(1) = CX*EPS | |
DO 90 I=2,N | |
90 ER(I) = ER(I-1)/R | |
100 SF = SC(KL) | |
DO 110 I=1,IX | |
A(I) = MULT(SF,A(I)) | |
110 SF = SF/R | |
L1 = L2 | |
L2 = 1 | |
IF(INR.GT.IW(KM)) GOTO 150 | |
IF(LF.EQ.1) GOTO 140 | |
C | |
C TEST IF THE SERIES IS A TAYLOR OR A LAURENT SERIES. | |
C | |
SR = 0. | |
SP = 0. | |
DO 130 J=1,3 | |
RQ = MULT(R,RV(J)) | |
S = A(IX) | |
DO 120 I=2,IX | |
IA = IX+1-I | |
120 S = S*RQ+A(IA) | |
CP = CABS(S) | |
IF(CP.GT.SP) SP=CP | |
CM = CABS(S-F(Z+RQ)) | |
130 IF(CM.GT.SR) SR=CM | |
IF(SR.GT.1.E-3*SP) GOTO 150 | |
LF = 1 | |
140 L2 = -1 | |
C | |
C DETERMINATION OF THE NEXT RADIUS TO BE USED. | |
C | |
150 IF(NR.GE.0) GOTO 160 | |
FACT = 2. | |
IF(L2.EQ.1) FACT=.5 | |
L1 = L2 | |
NR = 0 | |
160 IF(L1.NE.L2) GOTO 180 | |
IF(NR.GT.0) GOTO 170 | |
NP = NP+1 | |
IF(NP-15) 190,190,260 | |
170 FACT = 1./FACT | |
180 FACT = 1./SQRT(FACT) | |
NR = NR+1 | |
190 R = R*FACT | |
M = NR-KL-1 | |
IF(M.LE.0) GOTO 30 | |
C | |
C THE RESULTS FOR THE LAST THREE RADII ARE STORED. | |
C | |
DO 200 I=1,N | |
200 RT(I,M) = A(I) | |
IF(M.EQ.1) GOTO 220 | |
C | |
C EXTRAPOLATION. | |
C | |
DO 210 I=1,N | |
XK = RT(I,M-1)-RT(I,M) | |
210 RT(I,M-1) = RT(I,M)-MULT(FC(M-1),XK) | |
IF(M.EQ.3) GOTO 230 | |
C | |
C CALCULATION OF THE EXTRAPOLATION CONSTANTS. | |
C | |
220 FC(M) = 1.5+SIGN(.5,FACT-1.) | |
IF(M.EQ.2) FC(M)=FC(M)+RA | |
IF(FACT.GT.1.) FC(M)=-FC(M) | |
GOTO 30 | |
230 FC(3) = FC(1)*FC(2)/(FC(1)+FC(2)+1.) | |
C | |
C FINAL EXTRAPOLATION AND ERROR ESTIMATE. | |
C | |
DO 240 I=1,N | |
XK = RT(I,1)-RT(I,2) | |
ER(I) = ER(I)+EP6*CABS(XK) | |
240 RS(I) = RT(I,2)-MULT(FC(3),XK) | |
C | |
C MULTIPLY POWER SERIES COEFFICIENTS AND ERROR ESTIMATE BY FACTORIALS | |
C IF DERIVATIVES WANTED. | |
C | |
IF(IC.EQ.0) RETURN | |
FAC = 0. | |
FACT = 1. | |
DO 250 I=1,N | |
RS(I) = MULT(FACT,RS(I)) | |
ER(I) = FACT*ER(I) | |
FAC = FAC+1. | |
250 FACT = FACT*FAC | |
RETURN | |
C | |
C ERROR RETURN. | |
C | |
260 DO 270 I=1,N | |
RS(I) = (0.,0.) | |
270 ER(I) = 1.E10 | |
RETURN | |
END |
$ gfortran -Wall -Wextra -Wimplicit-interface -fPIC -g -fcheck=all -fbacktrace -fdefault-real-8 a.f | |
a.f:131.16: | |
A(IQ) = F(Z+V) | |
1 | |
Warning: Procedure 'f' called with an implicit interface at (1) | |
a.f:132.18: | |
40 A(IQ-1) = F(Z-V) | |
1 | |
Warning: Procedure 'f' called with an implicit interface at (1) | |
a.f:189.21: | |
CM = CABS(S-F(Z+RQ)) | |
1 | |
Warning: Procedure 'f' called with an implicit interface at (1) | |
a.f:19.72: | |
CALL CPSC(F,(0.,0.),IR,1,R,A,ER) | |
1 | |
Warning: Procedure 'cpsc' called with an implicit interface at (1) | |
a.f:36:0: warning: unused parameter ‘f’ [-Wunused-parameter] |
$ ./a.out | |
DERIVATIVES , RANGE 1 | |
0.1000000000E+01 0.6E-16 0.9999999999E+00 0.1E-15 0.4000000000E+01 -0.3E-14 0.3999999999E+01 -0.3E-13 | |
0.2800000000E+02 -0.2E-12 -0.1640000000E+03 -0.3E-11 | |
ESTIMATED ERRORS | |
0.4E-10 0.5E-10 0.1E-09 0.5E-09 0.3E-08 0.2E-07 | |
DERIVATIVES , RANGE 2 | |
0.9999819306E+00 -0.5E-16 0.1000025657E+01 0.2E-15 0.3999936778E+01 0.3E-15 0.4000235832E+01 0.9E-15 | |
0.2799874093E+02 0.1E-13 -0.1639923104E+03 -0.5E-12 0.6393998808E+02 -0.6E-11 -0.1337546537E+05 -0.1E-09 | |
0.4724261903E+05 -0.9E-09 -0.8581614514E+06 -0.1E-08 0.1382903533E+08 0.5E-06 -0.1126947760E+09 0.9E-05 | |
ESTIMATED ERRORS | |
0.2E-12 0.5E-12 0.2E-11 0.1E-10 0.1E-09 0.1E-08 0.1E-07 0.2E-06 0.3E-05 0.5E-04 0.1E-02 0.3E-01 | |
DERIVATIVES , RANGE 3 | |
-Infinity-Infinity Infinity-Infinity -Infinity-Infinity Infinity Infinity | |
-Infinity Infinity Infinity Infinity -Infinity Infinity Infinity Infinity | |
-Infinity-Infinity Infinity-Infinity -Infinity-Infinity Infinity-Infinity | |
-Infinity Infinity Infinity Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity-Infinity Infinity Infinity -Infinity Infinity Infinity Infinity | |
-Infinity Infinity Infinity-Infinity -Infinity-Infinity Infinity-Infinity | |
-Infinity Infinity | |
ESTIMATED ERRORS | |
0.4E-12 0.7E-12 0.2E-11 0.1E-10 0.7E-10 0.6E-09 0.6E-08 0.7E-07 0.9E-06 0.1E-04 0.2E-03 0.4E-02 0.7E-01 0.2E+01 0.3E+02 0.8E+03 | |
0.2E+05 0.6E+06 0.2E+08 0.5E+09 0.2E+11 0.6E+12 0.2E+14 0.8E+15 0.3E+17 | |
DERIVATIVES , RANGE 4 | |
-Infinity Infinity Infinity Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity Infinity Infinity-Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity Infinity Infinity-Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity Infinity Infinity-Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity Infinity Infinity-Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity Infinity Infinity Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity Infinity Infinity-Infinity -Infinity Infinity Infinity Infinity | |
-Infinity Infinity Infinity Infinity -Infinity Infinity Infinity Infinity | |
-Infinity Infinity Infinity Infinity -Infinity-Infinity Infinity-Infinity | |
-Infinity-Infinity Infinity-Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity-Infinity Infinity-Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity Infinity Infinity Infinity -Infinity Infinity Infinity-Infinity | |
-Infinity Infinity Infinity-Infinity -Infinity Infinity | |
ESTIMATED ERRORS | |
0.8E-12 0.1E-11 0.3E-11 0.1E-10 0.8E-10 0.6E-09 0.5E-08 0.5E-07 0.6E-06 0.8E-05 0.1E-03 0.2E-02 0.3E-01 0.6E+00 0.1E+02 0.2E+03 | |
0.6E+04 0.1E+06 0.4E+07 0.1E+09 0.3E+10 0.9E+11 0.3E+13 0.9E+14 0.3E+16 0.1E+18 0.4E+19 0.2E+21 0.7E+22 0.3E+24 0.1E+26 0.5E+27 | |
0.2E+29 0.1E+31 0.6E+32 0.3E+34 0.1E+36 0.8E+37 0.4E+39 0.2E+41 0.1E+43 0.8E+44 0.5E+46 0.3E+48 0.2E+50 0.1E+52 0.8E+53 0.5E+55 | |
0.4E+57 0.3E+59 0.2E+61 |