This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / elfun64.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:01  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10 #if defined(CERNLIB_DOUBLE)
11       SUBROUTINE DELFUN(U,AK2,SN,CN,DN)
12 #include "gen/imp64.inc"
13 #endif
14 #if defined(CERNLIB_SINGLE)
15       SUBROUTINE RELFUN(U,AK2,SN,CN,DN)
16 #endif
17       DIMENSION C(4)
18
19 C     Machine-dependent: EPS1=2**-(MB/2), EPS2=2**-(MB+3)
20 C     Where M = Number of bits in mantissa
21
22       PARAMETER (MB = 64)
23       PARAMETER (Z0 = 0, Z1 = 1, Z2 = 2, HF = Z1/2, QU = Z1/4)
24       PARAMETER (PI = 3.14159 26535 89793 24D0, PIH = PI/2)
25       PARAMETER (EPS1 = Z2**(-MB/2), EPS2 = Z2**(-(MB+3)))
26
27       DATA AM0 /-1D20/
28
29       SAVE AM0,C,A,L,BIGK
30
31 #if defined(CERNLIB_SINGLE)
32       ENTRY ELFUN(U,AK2,SN,CN,DN)
33 #endif
34
35       XM=ABS(AK2)
36       IF(U .EQ. 0) THEN
37        SN=0
38        DN=1
39        CN=1
40       ELSEIF(XM .EQ. 0) THEN
41        SN=SIN(U)
42        DN=1
43        CN=COS(U)
44       ELSEIF(XM .EQ. 1) THEN
45        SN=TANH(U)
46        DN=1/COSH(U)
47        CN=DN
48       ELSE
49        IF(XM .LE. 1) THEN
50         U1=U
51         AM=XM
52        ELSE
53         W=SQRT(XM)
54         U1=W*U
55         AM=1/XM
56        ENDIF
57        IF(AM .LE. HF) THEN
58         IF(AM .EQ. AM0) GO TO 1
59         AM0=AM
60         C(4)=QU*AM
61         B=SQRT(1-AM)
62        ELSE
63         AMC=1-AM
64         IF(AMC .EQ. AM0) GO TO 1
65         AM0=AMC
66         C(4)=QU*AMC
67         B=SQRT(AM)
68        ENDIF
69
70 C     Gauss arithmetic-geometric mean. Skipped if previous modulus.
71
72        A=1
73        L=4
74     2  IF(C(L) .GE. EPS1) THEN
75         L=L-1
76         C(L)=(QU*(A-B))**2
77         A1=HF*(A+B)
78         B=SQRT(A*B)
79         A=A1
80         GO TO 2
81        ENDIF
82        A=HF*(A+B)
83        BIGK=PIH/A
84
85 C     Descending Landen-Gauss transformation for real argument
86
87     1  IF(AM .LE. HF) THEN
88         X=SIN(A*U1)
89         IF(X .EQ. 0) THEN
90          SN=0
91          DN=1
92          CN=1
93         ELSE
94          X=A/X
95          DO 3 J = L,4
96          X1=C(J)/X
97     3    X=X1+X
98          H=1/X
99          SN=H
100          DN=1-2*X1*H
101          CN=SIGN(SQRT(ABS(1-H**2)),BIGK-ABS(U1))
102         ENDIF
103        ELSE
104
105 C     Descending Landen-Gauss Transformation for imaginary argument
106
107         Y=A/SINH(A*U1)
108         DO 4 J=L,4
109         Y1=C(J)/Y
110         Y=Y-Y1
111     4   IF(Y .EQ. 0) Y=EPS2
112         H=1/Y
113         Y1=2*Y1*H
114         CN=SIGN(SQRT(Y/(H+Y)),Y1)
115         DN=CN*(1+Y1)
116         SN=CN*H
117        ENDIF
118        IF(XM .GT. 1) THEN
119         SN=SN/W
120         H=DN
121         DN=CN
122         CN=H
123        ENDIF
124       ENDIF
125       RETURN
126       END