*Bug #13655 fixed - acos([2 %nan]) and asin([2 %nan]) returned wrong results (0 inste...
[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       integer          isanan
48       external         isanan
49
50 *     CONSTANTS
51       double precision LN2, PI, HALFPI, Across, Bcross
52       parameter       (LN2    = 0.6931471805599453094172321d0,
53      $                 HALFPI = 1.5707963267948966192313216d0,
54      $                     PI = 3.1415926535897932384626433d0,
55      $                 Across = 1.5d0,
56      $                 Bcross = 0.6417d0)  
57 *     LOCAL VARIABLES
58       double precision x, y, A, B, R, S, Am1, szr, szi
59
60
61 *     STATIC VARIABLES
62       double precision LSUP, LINF, EPSM
63       save             LSUP, LINF, EPSM
64       logical          first
65       save             first
66       data             first /.true./
67
68 *     TEXT
69 *     got f.p. parameters used by the algorithm
70       if (first) then
71           LSUP = sqrt(dlamch('o'))/8.d0
72           LINF = 4.d0*sqrt(dlamch('u'))
73           EPSM = sqrt(dlamch('e'))
74           first = .false.
75       endif
76
77 *     avoid memory pb ...
78       x = abs(zr)
79       y = abs(zi)
80       szr = sign(1.d0,zr)
81       szi = sign(1.d0,zi)
82             
83
84       if (LINF .le. min(x,y) .and. max(x,y) .le. LSUP ) then
85 *        we are in the safe region
86          R = sqrt((x+1.d0)**2 + y**2)
87          S = sqrt((x-1.d0)**2 + y**2)
88          A = 0.5d0*(R + S)
89          B = x/A
90
91 *        compute the real part
92          if ( B .le. Bcross ) then
93             ar = acos(B)
94          elseif ( x .le. 1.d0 ) then
95             ar = atan( sqrt( 0.5d0*(A+x) *
96      $                 ( (y**2)/(R+(x+1.d0)) + (S+(1.d0-x)) ) ) / x )
97          else
98             ar = atan((y*sqrt(0.5d0*((A+x)/(R+(x+1.d0))
99      $                  +(A+x)/(S+(x-1.d0))))) / x)
100          endif
101
102 *        compute the imaginary part
103          if ( A .le. Across ) then
104             if ( x .lt. 1.d0 ) then
105                Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x)))
106             else
107                Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0)))
108             endif
109             ai = logp1(Am1 + sqrt(Am1*(A+1.d0)))
110          else
111             ai = log(A + sqrt(A**2 - 1.d0))
112          endif
113
114       else
115 *        HANDLE BLOC : evaluation in the special regions ...
116          if ( y .le. EPSM*abs(x-1.d0) ) then
117             if (x .lt. 1.d0 ) then
118                ar = acos(x)
119                ai = y/sqrt((1.d0+x)*(1.d0-x))
120             else
121                ar = 0.d0
122                if ( x .le. LSUP ) then
123                   ai = logp1((x-1.d0) + sqrt((x-1.d0)*(x+1.d0)))
124                else
125                   ai = LN2 + log(x)
126                endif
127             endif
128          elseif (y .lt. LINF) then
129              if (isanan(x).eq.1) then
130                  ar = x
131                  ai = y
132              else
133                  ar = sqrt(y)
134                  ai = ar
135              endif
136          elseif (EPSM*y - 1.d0 .ge. x) then
137             ar = HALFPI
138             ai = LN2 + log(y)
139          elseif (x .gt. 1.d0) then
140             ar = atan(y/x)
141             ai = LN2 + log(y) + 0.5d0*logp1((x/y)**2)
142          else
143              if (isanan(x).eq.1) then
144                  ar = x
145              else
146                  ar = HALFPI
147              endif
148             A = sqrt(1.d0 + y**2)
149             ai = 0.5d0*logp1(2.d0*y*(y+A))
150          endif
151       endif
152
153 *     recover the signs
154       if (szr .lt. 0.d0) then
155          ar = PI - ar
156       endif
157
158       if (y.ne.0.d0 .or. szr.lt.0.d0) then
159          ai = - szi * ai
160       endif
161
162       end
163
164