diff options
Diffstat (limited to '2.3-1/src/c/elementaryFunctions/atan/catans.c')
| -rw-r--r-- | 2.3-1/src/c/elementaryFunctions/atan/catans.c | 249 |
1 files changed, 249 insertions, 0 deletions
diff --git a/2.3-1/src/c/elementaryFunctions/atan/catans.c b/2.3-1/src/c/elementaryFunctions/atan/catans.c new file mode 100644 index 00000000..d2081181 --- /dev/null +++ b/2.3-1/src/c/elementaryFunctions/atan/catans.c @@ -0,0 +1,249 @@ +/*
+ * 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 "atan.h"
+#include "abs.h"
+#include "lnp1m1.h"
+#include "lapack.h"
+
+#define _sign(a, b) b >=0 ? a : -a
+
+floatComplex catans(floatComplex z) {
+ static float sSlim = 0.2f;
+ /* .
+ ** / \ WARNING : this algorithm was based on double precision
+ ** / ! \ using float truncate the value to 0.
+ ** `----'
+ **
+ ** static float sAlim = 1E-150f;
+ */
+ static float sAlim = 0.0f;
+ static float sTol = 0.3f;
+ static float sLn2 = 0.6931471805599453094172321f;
+
+ float RMax = (float) getOverflowThreshold();
+ float Pi_2 = 2.0f * satans(1);
+
+ float _inReal = creals(z);
+ float _inImg = cimags(z);
+ float _outReal = 0;
+ float _outImg = 0;
+
+ /* Temporary variables */
+ float R2 = 0;
+ float S = 0;
+
+
+ if(_inImg == 0)
+ {
+ _outReal = satans(_inReal);
+ _outImg = 0;
+ }
+ else
+ {
+ R2 = _inReal * _inReal + _inImg * _inImg; /* Oo */
+ if(R2 > RMax)
+ {
+ if( dabss(_inImg) > RMax)
+ S = 0;
+ else
+ S = 1.0f / (((0.5f * _inReal) / _inImg) * _inReal + 0.5f * _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 = slnp1m1s(S) * 0.25f;
+ }
+ else
+ {
+ if(sabss(S) == 1 && sabss(_inReal) <= sAlim)
+ {
+ /* |s| >= SLIM => |z| is inside D+ or D- */
+ _outImg = _sign(0.5f,_inImg) * ( sLn2 - logf(sabss(_inReal)));
+ }
+ else
+ {
+ _outImg = 0.25f * logf((powf(_inReal,2) + powf((_inImg + 1.0f),2)) / (powf(_inReal,2) + powf((_inImg - 1.0f),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(sabss(1 - R2) + sabss(_inReal) <= sTol)
+ {/* |b| is very near 1 (and a is near 0) some cancellation occur in the (next) generic formula */
+ _outReal = 0.5f * atan2f(2.0f * _inReal, (1.0f - _inImg) * (1.0f + _inImg) - powf(_inReal,2.0f));
+ }
+ else
+ _outReal = 0.5f * atan2f(2.0f * _inReal, 1.0f - R2);
+ }
+
+ return FloatComplex(_outReal, _outImg);
+}
|