summaryrefslogtreecommitdiff
path: root/src/c/elementaryFunctions/atan/zatans.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/c/elementaryFunctions/atan/zatans.c')
-rw-r--r--src/c/elementaryFunctions/atan/zatans.c483
1 files changed, 242 insertions, 241 deletions
diff --git a/src/c/elementaryFunctions/atan/zatans.c b/src/c/elementaryFunctions/atan/zatans.c
index c511d790..4b8e9640 100644
--- a/src/c/elementaryFunctions/atan/zatans.c
+++ b/src/c/elementaryFunctions/atan/zatans.c
@@ -1,241 +1,242 @@
-/*
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2006-2008 - INRIA - Bruno JOFRET
- *
- * 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);
-}
+/*
+ * 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);
+}
generated by cgit (git 2.34.1) at 2025-03-13 12:25:40 +0000