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