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