convert to unix format

1 files changed, 242 insertions, 242 deletions

diff --git a/src/c/elementaryFunctions/atan/zatans.c b/src/c/elementaryFunctions/atan/zatans.c
index 4b8e9640..ab1354f7 100644
--- a/src/c/elementaryFunctions/atan/zatans.c
+++ b/src/c/elementaryFunctions/atan/zatans.c
@@ -1,242 +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 <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);
+}