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
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
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