This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / rmullz64.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:01:52  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10 #if defined(CERNLIB_DOUBLE)
11       SUBROUTINE DMULLZ(A,N,MAXITR,Z0)
12 #endif
13 #if !defined(CERNLIB_DOUBLE)
14       SUBROUTINE RMULLZ(A,N,MAXITR,Z0)
15 #endif
16 #include "gen/imp64.inc"
17 #include "gen/defc64.inc"
18      +  Z0,Z,DX,X,X3,Y1,Y2,Y,TE(7)
19       LOGICAL LSW
20       DIMENSION A(0:*),Z0(*)
21       CHARACTER*(*) NAME
22 #if defined(CERNLIB_DOUBLE)
23       PARAMETER( NAME='DMULLZ')
24 #endif
25 #if !defined(CERNLIB_DOUBLE)
26       PARAMETER( NAME='RMULLZ')
27 #endif
28       PARAMETER (Z1 = 1, HF = Z1/2)
29       PARAMETER (ETA1 = 1D-14, ETA2 = 6D-8, BIG = 1D20)
30       PARAMETER (C1 = 0.9D0, C2 = 0.002D0, C3 = 0.1D0)
31
32 #include "gen/gcmpfun.inc"
33       SUMABS(Z)=ABS(GREAL(Z))+ABS(GIMAG(Z))
34
35       IF(N .LE. 0) RETURN
36       IF(A(0) .EQ. 0) THEN
37        CALL MTLPRT(NAME,'C202.1','A(0) = 0')
38        RETURN
39       ENDIF
40
41       N1=N
42     2 IF(N1 .EQ. 1) THEN
43        Z0(1)=-A(1)/A(0)
44        RETURN
45       ENDIF
46       IF(A(N1) .EQ. 0) THEN
47        Z0(N1)=0
48        N1=N1-1
49        GO TO 2
50       ENDIF
51
52       SCALE=ABS(A(N1)/A(0))**(Z1/N1)
53       B=A(0)
54       DO 6 I = 1,N1
55       B=B*SCALE
56     6 Z0(I)=A(I)/B
57       IF(N1 .EQ. 2) GO TO 1
58
59    10 LSW=.TRUE.
60       Y1=Z0(1)+1
61       Y2=Z0(1)-1
62       DO 11 I = 2,N1
63       Y1=Z0(I)+Y1
64    11 Y2=Z0(I)-Y2
65       Y=Z0(N1)
66       X=0
67       DX=1
68
69       TE(1)=-2
70    12 TE(2)=Y2/Y
71       TE(3)=(Y1-Y2)/(Y*TE(1))
72       DO 17 ITER = 1,MAXITR
73       TE(4)=TE(2)-1
74       TE(5)=(TE(4)-TE(3))/(TE(1)+1)
75       TE(6)=HF*(TE(5)+TE(4))
76       TE(7)=SQRT(TE(6)**2+TE(5))
77       TE(1)=TE(6)+TE(7)
78       TE(7)=TE(6)-TE(7)
79       IF(ABS(TE(1)) .LE. ABS(TE(7))) THEN
80        IF(TE(7) .EQ. 0) TE(7)=C1
81        TE(1)=TE(7)
82       ENDIF
83       DX=DX/TE(1)
84       X=DX+X
85       EPSI=SUMABS(X)*ETA1
86       IF(SUMABS(DX) .LT. EPSI .AND. SUMABS(Y) .LT. C2) GO TO 18
87       Y2=Y
88       Y=X+Z0(1)
89       DO 21 I = 2,N1
90    21 Y=Y*X+Z0(I)
91       IF(Y .EQ. 0) GO TO 18
92
93    15 IF(SUMABS(Y) .GE. 100*SUMABS(Y2) .AND. SUMABS(DX) .GE. EPSI) THEN
94        TE(1)=2*TE(1)
95        DX=HF*DX
96        X=X-DX
97        Y=X+Z0(1)
98        DO 22 I = 2,N1
99    22  Y=Y*X+Z0(I)
100        IF(Y .EQ. 0) GO TO 18
101        GOTO 15
102       ENDIF
103
104       TE(2)=Y2/Y
105       TE(3)=(TE(2)/TE(1))*TE(4)
106    17 CONTINUE
107
108       CN=ABS(Z0(N1))
109       IF(ABS(CN-1) .LT. C3) THEN
110        DO 40 I = 1,N1
111    40  Z0(I)=BIG
112        CALL MTLPRT(NAME,'C202.2','TOO MANY ITERATIONS')
113        RETURN
114       ENDIF
115       S=CN**(-Z1/N1)
116       SCALE=SCALE/S
117       B=1
118       DO 30 I = 1,N1
119       B=B*S
120    30 Z0(I)=Z0(I)*B
121       GO TO 10
122
123    20 IF(ABS(GIMAG(X)) .LT. ABS(GREAL(X))*ETA2) GO TO 10
124       LSW=.FALSE.
125       X3=GCONJG(X)
126       DX=GCONJG(DX)
127       TE(1)=GCONJG(TE(1))
128       X=X3-DX
129       Y=X+Z0(1)
130       DO 51 I = 2,N1
131    51 Y=Y*X+Z0(I)
132       IF(Y .EQ. 0) GO TO 18
133       Y2=Y
134       X=X-DX*TE(1)
135       Y=X+Z0(1)
136       DO 52 I = 2,N1
137    52 Y=Y*X+Z0(I)
138       IF(Y .EQ. 0) GO TO 18
139       Y1=Y
140       X=X3
141       Y=X+Z0(1)
142       DO 53 I = 2,N1
143    53 Y=Y*X+Z0(I)
144       IF(Y .NE. 0) GO TO 12
145
146    18 Z0(N1)=X*SCALE
147       N1=N1-1
148       Z0(1)=X+Z0(1)
149       DO 19 I = 2,N1
150    19 Z0(I)=Z0(I-1)*X+Z0(I)
151       IF(N1 .GT. 2) THEN
152        IF(LSW) GO TO 20
153        GO TO 10
154       ENDIF
155
156     1 Z0(2)=(SQRT((HF*Z0(1))**2-Z0(2))-HF*Z0(1))*SCALE
157       Z0(1)=-Z0(1)*SCALE-Z0(2)
158       RETURN
159       END