3f50accb5fdc077df4d3d5305f6a350d227b592b
[scilab.git] / scilab / modules / elementary_functions / src / fortran / wacos.f
1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) Bruno Pincon
3
4 c This file must be used under the terms of the CeCILL.
5 c This source file is licensed as described in the file COPYING, which
6 c you should have received as part of this distribution.  The terms
7 c are also available at    
8 c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
9
10       subroutine wacos(zr, zi, ar, ai)
11 *
12 *     PURPOSE
13 *        Compute the arccosine of a complex number
14 *         a = ar + i ai = acos(z), z = zr + i zi
15 *       
16 *     CALLING LIST / PARAMETERS
17 *        subroutine wacos(zr,zi,ar,ai)
18 *        double precision zr,zi,ar,ai
19 *
20 *        zr,zi: real and imaginary parts of the complex number
21 *        ar,ai: real and imaginary parts of the result
22 *               ar,ai may have the same memory cases than zr et zi
23 *
24 *     REFERENCE
25 *        This is a Fortran-77 translation of an algorithm by 
26 *        T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which 
27 *        appears in their article :
28 *          "Implementing the Complex Arcsine and Arccosine 
29 *           Functions Using Exception Handling", ACM, TOMS, 
30 *           Vol 23, No. 3, Sept 1997, p. 299-335
31 *
32 *        with some modifications so as don't rely on ieee handle
33 *        trap functions.
34 *
35 *     AUTHOR
36 *        Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>
37 *        Thanks to Tom Fairgrieve
38 *
39       implicit none
40
41 *     PARAMETERS
42       double precision zr, zi, ar, ai
43
44 *     EXTERNAL FUNCTIONS
45       double precision dlamch, logp1
46       external         dlamch, logp1
47
48 *     CONSTANTS
49       double precision LN2, PI, HALFPI, Across, Bcross
50       parameter       (LN2    = 0.6931471805599453094172321d0,
51      $                 HALFPI = 1.5707963267948966192313216d0,
52      $                     PI = 3.1415926535897932384626433d0,
53      $                 Across = 1.5d0,
54      $                 Bcross = 0.6417d0)  
55 *     LOCAL VARIABLES
56       double precision x, y, A, B, R, S, Am1, szr, szi
57
58
59 *     STATIC VARIABLES
60       double precision LSUP, LINF, EPSM
61       save             LSUP, LINF, EPSM
62       logical          first
63       save             first
64       data             first /.true./
65
66 *     TEXT
67 *     got f.p. parameters used by the algorithm
68       if (first) then
69           LSUP = sqrt(dlamch('o'))/8.d0
70           LINF = 4.d0*sqrt(dlamch('u'))
71           EPSM = sqrt(dlamch('e'))
72           first = .false.
73       endif
74
75 *     avoid memory pb ...
76       x = abs(zr)
77       y = abs(zi)
78       szr = sign(1.d0,zr)
79       szi = sign(1.d0,zi)
80
81
82       if (LINF .le. min(x,y) .and. max(x,y) .le. LSUP ) then
83 *        we are in the safe region
84          R = sqrt((x+1.d0)**2 + y**2)
85          S = sqrt((x-1.d0)**2 + y**2)
86          A = 0.5d0*(R + S)
87          B = x/A
88
89 *        compute the real part
90          if ( B .le. Bcross ) then
91             ar = acos(B)
92          elseif ( x .le. 1.d0 ) then
93             ar = atan( sqrt( 0.5d0*(A+x) *
94      $                 ( (y**2)/(R+(x+1.d0)) + (S+(1.d0-x)) ) ) / x )
95          else
96             ar = atan((y*sqrt(0.5d0*((A+x)/(R+(x+1.d0))
97      $                  +(A+x)/(S+(x-1.d0))))) / x)
98          endif
99
100 *        compute the imaginary part
101          if ( A .le. Across ) then
102             if ( x .lt. 1.d0 ) then
103                Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x)))
104             else
105                Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0)))
106             endif
107             ai = logp1(Am1 + sqrt(Am1*(A+1.d0)))
108          else
109             ai = log(A + sqrt(A**2 - 1.d0))
110          endif
111
112       else
113 *        HANDLE BLOC : evaluation in the special regions ...
114          if ( y .le. EPSM*abs(x-1.d0) ) then
115             if (x .lt. 1.d0 ) then
116                ar = acos(x)
117                ai = y/sqrt((1.d0+x)*(1.d0-x))
118             else
119                ar = 0.d0
120                if ( x .le. LSUP ) then
121                   ai = logp1((x-1.d0) + sqrt((x-1.d0)*(x+1.d0)))
122                else
123                   ai = LN2 + log(x)
124                endif
125             endif
126          elseif (y .lt. LINF) then
127             ar = sqrt(y)
128             ai = ar
129          elseif (EPSM*y - 1.d0 .ge. x) then
130             ar = HALFPI
131             ai = LN2 + log(y)
132          elseif (x .gt. 1.d0) then
133             ar = atan(y/x)
134             ai = LN2 + log(y) + 0.5d0*logp1((x/y)**2)
135          else
136             ar = HALFPI
137             A = sqrt(1.d0 + y**2)
138             ai = 0.5d0*logp1(2.d0*y*(y+A))
139          endif
140       endif
141
142 *     recover the signs
143       if (szr .lt. 0.d0) then
144          ar = PI - ar
145       endif
146
147       if (y.ne.0.d0 .or. szr.lt.0.d0) then
148          ai = - szi * ai
149       endif
150
151       end
152
153