Skip to content

Instantly share code, notes, and snippets.

@certik
Created June 17, 2012 00:22
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save certik/2942970 to your computer and use it in GitHub Desktop.
Save certik/2942970 to your computer and use it in GitHub Desktop.
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment