1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) Bruno Pincon
3 c
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
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
67 *     PARAMETER
68       double precision xr, xi, yr, yi
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
78 *     STATIC VAR
79       logical first
80         double precision RMAX, BRMIN
81       save    first
82       data    first /.true./
83       save             RMAX, BRMIN
86       if (first) then
87          RMAX = dlamch('O')
88          BRMIN = 2.d0*dlamch('U')
89          first = .false.
90       endif
92       a = xr
93       b = xi
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
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             if ( b .ge. 0 ) then
116                yi = 0.5d0*t
117             else
118                yi = -0.5d0*t
119             endif
120          endif
121       else
122 *        Here we treat the special cases where a and b are +- 00 or NaN.
123 *        The following is the treatment recommended by the C99 standard
124 *        with the simplification of returning NaN + i NaN if the
125 *        the real part or the imaginary part is NaN (C99 recommends
126 *        something more complicated)
127          if ( (isanan(a) .eq. 1) .or. (isanan(b) .eq. 1) ) then
128 *           got NaN + i NaN
129             yr = a + b
130             yi = yr
131          elseif ( abs(b) .gt. RMAX ) then
132 *           case a +- i oo -> result must be +oo +- i oo  for all a (finite or not)
133             yr = abs(b)
134             yi = b
135          elseif ( a .lt. -RMAX ) then
136 *           here a is -Inf and b is finite
137             yr = 0.d0
138             if ( b .ge. 0 ) then
139                yi = 1.d0*abs(a)
140             else
141                yi = -1.d0*abs(a)
142             endif
143          else
144 *           here a is +Inf and b is finite
145             yr = a
146             yi = 0.d0
147          endif
148       endif
150       return
152 *     handle (spurious) overflow by scaling a and b
153  100  continue
154       a = a/16.d0
155       b = b/16.d0
156       t = sqrt(2.d0*(abs(a) + pythag(a,b)))
157       if ( a .ge. 0.d0 ) then
158          yr = 2.d0*t
159          yi = 4.d0*b/t
160       else
161          yr = 4.d0*abs(b)/t
162          if ( b .ge. 0 ) then
163             yi = 2.d0*t
164          else
165             yi = -2.d0*t
166          endif
167       endif
169       end