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
|
subroutine wlog(xr,xi,yr,yi)
*
* PURPOSE
* wlog compute the logarithm of a complex number
* y = yr + i yi = log(x), x = xr + i xi
*
* CALLING LIST / PARAMETERS
* subroutine wlog(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
*
* METHOD
* adapted with some modifications from Hull,
* Fairgrieve, Tang, "Implementing Complex
* Elementary Functions Using Exception Handling",
* ACM TOMS, Vol. 20 (1994), pp 215-244
*
* y = yr + i yi = log(x)
* yr = log(|x|) = various formulae depending where x is ...
* yi = Arg(x) = atan2(xi, xr)
*
* 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, r
* CONSTANTS
double precision R2
parameter (R2 = 1.41421356237309504d0)
* EXTERNAL
double precision dlamch, logp1, pythag
external dlamch, logp1, pythag
* STATIC VAR
logical first
double precision RMAX, LSUP, LINF
save first
data first /.true./
save RMAX, LSUP, LINF
if (first) then
RMAX = dlamch('O')
LINF = sqrt(dlamch('U'))
LSUP = sqrt(0.5d0*RMAX)
first = .false.
endif
* (0) avoid memory pb ...
a = xr
b = xi
* (1) compute the imaginary part
yi = atan2(b, a)
* (2) compute the real part
a = abs(a)
b = abs(b)
* Order a and b such that 0 <= b <= a
if (b .gt. a) then
t = b
b = a
a = t
endif
if ( (0.5d0 .le. a) .and. (a .le. R2) ) then
yr = 0.5d0*logp1((a-1.d0)*(a+1.d0) + b*b)
elseif (LINF .lt. b .and. a .lt. LSUP) then
* no overflow or underflow can occur in computing a*a + b*b
yr = 0.5d0*log(a*a + b*b)
elseif (a .gt. RMAX) then
* overflow
yr = a
else
t = pythag(a,b)
if (t .le. RMAX) then
yr = log(t)
else
* handle rare spurious overflow with :
r = b/a
yr = log(a) + 0.5d0*logp1(r*r)
endif
endif
end
|
