1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2c Copyright (C) INRIA
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.
12c
13      subroutine wasin(zr, zi, ar, ai)
14*
15*     PURPOSE
16*        Compute the arcsin of a complex number
17*         a = ar + i ai = asin(z), z = zr + i zi
18*
19*     CALLING LIST / PARAMETERS
20*        subroutine wasin(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, HALFPI, Across, Bcross
55      parameter       (LN2    = 0.6931471805599453094172321d0,
56     $                 HALFPI = 1.5707963267948966192313216d0,
57     $                 Across = 1.5d0,
58     $                 Bcross = 0.6417d0)
59
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*     TEXT
71*     got f.p. parameters used by the algorithm
72      if (first) then
73          LSUP = sqrt(dlamch('o'))/8.d0
74          LINF = 4.d0*sqrt(dlamch('u'))
75          EPSM = sqrt(dlamch('e'))
76          first = .false.
77      endif
78
79*     avoid memory pb ...
80      x = abs(zr)
81      y = abs(zi)
82      szr = sign(1.d0,zr)
83      szi = sign(1.d0,zi)
84
85
86      if (LINF .le. min(x,y) .and. max(x,y) .le. LSUP ) then
87*        we are in the safe region
88         R = sqrt((x+1.d0)**2 + y**2)
89         S = sqrt((x-1.d0)**2 + y**2)
90         A = 0.5d0*(R + S)
91         B = x/A
92
93*        compute the real part
94         if ( B .le. Bcross ) then
95            ar = asin(B)
96         elseif ( x .le. 1.d0 ) then
97            ar = atan( x / sqrt(
98     $              0.5d0*(A+x)*((y**2)/(R+(x+1.d0))+(S+(1.d0-x)))) )
99         else
100            ar = atan( x /
101     $          (y*sqrt(0.5d0*((A+x)/(R+(x+1.d0))+(A+x)/(S+(x-1.d0))))))
102         endif
103
104*        compute the imaginary part
105         if ( A .le. Across ) then
106            if ( x .lt. 1.d0 ) then
107               Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x)))
108            else
109               Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0)))
110            endif
111            ai = logp1(Am1 + sqrt(Am1*(A+1.d0)))
112         else
113            ai = log(A + sqrt(A**2 - 1.d0))
114         endif
115
116      else
117*        HANDLE BLOC : evaluation in the special regions ...
118         if ( y .le. EPSM*abs(x-1.d0) ) then
119            if (x .lt. 1.d0 ) then
120               ar = asin(x)
121               ai = y/sqrt((1.d0+x)*(1.d0-x))
122            else
123               ar = HALFPI
124               if ( x .le. LSUP ) then
125                  ai = logp1((x-1.d0) + sqrt((x-1.d0)*(x+1.d0)))
126               else
127                  ai = LN2 + log(x)
128               endif
129            endif
130
131         elseif (y .lt. LINF) then
132            if (isanan(x).eq.1) then
133               ar = x
134            else
135               ar = HALFPI - sqrt(y)
136            endif
137            ai = sqrt(y)
138
139         elseif (EPSM*y - 1.d0 .ge. x) then
140            ar = x/y
141            ai = LN2 + log(y)
142
143         elseif (x .gt. 1.d0) then
144            ar = atan(x/y)
145            ai = LN2 + log(y) + 0.5d0*logp1((x/y)**2)
146
147         else
148            A = sqrt(1 + y**2)
149            ar = x/A
150            ai = 0.5*logp1(2.d0*y*(y+A))
151         endif
152      endif
153
154*     recover the signs
155      ar = szr * ar
156      if (y.eq.0d00 .and. szr.gt.0d00) then
157          szi = - szi
158      endif
159      ai = szi * ai
160
161      end
162
163
164