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