e9d797ad0edd2bb9a5a993f1b89c4c0a11a14e30
[scilab.git] / scilab / modules / elementary_functions / src / fortran / wsqrt.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 wsqrt(xr,xi,yr,yi)
11 *
12 *     PURPOSE
13 *        wsqrt compute the square root of a complex number
14 *        y = yr + i yi = sqrt(x), x = xr + i xi
15 *
16 *     CALLING LIST / PARAMETERS
17 *        subroutine wsqrt(xr,xi,yr,yi)
18 *        double precision xr,xi,yr,yi
19 *
20 *        xr,xi: real and imaginary parts of the complex number
21 *        yr,yi: real and imaginary parts of the result
22 *               yr,yi may have the same memory cases than xr et xi
23 *
24 *     ALGORITHM
25 *        essentially the classic one which consists in
26 *        choosing the good formula such as avoid cancellation ; 
27 *        Also rare spurious overflow are treated with a
28 *        "manual" method. For some more "automated" methods
29 *        (but currently difficult to implement in a portable
30 *        way) see :
31 *
32 *          Hull, Fairgrieve, Tang,
33 *          "Implementing Complex Elementary Functions Using
34 *          Exception Handling", ACM TOMS, Vol. 20 (1994), pp 215-244
35 *
36 *        for xr > 0 :                            
37 *          yr = sqrt(2( xr + sqrt( xr^2 + xi^2)) )/ 2
38 *          yi = xi / sqrt(2(xr + sqrt(xr^2 + xi^2)))
39 *
40 *        and for xr < 0 :                           
41 *          yr = |xi| / sqrt( 2( -xr + sqrt( xr^2 + xi^2 )) )
42 *          yi = sign(xi) sqrt(2(-xr + sqrt( xr^2 + xi^2))) / 2
43 *     
44 *        for xr = 0 use          
45 *          yr = sqrt(0.5)*sqrt(|xi|)  when |xi| is such that 0.5*|xi| may underflow
46 *             = sqrt(0.5*|xi|)        else
47 *          yi = sign(xi) yr
48 *
49 *        Noting t = sqrt( 2( |xr| + sqrt( xr^2 + yr^2)) ) 
50 *                 = sqrt( 2( |xr| + pythag(xr,xi) ) )
51 *        it comes :
52 *
53 *          for xr > 0   |  for xr < 0
54 *         --------------+---------------------
55 *           yr = 0.5*t  |  yr =  |xi| / t
56 *           yi = xi / t |  yi = sign(xi)*0.5* t
57 *    
58 *        as the function pythag must not underflow (and overflow only 
59 *        if sqrt(x^2+y^2) > rmax) only spurious (rare) case of overflow 
60 *        occurs in which case a scaling is done.
61 *
62 *     AUTHOR
63 *        Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>
64 *
65       implicit none
66
67 *     PARAMETER
68       double precision xr, xi, yr, yi
69
70 *     LOCAL VAR
71       double precision a, b, t
72 *     EXTERNAL
73       double precision pythag, dlamch
74       external         pythag, dlamch 
75       integer          isanan
76       external         isanan
77
78 *     STATIC VAR
79       logical first
80         double precision RMAX, BRMIN
81       save    first
82       data    first /.true./
83       save             RMAX, BRMIN
84       
85
86       if (first) then
87          RMAX = dlamch('O')
88          BRMIN = 2.d0*dlamch('U')
89          first = .false.
90       endif
91
92       a = xr
93       b = xi
94
95       if ( a .eq. 0.d0 ) then
96 *        pure imaginary case
97          if ( abs(b) .ge. BRMIN ) then
98             yr = sqrt(0.5d0*abs(b))
99          else
100             yr = sqrt(abs(b))*sqrt(0.5d0)
101          endif
102          yi = sign(1.d0, b)*yr
103
104       elseif ( abs(a).le.RMAX .and. abs(b).le.RMAX ) then
105 *        standard case : a (not zero) and b are finite
106          t = sqrt(2.d0*(abs(a) + pythag(a,b)))
107 *        overflow test
108          if ( t .gt. RMAX ) goto 100
109 *        classic switch to get the stable formulas
110          if ( a .ge. 0.d0 ) then
111             yr = 0.5d0*t
112             yi = b/t
113          else
114             yr = abs(b)/t
115             yi = sign(0.5d0,b)*t
116          endif
117       else
118 *        Here we treat the special cases where a and b are +- 00 or NaN. 
119 *        The following is the treatment recommended by the C99 standard
120 *        with the simplification of returning NaN + i NaN if the
121 *        the real part or the imaginary part is NaN (C99 recommends
122 *        something more complicated)
123          if ( (isanan(a) .eq. 1) .or. (isanan(b) .eq. 1) ) then
124 *           got NaN + i NaN
125             yr = a + b
126             yi = yr
127          elseif ( abs(b) .gt. RMAX ) then
128 *           case a +- i oo -> result must be +oo +- i oo  for all a (finite or not)
129             yr = abs(b)
130             yi = b
131          elseif ( a .lt. -RMAX ) then
132 *           here a is -Inf and b is finite
133             yr = 0.d0
134             yi = sign(1.d0,b)*abs(a)
135          else
136 *           here a is +Inf and b is finite
137             yr = a
138             yi = 0.d0
139          endif
140       endif
141
142       return
143
144 *     handle (spurious) overflow by scaling a and b
145  100  continue
146       a = a/16.d0
147       b = b/16.d0
148       t = sqrt(2.d0*(abs(a) + pythag(a,b)))
149       if ( a .ge. 0.d0 ) then
150          yr = 2.d0*t
151          yi = 4.d0*b/t
152       else
153          yr = 4.d0*abs(b)/t
154          yi = sign(2.d0,b)*t
155       endif
156
157       end
158
159
160