This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / b / rsrtnt64.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:01:51  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10 #if defined(CERNLIB_DOUBLE)
11       SUBROUTINE DSRTNT(K,N,A,B,C,U1,V1,RES,LRL)
12 #endif
13 #if !defined(CERNLIB_DOUBLE)
14       SUBROUTINE RSRTNT(K,N,A,B,C,U1,V1,RES,LRL)
15 #endif
16 #include "gen/imp64.inc"
17       LOGICAL LRL,LLL
18       CHARACTER NAME*(*)
19       CHARACTER*80 ERRTXT
20 #if defined(CERNLIB_DOUBLE)
21       PARAMETER (NAME = 'DSRTNT')
22 #endif
23 #if !defined(CERNLIB_DOUBLE)
24       PARAMETER (NAME = 'RSRTNT')
25 #endif
26       PARAMETER (R1 = 1, HF = R1/2)
27
28       DIMENSION BK(0:4,0:4),SGN(0:4)
29
30       DATA BK(0,0) /1/, (BK(1,J),J=0,1) /1,1/, (BK(2,J),J=0,2) /1,2,1/
31       DATA (BK(3,J),J=0,3) /1,3,3,1/, (BK(4,J),J=0,4) /1,4,6,4,1/
32       DATA (SGN(J),J=0,4) /1,-1,1,-1,1/
33
34       P(X)=A+B*X+C*X**2
35       RT(X)=SQRT(P(X))
36
37       H=0
38       LLL=ABS(A)+ABS(B)+ABS(C) .GT. 0
39       IF(.NOT.LLL) THEN
40        CALL MTLPRT(NAME,'B300.1','A = B = C = 0')
41        GO TO 9
42       ENDIF
43       LLL=ABS(K) .LE. 3 .AND. (N .EQ. 1 .OR. N .EQ. 3)
44       IF(.NOT.LLL) THEN
45        WRITE(ERRTXT,102) K,N
46        CALL MTLPRT(NAME,'B300.2',ERRTXT)
47        GO TO 9
48       ENDIF
49
50       U=MIN(U1,V1)
51       V=MAX(U1,V1)
52       LLL=U .EQ. V
53       IF(LLL) GO TO 9
54
55       LLL=K .GE. 0 .OR. K .EQ. -1 .AND. U*V .NE. 0
56      1             .OR. K .LE. -2 .AND. U*V .GT. 0
57       IF(.NOT.LLL) GO TO 9
58
59       DELTA=4*A*C-B**2
60       RTD=SQRT(ABS(DELTA))
61       IF(C .EQ. 0) THEN
62        IF(B .LT. 0) THEN
63         X0=-A/B
64         LLL=U .LE. X0 .AND. V .LE. X0
65        ELSEIF(B .EQ. 0) THEN
66         LLL=A .GT. 0
67        ELSE
68         X0=-A/B
69         LLL=U .GE. X0 .AND. V .GE. X0
70        ENDIF
71       ELSE
72        IF(DELTA .GT. 0) THEN
73         LLL=C .GT. 0
74        ELSEIF(DELTA .EQ. 0) THEN
75         IF(C .LT. 0) THEN
76          LLL=.FALSE.
77         ELSE
78          X0=-B/(2*C)
79          LLL=U .LT. X0 .AND. V .LT. X0 .OR. U .GT. X0 .AND. V .GT. X0
80         ENDIF
81        ELSE
82         A1=(-B+RTD)/(2*C)
83         A2=(-B-RTD)/(2*C)
84         W1=MIN(A1,A2)
85         W2=MAX(A1,A2)
86         IF(C .LT. 0) THEN
87          LLL=W1 .LE. U .AND. U .LE. W2 .AND. W1 .LE. V .AND. V .LE. W2
88         ELSE
89          LLL=U .LE. W1 .AND. V .LE. W1 .OR. U .GE. W2 .AND. V .GE. W2
90         ENDIF
91        ENDIF
92       ENDIF
93       IF(.NOT.LLL) GO TO 9
94
95       IF(K .GE. 0) THEN
96        IF(C .EQ. 0) THEN
97         IF(B .EQ. 0) THEN
98          K1=K+1
99          H=(V**K1-U**K1)/(K1*SQRT(A)**N)
100         ELSE
101          N1=2-N
102          XV=A+B*V
103          XU=A+B*U
104          IF(A .EQ. 0) THEN
105           K1=2*K+N1
106           H=2*(SQRT(XV)**K1-SQRT(XU)**K1)/(K1*B**(K+1))
107          ELSE
108           HV=SQRT(XV)**N1
109           HU=SQRT(XU)**N1
110           H1=-XV/A
111           H2=-XU/A
112           S=N1*(HV-HU)
113           DO 1 J = 1,K
114     1     S=S+BK(K,J)*(H1**J*HV-H2**J*HU)/(2*J+N1)
115           H=2*(-A/B)**K*S/B
116          ENDIF
117         ENDIF
118        ELSE
119         IF(N .EQ. 1) THEN
120          ASSIGN 11 TO JMP1
121          GO TO 10
122    11    IF(K .EQ. 0) THEN
123           H=H
124          ELSEIF(K .EQ. 1) THEN
125           H=(RT(V)-RT(U)-HF*B*H)/C
126          ELSEIF(K .EQ. 2) THEN
127           H1=4*C
128           H2=6*B
129           H=((H1*V-H2)*RT(V)-(H1*U-H2)*RT(U)+(2*B**2-DELTA)*H)/(8*C**2)
130          ELSEIF(K .EQ. 3) THEN
131           H1=8*C**2
132           H2=10*B*C
133           G1=15*B**2
134           G2=A*C
135           H3=G1-16*G2
136           H=(((H1*V-H2)*V+H3)*RT(V)-((H1*U-H2)*U+H3)*RT(U)-
137      1       (HF*G1-18*G2)*B*H)/(24*C**3)
138          ENDIF
139         ELSE
140          IF(DELTA .EQ. 0) THEN
141           IF(B .EQ. 0) THEN
142            IF(K .EQ. 2) THEN
143             H=LOG(ABS(V/U))/SQRT(C)**3
144            ELSE
145             K1=K-2
146             H=(V**K1-U**K1)/(K1*SQRT(C)**3)
147            ENDIF
148            H=SIGN(R1,U)*H
149           ELSE
150            X0=B/(2*C)
151            XV=V+X0
152            XU=U+X0
153            HV=1/XV**2
154            HU=1/XU**2
155            H1=-XV/X0
156            H2=-XU/X0
157            S=HF*(HU-HV)
158            DO 2 J = 1,K
159            IF(J .NE. 2) THEN
160             S=S+BK(K,J)*(H1**J*HV-H2**J*HU)/(J-2)
161            ELSE
162             S=S+BK(K,2)*LOG(ABS(XV/XU))/X0**2
163            ENDIF
164     2      CONTINUE
165            H=(-X0)**K*S/SQRT(C)**3
166            IF(U .LT. -X0) H=-H
167           ENDIF
168          ELSE
169           IF(K .EQ. 0) THEN
170            H1=2*C
171            H=2*((H1*V+B)/RT(V)-(H1*U+B)/RT(U))/DELTA
172           ELSEIF(K .EQ. 1) THEN
173            H1=2*A
174            H=2*((H1+B*U)/RT(U)-(H1+B*V)/RT(V))/DELTA
175            LB1=11
176           ELSEIF(K .EQ. 2) THEN
177            ASSIGN 12 TO JMP1
178            GO TO 10
179    12      H1=DELTA-B**2
180            H2=2*A*B
181            H=(((H1*U-H2)/RT(U)-(H1*V-H2)/RT(V))/DELTA+H)/C
182           ELSEIF(K .EQ. 3) THEN
183            ASSIGN 13 TO JMP1
184            GO TO 10
185    13      H1=C*DELTA
186            G1=A*C
187            G2=3*B**2
188            H2=B*(10*G1-G2)
189            H3=A*(8*G1-G2)
190            H=(2*(((H1*V+H2)*V+H3)/RT(V)-((H1*U+H2)*U+H3)/RT(U))/
191      1       DELTA-3*B*H)/(2*C**2)
192           ENDIF
193          ENDIF
194         ENDIF
195        ENDIF
196       ELSE
197        IF(A .EQ. 0) THEN
198         IF(B .EQ. 0) THEN
199          K1=K-N+1
200          H=SIGN(R1,U)*(V**K1-U**K1)/(K1*SQRT(C)**N)
201         ELSE
202          IF(C .EQ. 0) THEN
203           K1=2*K-N+2
204           H=2*(SQRT(B*V)**K1-SQRT(B*U)**K1)/(K1*B**(K+1))
205          ELSE
206           XV=SQRT(C+B/V)
207           XU=SQRT(C+B/U)
208           N1=N-2
209           K1=-K+N1
210           S=0
211           DO 4 J = 0,K1
212           KJ=2*J-N1
213     4     S=S+SGN(J)*BK(K1,J)*(XU**KJ-XV**KJ)/(KJ*C**J)
214           H=2*(-C/B)**K1*S/B
215           IF(U .LT. 0 .AND. V .LT. 0) H=-H
216          ENDIF
217         ENDIF
218        ELSE
219         IF(N .EQ. 1) THEN
220          ASSIGN 21 TO JMP2
221          GO TO 20
222    21    IF(K .EQ. -1) THEN
223           H=H
224          ELSEIF(K .EQ. -2) THEN
225           H=(RT(U)/U-RT(V)/V-HF*B*H)/A
226          ELSEIF(K .EQ. -3) THEN
227           H1=6*B
228           H2=4*A
229           H=((H1*V-H2)*RT(V)/V**2-(H1*U-H2)*RT(U)/U**2+
230      1       (3*B**2-H2*C)*H)/(8*A**2)
231          ENDIF
232         ELSE
233          IF(DELTA .EQ. 0) THEN
234           IF(C .EQ. 0) THEN
235            IF(K .EQ. -1) THEN
236             H=LOG(ABS(V/U))/SQRT(A)**3
237            ELSE
238             K1=K+1
239             H=(V**K1-U**K1)/(K1*SQRT(A)**3)
240            ENDIF
241           ELSE
242            X0=B/(2*C)
243            XV=1+X0/V
244            XU=1+X0/U
245            K1=-K+1
246            K2=-K-1
247            S=0
248            DO 3 J = 0,K1
249            KJ=K2-J
250            IF(KJ .NE. 0) THEN
251             S=S+SGN(J)*BK(K1,J)*(XV**KJ-XU**KJ)/KJ
252            ELSE
253             S=S+SGN(K2)*BK(K1,K2)*LOG(ABS(XV/XU))
254            ENDIF
255     3      CONTINUE
256            H=-S/(SQRT(C)**3*X0**(K1+1))
257            IF(U. LT. -X0) H=-H
258           ENDIF
259          ELSE
260           ASSIGN 22 TO JMP2
261           GO TO 20
262    22     IF(K .EQ. -1) THEN
263            H1=B*C
264            H2=B**2-2*A*C
265            H=(2*((H1*U+H2)/RT(U)-(H1*V+H2)/RT(V))/DELTA+H)/A
266           ELSEIF(K .EQ. -2) THEN
267            G1=3*B**2
268            G2=A*C
269            H1=(G1-8*G2)*C
270            H2=(G1-10*G2)*B
271            H3=A*DELTA
272            H=(((H1*V+H2-H3/V)/RT(V)-(H1*U+H2-H3/U)/RT(U))/DELTA
273      1       -3*HF*B*H)/A**2
274           ELSEIF(K .EQ. -3) THEN
275            G1=A*DELTA
276            G2=A*C
277            G3=B**2
278            G4=15*G3
279            H1=2*A*G1
280            H2=5*B*G1
281            H3=(G4-62*G2)*G3+24*G2**2
282            H4=B*C*(G4-52*G2)
283            H=((((H2-H1/V)/V-H3-H4*V)/RT(V)-((H2-H1/U)/U-H3-H4*U)/RT(U))
284      1        /DELTA+HF*(G4-12*G2)*H)/(4*A**3)
285           ENDIF
286          ENDIF
287         ENDIF
288        ENDIF
289       ENDIF
290       GO TO 9
291
292    10 C2=2*C
293       IF(DELTA .GT. 0 .OR. DELTA .LT. 0 .AND. C .GT. 0) THEN
294        H=LOG(ABS((2*SQRT(C*P(V))+C2*V+B)/
295      1           (2*SQRT(C*P(U))+C2*U+B)))/SQRT(C)
296       ELSEIF(DELTA .EQ. 0) THEN
297        H=ABS(LOG(ABS((C2*V+B)/(C2*U+B))))/SQRT(C)
298       ELSE
299        H=(ASIN((C2*U+B)/RTD)-ASIN((C2*V+B)/RTD))/SQRT(-C)
300       ENDIF
301       GO TO JMP1, (11,12,13)
302
303    20 IF(C .EQ. 0) THEN
304        IF(B .EQ. 0) THEN
305         H=LOG(ABS(V/U))/SQRT(A)
306        ELSE
307         IF(A .LT. 0) THEN
308          H=2*(ATAN(SQRT(-(A+B*V)/A))-ATAN(SQRT(-(A+B*U)/A)))/SQRT(-A)
309         ELSE
310          WA=SQRT(A)
311          WU=SQRT(A+B*U)
312          WV=SQRT(A+B*V)
313          H=LOG(ABS((WV-WA)*(WU+WA)/((WV+WA)*(WU-WA))))/WA
314         ENDIF
315        ENDIF
316       ELSE
317        A2=2*A
318        IF(DELTA .GT. 0 .OR. DELTA .LT. 0 .AND. A .GT. 0) THEN
319         H=LOG(ABS((-2*SQRT(A*P(V))+B*V+A2)*U/
320      1           ((-2*SQRT(A*P(U))+B*U+A2)*V)))/SQRT(A)
321        ELSEIF(DELTA .EQ. 0) THEN
322         H=LOG(ABS((B*U+A2)*V/((B*V+A2)*U)))/SQRT(A)
323         IF(U*V .GT. 0) H=SIGN(H,U)
324        ELSE
325         H=(ASIN((B*V+A2)/(V*RTD))-ASIN((B*U+A2)/(U*RTD)))/SQRT(-A)
326         IF(U .LT. 0 .AND. V .LT. 0) H=-H
327        ENDIF
328       ENDIF
329       GO TO JMP2, (21,22)
330
331     9 RES=SIGN(R1,V1-U1)*H
332       LRL=LLL
333       RETURN
334   102 FORMAT('ILLEGAL VALUE(S) K = ',I5,', N = ',I5)
335       END