]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:01:57 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | #if defined(CERNLIB_SINGLE) | |
11 | SUBROUTINE C309R7(A,B,IBEG,INUM,XX,EPS) | |
12 | C | |
13 | C******************************************************************* | |
14 | C | |
15 | C RCF converts polynomial A to the corresponding continued | |
16 | C fraction, in 'normal' form with coefficients B | |
17 | C by the 'P algorithm' of Patry & Gupta | |
18 | C | |
19 | C A(z) = A1/z + A2/z**3 + A3/z**5 + ... + An/z**(2n-1) | |
20 | C | |
21 | C B(z) = B1/z+ B2/z+ B3/z+ .../(z+ Bn/z) | |
22 | C | |
23 | C data: | |
24 | C A vector A(k), k=1,INUM input | |
25 | C B vector B(k), k=IBEG,INUM output | |
26 | C IBEG order of first coef. calc. input | |
27 | C INUM order of A, even or odd input | |
28 | C XX auxiliary vector of length .ge. length of vector B | |
29 | C caller provides space for A,B,XX | |
30 | C Note that neither of the first two terms A(1) A(2) should be zero | |
31 | C & the user can start the calculation with any value of | |
32 | C IBEG provided the c.f. coefs have been already | |
33 | C calculated up to INUM = IBEG-1 | |
34 | C & the method breaks down as soon as the absolute value | |
35 | C of a c.f. coef. is less than EPS. At the time of the | |
36 | C break up INUM has been replaced by minus times the number | |
37 | C of this coefficient. | |
38 | C algorithm: J. Patry & S. Gupta, EIR-Bericht 247, November 1973 | |
39 | C Eidg. Institut fur Reaktorforschung | |
40 | C Wuerenlingen, Switzerland | |
41 | C see also: Haenggi, Roesel & Trautmann, | |
42 | C J. Comput. Phys., v. 137, (1980) 242-258 | |
43 | C note: restart procedure modified by I.J.Thompson | |
44 | C | |
45 | C******************************************************************* | |
46 | C | |
47 | IMPLICIT COMPLEX(A-H,O-Z) | |
48 | DIMENSION A(100),B(100),XX(2,100) | |
49 | LOGICAL EVEN | |
50 | REAL EPS | |
51 | ||
52 | IBN=INUM | |
53 | IF(IBEG .GT. 4) GO TO 50 | |
54 | IF(IBEG .EQ. 4) GO TO 20 | |
55 | B(1)=A(1) | |
56 | IF(IBN .GE. 2) B(2)=-A(2)/A(1) | |
57 | IF(IBN .LT. 3) RETURN | |
58 | X0=A(3)/A(2) | |
59 | XX(2,1)=B(2) | |
60 | XX(1,1)=-X0 | |
61 | XX(1,2)=0 | |
62 | B(3)=-X0-B(2) | |
63 | X0=-B(3)*A(2) | |
64 | M=3 | |
65 | MP12=2 | |
66 | EVEN=.TRUE. | |
67 | IF(IBN .LE. 3) RETURN | |
68 | 20 IF(ABS(B(3)) .LT. EPS*ABS(X0)) THEN | |
69 | INUM=-M | |
70 | RETURN | |
71 | END IF | |
72 | M=4 | |
73 | 30 X1=A(M) | |
74 | M2M1=MP12 | |
75 | MP12=M2M1+1 | |
76 | IF(EVEN) MP12=M2M1 | |
77 | DO 40 K = 2,MP12 | |
78 | 40 X1=X1+A(M-K+1)*XX(1,K-1) | |
79 | B(M)=-X1/X0 | |
80 | IF(M .GE. IBN) RETURN | |
81 | 50 IF(ABS(B(M)) .LT. EPS*ABS(X0)) THEN | |
82 | INUM=-M | |
83 | RETURN | |
84 | END IF | |
85 | DO 60 K = M2M1,2,-1 | |
86 | 60 XX(2,K)=XX(1,K)+B(M)*XX(2,K-1) | |
87 | XX(2,1)=XX(1,1)+B(M) | |
88 | DO 70 K = 1,M2M1 | |
89 | X0=XX(2,K) | |
90 | XX(2,K)=XX(1,K) | |
91 | 70 XX(1,K)=X0 | |
92 | X0=X1 | |
93 | XX(1,M2M1+1)=0 | |
94 | M=M+1 | |
95 | EVEN=.NOT.EVEN | |
96 | GO TO 30 | |
97 | END | |
98 | FUNCTION C309R8(Z,ACC) | |
99 | COMPLEX C309R8,Z | |
100 | REAL ACC | |
101 | REAL X,Y,AX,AY,A | |
102 | ||
103 | X=REAL(Z) | |
104 | Y=AIMAG(Z) | |
105 | AX=ABS(X) | |
106 | AY=ABS(Y) | |
107 | A=5*ACC*(AX+AY) | |
108 | IF(AX .LT. A) X=0 | |
109 | IF(AY .LT. A) Y=0 | |
110 | C309R8=CMPLX(X,Y) | |
111 | RETURN | |
112 | END | |
113 | #endif |