[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / reli364.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:11  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
#if defined(CERNLIB_DOUBLE)
      FUNCTION DELI3(X,AKP,P)
C
#include "gen/imp64.inc"
C
      CHARACTER*(*) NAME
      PARAMETER(NAME='RELI3/DELI3')
#endif
#if !defined(CERNLIB_DOUBLE)
      FUNCTION RELI3(X,AKP,P)
C
      CHARACTER*(*) NAME
      PARAMETER(NAME='RELI3')
#endif
C
C     Translation of Algol procedure el3(x,kc,p) in
C      R. BULIRSCH Numerical Calculation of Elliptic Integrals and
C      Elliptic Functions III.,  Numer. Math. 13 (1969) 305-315
C
      LOGICAL LBO,LBK

      PARAMETER (ID = 16, IB = 4)
      PARAMETER (PI  = 3.14159 26535 89793 24D0)
      PARAMETER (AL2 = 0.69314 71805 59945 31D0)
      PARAMETER (ALB = IB*AL2, RLB = 1/ALB)
      PARAMETER (Z1 = 1, Z10 = 10, HF = Z1/2, C1 = Z1/10)
      PARAMETER (ND = ID-2, CA = Z10**(-ID/2), CB = Z10**(-(ID+2)))
      PARAMETER (ZD = HF/(ND+1))

      CHARACTER*80 ERRTXT

      DIMENSION RA(2:ND),RB(2:ND),RR(2:ND)

      DATA (RB(K),RA(K),K=2,ND)
     2/2.50000000000000000D-01, 7.50000000000000000D-01,
     3 1.66666666666666667D-01, 8.33333333333333333D-01,
     4 1.25000000000000000D-01, 8.75000000000000000D-01,
     5 1.00000000000000000D-01, 9.00000000000000000D-01,
     6 8.33333333333333333D-02, 9.16666666666666667D-01,
     7 7.14285714285714286D-02, 9.28571428571428571D-01,
     8 6.25000000000000000D-02, 9.37500000000000000D-01,
     9 5.55555555555555556D-02, 9.44444444444444445D-01,
     A 5.00000000000000000D-02, 9.50000000000000000D-01,
     B 4.54545454545454545D-02, 9.54545454545454545D-01,
     C 4.16666666666666667D-02, 9.58333333333333333D-01,
     D 3.84615384615384615D-02, 9.61538461538461538D-01,
     E 3.57142857142857143D-02, 9.64285714285714286D-01/

      IF(X .EQ. 0) THEN
       H=0
       GO TO 9
      ENDIF
      HH=X**2
      F=P*HH
      S=AKP
      IF(S .EQ. 0) S=CA/(1+ABS(X))
      T=S**2
      PM=HF*T
      E=HH*T
      Z=ABS(F)
      R=ABS(P)
      H1=HH+1
      IF(E .LT. C1 .AND. Z .LT. C1 .AND. T .LT. 1 .AND. R .LT. 1) THEN
       S=P+PM
       DO 1 K = 2,ND
       RR(K)=S
       PM=PM*T*RA(K)
    1  S=S*P+PM
       U=S*ZD
       S=U
       LBO=.FALSE.
       DO 2 K = ND,2,-1
       U=U+(RR(K)-U)*RB(K)
       LBO=.NOT.LBO
       V=U
       IF(LBO) V=-U
    2  S=S*HH+V
       IF(LBO) S=-S
       U=HF*(U+1)
#if defined(CERNLIB_DOUBLE)
       H=(U-S*H1)*SQRT(H1)*X+U*DASINH(X)
#endif
#if !defined(CERNLIB_DOUBLE)
       H=(U-S*H1)*SQRT(H1)*X+U*ASINH(X)
#endif
      ELSE
       W=1+F
       IF(W .EQ. 0) THEN
        H=0
        WRITE(ERRTXT,103) X,P
        CALL MTLPRT(NAME,'C346.3',ERRTXT)
        GO TO 9
       ENDIF
       PP=P
       IF(PP .EQ. 0) PP=CB/HH
       S=ABS(S)
       Y=ABS(X)
       G=PP-1
       IF(G .EQ. 0) G=CB
       F=PP-T
       IF(F .EQ. 0) F=CB*T
       AM=1-T
       AP=1+E
       R=PP*H1
       FA=G/(F*PP)
       LBO=FA .GT. 0
       FA=ABS(FA)
       PZ=ABS(G*F)
       DE=SQRT(PZ)
       Q=SQRT(ABS(PP))
       PM=PP-MIN(PM,HF)
       IF(PM .GE. 0) THEN
        U=SQRT(R*AP)
        V=Y*DE
        IF(G .LT. 0) V=-V
        D=1/Q
        C=1
       ELSE
        U=SQRT(H1*AP*PZ)
        YE=Y*Q
        V=AM*YE
        Q=-DE/G
        D=-AM/DE
        C=0
        PZ=AP-R
       ENDIF
       IF(LBO) THEN
        R=V/U
        Z=1
        K=1
        IF(PM .LT. 0) THEN
         H1=Y*SQRT(H1/(AP*FA))
         H1=1/H1-H1
         Z=H1-2*R
         R=2+R*H1
         IF(R .EQ. 0) R=CB
         IF(Z .EQ. 0) Z=H1*CB
         R=R/Z
         Z=R
         W=PZ
        ENDIF
        U=U/W
        V=V/W
       ELSE
        T=U+ABS(V)
        LBK=.TRUE.
        IF(PP .LT. 0) THEN
         DE=V/PZ
         YE=2*U*YE
         U=T/PZ
         V=-(F+G*E)/T
         T=PZ*ABS(W)
         Z=(HH*R*F-G*AP+YE)/T
         YE=YE/T
        ELSE
         DE=V/W
         YE=0
         U=(E+PP)/T
         V=T/W
         Z=1
        ENDIF
        IF(S .GT. 1) THEN
         H1=U
         U=V
         V=H1
        ENDIF
       ENDIF
       Y=1/Y
       E=S
       T=1
       N=1
       L=0
       M=0
    3  Y=Y-E/Y
       IF(Y .EQ. 0) Y =SQRT(E)*CB
       F=C
       C=D/Q+C
       G=E/Q
       D=2*(F*G+D)
       Q=G+Q
       G=T
       T=S+T
       N=2*N
       M=2*M
       IF(LBO) THEN
        IF(Z .LT. 0) M=K+M
        K=0
        IF(R .GT. 0) K=1
        IF(R .LT. 0) K=-1
        H1=E/(U**2+V**2)
        U=U*(1+H1)
        V=V*(1-H1)
       ELSE
        R=U/V
        H1=Z*R
        Z=H1*Z
        HH=E/V
        IF(LBK) THEN
         DE=DE/U
         YE=YE*(H1+1/H1)+DE*(1+R)
         DE=DE*(U-HH)
         LBK=ABS(YE) .LT. 1
        ELSE
         A=LOG(Z)
         K=RLB*A
         IF(Z .GE. 1) K=K+1
         Z=EXP(A-K*ALB)
         M=M+K
        ENDIF
       ENDIF
       IF(ABS(G-S) .GT. CA*G) THEN
        IF(LBO) THEN
         G=HF*(1/R-R)
         HH=U+V*G
         H1=G*U-V
         IF(HH .EQ. 0) HH=U*CB
         IF(H1 .EQ. 0) H1=V*CB
         Z=R*H1
         R=HH/H1
        ELSE
         U=U+E/U
         V=V+HH
        ENDIF
        S=2*SQRT(E)
        E=S*T
        L=2*L
        IF(Y .LT. 0) L=L+1
        GO TO 3
       ENDIF
       IF(Y .LT. 0) L=L+1
       E=(ATAN(T/Y)+PI*L)*(C*T+D)/(T*(T+Q))
       IF(LBO) THEN
        H1=V/(T+U)
        Z=1-R*H1
        H1=R+H1
        IF(Z .EQ. 0) Z=CB
        IF(Z .LT. 0) THEN
         IF(H1 .GT. 0) M=M+1
         IF(H1 .LT. 0) M=M-1
        ENDIF
        S=ATAN(H1/Z)+M*PI
       ELSE
        IF(LBK) THEN
#if defined(CERNLIB_DOUBLE)
         S=HF*DASINH(YE)
#endif
#if !defined(CERNLIB_DOUBLE)
         S=HF*ASINH(YE)
#endif
        ELSE
         S=HF*(LOG(Z)+M*ALB)
        ENDIF
       ENDIF
       H=(E+SQRT(FA)*S)/N
       IF(X .LT. 0) H=-H
      ENDIF
#if defined(CERNLIB_DOUBLE)
    9 DELI3=H
#endif
#if !defined(CERNLIB_DOUBLE)
    9 RELI3=H
#endif
      RETURN
  103 FORMAT('FUNCTION SINGULAR FOR ',1P,
     1 'X = ',D15.8,3X,'P = ',D15.8,4X,'(P*X**2 = -1)')
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1996/04/01 15:02:11 mclareni
	6	* Mathlib gen
	7	*
	8	*
	9	#include "gen/pilot.h"
	10	#if defined(CERNLIB_DOUBLE)
	11	FUNCTION DELI3(X,AKP,P)
	12	C
	13	#include "gen/imp64.inc"
	14	C
	15	CHARACTER() NAME
	16	PARAMETER(NAME='RELI3/DELI3')
	17	#endif
	18	#if !defined(CERNLIB_DOUBLE)
	19	FUNCTION RELI3(X,AKP,P)
	20	C
	21	CHARACTER() NAME
	22	PARAMETER(NAME='RELI3')
	23	#endif
	24	C
	25	C Translation of Algol procedure el3(x,kc,p) in
	26	C R. BULIRSCH Numerical Calculation of Elliptic Integrals and
	27	C Elliptic Functions III., Numer. Math. 13 (1969) 305-315
	28	C
	29	LOGICAL LBO,LBK
	30
	31	PARAMETER (ID = 16, IB = 4)
	32	PARAMETER (PI = 3.14159 26535 89793 24D0)
	33	PARAMETER (AL2 = 0.69314 71805 59945 31D0)
	34	PARAMETER (ALB = IB*AL2, RLB = 1/ALB)
	35	PARAMETER (Z1 = 1, Z10 = 10, HF = Z1/2, C1 = Z1/10)
	36	PARAMETER (ND = ID-2, CA = Z10(-ID/2), CB = Z10(-(ID+2)))
	37	PARAMETER (ZD = HF/(ND+1))
	38
	39	CHARACTER*80 ERRTXT
	40
	41	DIMENSION RA(2:ND),RB(2:ND),RR(2:ND)
	42
	43	DATA (RB(K),RA(K),K=2,ND)
	44	2/2.50000000000000000D-01, 7.50000000000000000D-01,
	45	3 1.66666666666666667D-01, 8.33333333333333333D-01,
	46	4 1.25000000000000000D-01, 8.75000000000000000D-01,
	47	5 1.00000000000000000D-01, 9.00000000000000000D-01,
	48	6 8.33333333333333333D-02, 9.16666666666666667D-01,
	49	7 7.14285714285714286D-02, 9.28571428571428571D-01,
	50	8 6.25000000000000000D-02, 9.37500000000000000D-01,
	51	9 5.55555555555555556D-02, 9.44444444444444445D-01,
	52	A 5.00000000000000000D-02, 9.50000000000000000D-01,
	53	B 4.54545454545454545D-02, 9.54545454545454545D-01,
	54	C 4.16666666666666667D-02, 9.58333333333333333D-01,
	55	D 3.84615384615384615D-02, 9.61538461538461538D-01,
	56	E 3.57142857142857143D-02, 9.64285714285714286D-01/
	57
	58	IF(X .EQ. 0) THEN
	59	H=0
	60	GO TO 9
	61	ENDIF
	62	HH=X**2
	63	F=P*HH
	64	S=AKP
65	IF(S .EQ. 0) S=CA/(1+ABS(X))
66	T=S**2
67	PM=HF*T
68	E=HH*T
69	Z=ABS(F)
70	R=ABS(P)
71	H1=HH+1
72	IF(E .LT. C1 .AND. Z .LT. C1 .AND. T .LT. 1 .AND. R .LT. 1) THEN
73	S=P+PM
74	DO 1 K = 2,ND
75	RR(K)=S
76	PM=PMTRA(K)
77	1 S=S*P+PM
78	U=S*ZD
79	S=U
80	LBO=.FALSE.
81	DO 2 K = ND,2,-1
82	U=U+(RR(K)-U)*RB(K)
83	LBO=.NOT.LBO
84	V=U
85	IF(LBO) V=-U
86	2 S=S*HH+V
87	IF(LBO) S=-S
88	U=HF*(U+1)
89	#if defined(CERNLIB_DOUBLE)
90	H=(U-SH1)SQRT(H1)X+UDASINH(X)
91	#endif
92	#if !defined(CERNLIB_DOUBLE)
93	H=(U-SH1)SQRT(H1)X+UASINH(X)
94	#endif
95	ELSE
96	W=1+F
97	IF(W .EQ. 0) THEN
98	H=0
99	WRITE(ERRTXT,103) X,P
100	CALL MTLPRT(NAME,'C346.3',ERRTXT)
101	GO TO 9
102	ENDIF
103	PP=P
104	IF(PP .EQ. 0) PP=CB/HH
105	S=ABS(S)
106	Y=ABS(X)
107	G=PP-1
108	IF(G .EQ. 0) G=CB
109	F=PP-T
110	IF(F .EQ. 0) F=CB*T
111	AM=1-T
112	AP=1+E
113	R=PP*H1
114	FA=G/(F*PP)
115	LBO=FA .GT. 0
116	FA=ABS(FA)
117	PZ=ABS(G*F)
118	DE=SQRT(PZ)
119	Q=SQRT(ABS(PP))
120	PM=PP-MIN(PM,HF)
121	IF(PM .GE. 0) THEN
122	U=SQRT(R*AP)
123	V=Y*DE
124	IF(G .LT. 0) V=-V
125	D=1/Q
126	C=1
127	ELSE
128	U=SQRT(H1APPZ)
129	YE=Y*Q
130	V=AM*YE
131	Q=-DE/G
132	D=-AM/DE
133	C=0
134	PZ=AP-R
135	ENDIF
136	IF(LBO) THEN
137	R=V/U
138	Z=1
139	K=1
140	IF(PM .LT. 0) THEN
141	H1=YSQRT(H1/(APFA))
142	H1=1/H1-H1
143	Z=H1-2*R
144	R=2+R*H1
145	IF(R .EQ. 0) R=CB
146	IF(Z .EQ. 0) Z=H1*CB
147	R=R/Z
148	Z=R
149	W=PZ
150	ENDIF
151	U=U/W
152	V=V/W
153	ELSE
154	T=U+ABS(V)
155	LBK=.TRUE.
156	IF(PP .LT. 0) THEN
157	DE=V/PZ
158	YE=2UYE
159	U=T/PZ
160	V=-(F+G*E)/T
161	T=PZ*ABS(W)
162	Z=(HHRF-G*AP+YE)/T
163	YE=YE/T
164	ELSE
165	DE=V/W
166	YE=0
167	U=(E+PP)/T
168	V=T/W
169	Z=1
170	ENDIF
171	IF(S .GT. 1) THEN
172	H1=U
173	U=V
174	V=H1
175	ENDIF
176	ENDIF
177	Y=1/Y
178	E=S
179	T=1
180	N=1
181	L=0
182	M=0
183	3 Y=Y-E/Y
184	IF(Y .EQ. 0) Y =SQRT(E)*CB
185	F=C
186	C=D/Q+C
187	G=E/Q
188	D=2(FG+D)
189	Q=G+Q
190	G=T
191	T=S+T
192	N=2*N
193	M=2*M
194	IF(LBO) THEN
195	IF(Z .LT. 0) M=K+M
196	K=0
197	IF(R .GT. 0) K=1
198	IF(R .LT. 0) K=-1
199	H1=E/(U2+V2)
200	U=U*(1+H1)
201	V=V*(1-H1)
202	ELSE
203	R=U/V
204	H1=Z*R
205	Z=H1*Z
206	HH=E/V
207	IF(LBK) THEN
208	DE=DE/U
209	YE=YE(H1+1/H1)+DE(1+R)
210	DE=DE*(U-HH)
211	LBK=ABS(YE) .LT. 1
212	ELSE
213	A=LOG(Z)
214	K=RLB*A
215	IF(Z .GE. 1) K=K+1
216	Z=EXP(A-K*ALB)
217	M=M+K
218	ENDIF
219	ENDIF
220	IF(ABS(G-S) .GT. CA*G) THEN
221	IF(LBO) THEN
222	G=HF*(1/R-R)
223	HH=U+V*G
224	H1=G*U-V
225	IF(HH .EQ. 0) HH=U*CB
226	IF(H1 .EQ. 0) H1=V*CB
227	Z=R*H1
228	R=HH/H1
229	ELSE
230	U=U+E/U
231	V=V+HH
232	ENDIF
233	S=2*SQRT(E)
234	E=S*T
235	L=2*L
236	IF(Y .LT. 0) L=L+1
237	GO TO 3
238	ENDIF
239	IF(Y .LT. 0) L=L+1
240	E=(ATAN(T/Y)+PIL)(CT+D)/(T(T+Q))
241	IF(LBO) THEN
242	H1=V/(T+U)
243	Z=1-R*H1
244	H1=R+H1
245	IF(Z .EQ. 0) Z=CB
246	IF(Z .LT. 0) THEN
247	IF(H1 .GT. 0) M=M+1
248	IF(H1 .LT. 0) M=M-1
249	ENDIF
250	S=ATAN(H1/Z)+M*PI
251	ELSE
252	IF(LBK) THEN
253	#if defined(CERNLIB_DOUBLE)
254	S=HF*DASINH(YE)
255	#endif
256	#if !defined(CERNLIB_DOUBLE)
257	S=HF*ASINH(YE)
258	#endif
259	ELSE
260	S=HF(LOG(Z)+MALB)
261	ENDIF
262	ENDIF
263	H=(E+SQRT(FA)*S)/N
264	IF(X .LT. 0) H=-H
265	ENDIF
266	#if defined(CERNLIB_DOUBLE)
267	9 DELI3=H
268	#endif
269	#if !defined(CERNLIB_DOUBLE)
270	9 RELI3=H
271	#endif
272	RETURN
273	103 FORMAT('FUNCTION SINGULAR FOR ',1P,
274	1 'X = ',D15.8,3X,'P = ',D15.8,4X,'(PX*2 = -1)')
275	END