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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
|
/*
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
* Copyright (C) 2006-2008 - INRIA - Bruno JOFRET
* Copyright (C) Bruno Pincon
*
* This file must be used under the terms of the CeCILL.
* This source file is licensed as described in the file COPYING, which
* you should have received as part of this distribution. The terms
* are also available at
* http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
*
*/
/*
PURPOSE
watan compute the arctangent of a complex number
y = yr + i yi = atan(x), x = xr + i xi
CALLING LIST / PARAMETERS
subroutine watan(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
COPYRIGHT (C) 2001 Bruno Pincon and Lydia van Dijk
Written by Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr> so
as to get more precision. Also to fix the
behavior at the singular points and at the branch cuts.
Polished by Lydia van Dijk
<lvandijk@hammersmith-consulting.com>
CHANGES : - (Bruno on 2001 May 22) for ysptrk use a
minimax polynome to enlarge the special
evaluation zone |s| < SLIM. Also rename
this function as lnp1m1.
- (Bruno on 2001 June 7) better handling
of spurious over/underflow ; remove
the call to pythag ; better accuracy
in the real part for z near +-i
EXTERNALS FUNCTIONS
dlamch
lnp1m1 (at the end of this file)
ALGORITHM : noting z = a + i*b, we have:
Z = yr + yi*b = arctan(z) = (i/2) * log( (i+z)/(i-z) )
This function has two branch points at +i and -i and the
chosen branch cuts are the two half-straight lines
D1 = [i, i*oo) and D2 = (-i*oo, i]. The function is then
analytic in C \ (D1 U D2)).
From the definition it follows that:
yr = 0.5 Arg ( (i+z)/(i-z) ) (1)
yi = 0.5 log (|(i+z)/(i-z)|) (2)
so lim (z -> +- i) yr = undefined (and Nan is logical)
lim (z -> +i) yi = +oo
lim (z -> -i) yi = -oo
The real part of arctan(z) is discontinuous across D1 and D2
and we impose the following definitions:
if imag(z) > 1 then
Arg(arctan(z)) = pi/2 (=lim real(z) -> 0+)
if imag(z) < 1 then
Arg(arctan(z)) = -pi/2 (=lim real(z) -> 0-)
Basic evaluation: if we write (i+z)/(i-z) using
z = a + i*b, we get:
i+z 1-(a**2+b**2) + i*(2a)
--- = ----------------------
i-z a**2 + (1-b)**2
then, with r2 = |z|^2 = a**2 + b**2 :
yr = 0.5 * Arg(1-r2 + (2*a)*i)
= 0.5 * atan2(2a, (1-r2)) (3)
This formula is changed when r2 > RMAX (max pos float)
and also when |1-r2| and |a| are near 0 (see comments
in the code).
After some math:
yi = 0.25 * log( (a**2 + (b + 1)**2) /
(a**2 + (b - 1)**2) ) (4)
Evaluation for "big" |z|
------------------------
If |z| is "big", the direct evaluation of yi by (4) may
suffer of innaccuracies and of spurious overflow. Noting
that s = 2 b / (1 + |z|**2), we have:
yi = 0.25 log ( (1 + s)/(1 - s) ) (5)
3 5
yi = 0.25*( 2 * ( s + 1/3 s + 1/5 s + ... ))
yi = 0.25 * lnp1m1(s) if |s| < SLIM
So if |s| is less than SLIM we switch to a special
evaluation done by the function lnp1m1. The
threshold value SLIM is choosen by experiment
(with the Pari-gp software). For |s|
"very small" we used a truncated taylor dvp,
else a minimax polynome (see lnp1m1).
To avoid spurious overflows (which result in spurious
underflows for s) in computing s with s= 2 b / (1 + |z|**2)
when |z|^2 > RMAX (max positive float) we use :
s = 2d0 / ( (a/b)*a + b )
but if |b| = Inf this formula leads to NaN when
|a| is also Inf. As we have :
|s| <= 2 / |b|
we impose simply : s = 0 when |b| = Inf
Evaluation for z very near to i or -i:
--------------------------------------
Floating point numbers of the form a+i or a-i with 0 <
a**2 < tiny (approximately 1d-308) may lead to underflow
(i.e., a**2 = 0) and the logarithm will break formula (4).
So we switch to the following formulas:
If b = +-1 and |a| < sqrt(tiny) approximately 1d-150 (say)
then (by using that a**2 + 4 = 4 in machine for such a):
yi = 0.5 * log( 2/|a| ) for b=1
yi = 0.5 * log( |a|/2 ) for b=-1
finally: yi = 0.5 * sign(b) * log( 2/|a| )
yi = 0.5 * sign(b) * (log(2) - log(|a|)) (6)
The last trick is to avoid overflow for |a|=tiny! In fact
this formula may be used until a**2 + 4 = 4 so that the
threshold value may be larger.
*/
#include <math.h>
#include "lapack.h"
#include "atan.h"
#include "abs.h"
#include "lnp1m1.h"
#define _sign(a, b) b >=0 ? a : -a
doubleComplex zatans(doubleComplex z) {
static double sSlim = 0.2;
static double sAlim = 1E-150;
static double sTol = 0.3;
static double sLn2 = 0.6931471805599453094172321;
double RMax = getOverflowThreshold();
double Pi_2 = 2.0 * datans(1);
double _inReal = zreals(z);
double _inImg = zimags(z);
double _outReal = 0;
double _outImg = 0;
/* Temporary variables */
double R2 = 0;
double S = 0;
if(_inImg == 0)
{
_outReal = datans(_inReal);
_outImg = 0;
}
else
{
R2 = _inReal * _inReal + _inImg * _inImg; /* Oo */
if(R2 > RMax)
{
if( dabss(_inImg) > RMax)
S = 0;
else
S = 1 / (((0.5 * _inReal) / _inImg) * _inReal + 0.5 * _inImg );
}
else
S = (2 * _inImg) / (1+R2);
if(dabss(S) < sSlim)
{
/*
s is small: |s| < SLIM <=> |z| outside the following disks:
D+ = D(center = [0; 1/slim], radius = sqrt(1/slim**2 - 1)) if b > 0
D- = D(center = [0; -1/slim], radius = sqrt(1/slim**2 - 1)) if b < 0
use the special evaluation of log((1+s)/(1-s)) (5)
*/
_outImg = dlnp1m1s(S) * 0.25;
}
else
{
if(dabss(S) == 1 && dabss(_inReal) <= sAlim)
{
/* |s| >= SLIM => |z| is inside D+ or D- */
_outImg = _sign(0.5,_inImg) * ( sLn2 - log(dabss(_inReal)));
}
else
{
_outImg = 0.25 * log((pow(_inReal,2) + pow((_inImg + 1),2)) / (pow(_inReal,2) + pow((_inImg - 1),2)));
}
}
if(_inReal == 0)
{/* z is purely imaginary */
if( dabss(_inImg) > 1)
{/* got sign(b) * pi/2 */
_outReal = _sign(1, _inImg) * Pi_2;
}
else if( dabss(_inImg) == 1)
{/* got a Nan with 0/0 */
_outReal = (_inReal - _inReal) / (_inReal - _inReal); /* Oo */
}
else
_outReal = 0;
}
else if(R2 > RMax)
{/* _outImg is necessarily very near sign(a)* pi/2 */
_outReal = _sign(1, _inReal) * Pi_2;
}
else if(dabss(1 - R2) + dabss(_inReal) <= sTol)
{/* |b| is very near 1 (and a is near 0) some cancellation occur in the (next) generic formula */
_outReal = 0.5 * atan2(2 * _inReal, (1-_inImg) * (1 + _inImg) - pow(_inReal,2));
}
else
_outReal = 0.5 * atan2(2 * _inReal, 1 - R2);
}
return DoubleComplex(_outReal, _outImg);
}
|
