From db464f35f5a10b58d9ed1085e0b462689adee583 Mon Sep 17 00:00:00 2001 From: Siddhesh Wani Date: Mon, 25 May 2015 14:46:31 +0530 Subject: Original Version --- src/c/elementaryFunctions/atan/zatans.c | 242 ++++++++++++++++++++++++++++++++ 1 file changed, 242 insertions(+) create mode 100644 src/c/elementaryFunctions/atan/zatans.c (limited to 'src/c/elementaryFunctions/atan/zatans.c') diff --git a/src/c/elementaryFunctions/atan/zatans.c b/src/c/elementaryFunctions/atan/zatans.c new file mode 100644 index 0000000..4b8e964 --- /dev/null +++ b/src/c/elementaryFunctions/atan/zatans.c @@ -0,0 +1,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 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 + + + 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 +#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); +} -- cgit