Fixed missing copyright notice.

15 files changed, 1536 insertions, 1518 deletions

diff --git a/src/c/elementaryFunctions/acos/cacoss.c b/src/c/elementaryFunctions/acos/cacoss.c
index 1059c8c3..97420313 100644
--- a/src/c/elementaryFunctions/acos/cacoss.c
+++ b/src/c/elementaryFunctions/acos/cacoss.c
@@ -1,146 +1,147 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2007-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
- *
- */
-
-/*
- * This fonction is a translation of fortran wacos write by Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>
- *     REFERENCE
- *        This is a Fortran-77 translation of an algorithm by
- *        T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which
- *        appears in their article :
- *          "Implementing the Complex Arcsine and Arccosine
- *           Functions Using Exception Handling", ACM, TOMS,
- *           Vol 23, No. 3, Sept 1997, p. 299-335
- */
-
-#include "acos.h"
-#include "atan.h"
-#include "log.h"
-#include "log1p.h"
-#include "sqrt.h"
-#include "abs.h"
-#include "lapack.h"
-#include "min.h"
-#include "max.h"
-
-#define localSign(x) (x>0 ? 1.0f : -1.0f)
-
-floatComplex		cacoss(floatComplex z) {
-	static float sfltPi		= 3.1415926535897932384626433f;
-	static float sfltPi_2		= 1.5707963267948966192313216f;
-	static float sfltLn2		= 0.6931471805599453094172321f;
-	static float sfltAcross	= 1.5f;
-	static float sfltBcross	= 0.6417f;
-
-	float fltLsup = ssqrts((float) getOverflowThreshold())/8.0f;
-	float fltLinf = 4.0f * ssqrts((float) getUnderflowThreshold());
-	float fltEpsm = ssqrts((float) getRelativeMachinePrecision());
-
-	float fltAbsReal	= sabss(creals(z));
-	float fltAbsImg		= sabss(cimags(z));
-	float fltSignReal	= localSign(creals(z));
-	float fltSignImg	= localSign(cimags(z));
-
-	float fltR = 0, fltS = 0, fltA = 0, fltB = 0;
-
-	float fltTemp = 0;
-
-	float _pfltReal = 0;
-	float _pfltImg = 0;
-
-	if( min(fltAbsReal, fltAbsImg) > fltLinf && max(fltAbsReal, fltAbsImg) <= fltLsup)
-	  {/* we are in the safe region */
-		fltR = ssqrts( (fltAbsReal + 1 )*(fltAbsReal + 1 ) + fltAbsImg*fltAbsImg);
-		fltS = ssqrts( (fltAbsReal - 1 )*(fltAbsReal - 1 ) + fltAbsImg*fltAbsImg);
-		fltA = 0.5f * ( fltR + fltS );
-		fltB = fltAbsReal / fltA;
-
-
-		/* compute the real part */
-		if(fltB <= sfltBcross)
-			_pfltReal = sacoss(fltB);
-		else if( fltAbsReal <= 1)
-			_pfltReal = satans(ssqrts(0.5f * (fltA + fltAbsReal) * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + (fltS + (1 - fltAbsReal)))) / fltAbsReal);
-		else
-			_pfltReal = satans((fltAbsImg * ssqrts(0.5f * ((fltA + fltAbsReal) / (fltR + (fltAbsReal + 1)) + (fltA + fltAbsReal) / (fltS + (fltAbsReal - 1))))) / fltAbsReal);
-
-		/* compute the imaginary part */
-		if(fltA <= sfltAcross)
-		{
-			float fltImg1 = 0;
-
-			if(fltAbsReal < 1)
-			  /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */
-				fltImg1 = 0.5f * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + fltAbsImg*fltAbsImg / (fltS + (1 - fltAbsReal)));
-			else
-			  /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */
-				fltImg1 = 0.5f * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + (fltS + (fltAbsReal - 1)));
-			/* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */
-			fltTemp = fltImg1 + ssqrts(fltImg1 *( fltA + 1));
-			_pfltImg = slog1ps(fltTemp);
-		}
-		else
-		  /* ai = log(A + sqrt(A**2 - 1.d0)) */
-			_pfltImg = slogs(fltA + ssqrts(fltA*fltA - 1));
-	}
-	else
-	  {/* evaluation in the special regions ... */
-		if(fltAbsImg <= fltEpsm * sabss(fltAbsReal - 1))
-		{
-			if(fltAbsReal < 1)
-			{
-				_pfltReal	= sacoss(fltAbsReal);
-				_pfltImg	= fltAbsImg / ssqrts((1 + fltAbsReal) * (1 - fltAbsReal));
-			}
-			else
-			{
-				_pfltReal = 0;
-				if(fltAbsReal <= fltLsup)
-				{
-					fltTemp		= (fltAbsReal - 1) + ssqrts((fltAbsReal - 1) * (fltAbsReal + 1));
-					_pfltImg	= slog1ps(fltTemp);
-				}
-				else
-					_pfltImg	= sfltLn2 + slogs(fltAbsReal);
-			}
-		}
-		else if(fltAbsImg < fltLinf)
-		{
-			_pfltReal	= ssqrts(fltAbsImg);
-			_pfltImg	= _pfltReal;
-		}
-		else if((fltEpsm * fltAbsImg - 1 >= fltAbsReal))
-		{
-			_pfltReal	= sfltPi_2;
-			_pfltImg	= sfltLn2 + slogs(fltAbsImg);
-		}
-		else if(fltAbsReal > 1)
-		{
-			_pfltReal	= satans(fltAbsImg / fltAbsReal);
-			fltTemp		= (fltAbsReal / fltAbsImg)*(fltAbsReal / fltAbsImg);
-			_pfltImg	= sfltLn2 + slogs(fltAbsImg) + 0.5f * slog1ps(fltTemp);
-		}
-		else
-		{
-			float fltTemp2 = ssqrts(1 + fltAbsImg*fltAbsImg);
-			_pfltReal	= sfltPi_2;
-			fltTemp		= 2 * fltAbsImg * (fltAbsImg + fltTemp2);
-			_pfltImg	= 0.5f * slog1ps(fltTemp);
-		}
-	}
-	if(fltSignReal < 0)
-		_pfltReal = sfltPi - _pfltReal;
-
-	if(fltAbsImg != 0 || fltSignReal < 0)
-		_pfltImg = - fltSignImg * _pfltImg;
-
-	return FloatComplex(_pfltReal, _pfltImg);
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2007-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

+ *

+ */

+

+/*

+ * This fonction is a translation of fortran wacos write by Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>

+ *     REFERENCE

+ *        This is a Fortran-77 translation of an algorithm by

+ *        T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which

+ *        appears in their article :

+ *          "Implementing the Complex Arcsine and Arccosine

+ *           Functions Using Exception Handling", ACM, TOMS,

+ *           Vol 23, No. 3, Sept 1997, p. 299-335

+ */

+

+#include "acos.h"

+#include "atan.h"

+#include "log.h"

+#include "log1p.h"

+#include "sqrt.h"

+#include "abs.h"

+#include "lapack.h"

+#include "min.h"

+#include "max.h"

+

+#define localSign(x) (x>0 ? 1.0f : -1.0f)

+

+floatComplex		cacoss(floatComplex z) {

+	static float sfltPi		= 3.1415926535897932384626433f;

+	static float sfltPi_2		= 1.5707963267948966192313216f;

+	static float sfltLn2		= 0.6931471805599453094172321f;

+	static float sfltAcross	= 1.5f;

+	static float sfltBcross	= 0.6417f;

+

+	float fltLsup = ssqrts((float) getOverflowThreshold())/8.0f;

+	float fltLinf = 4.0f * ssqrts((float) getUnderflowThreshold());

+	float fltEpsm = ssqrts((float) getRelativeMachinePrecision());

+

+	float fltAbsReal	= sabss(creals(z));

+	float fltAbsImg		= sabss(cimags(z));

+	float fltSignReal	= localSign(creals(z));

+	float fltSignImg	= localSign(cimags(z));

+

+	float fltR = 0, fltS = 0, fltA = 0, fltB = 0;

+

+	float fltTemp = 0;

+

+	float _pfltReal = 0;

+	float _pfltImg = 0;

+

+	if( min(fltAbsReal, fltAbsImg) > fltLinf && max(fltAbsReal, fltAbsImg) <= fltLsup)

+	  {/* we are in the safe region */

+		fltR = ssqrts( (fltAbsReal + 1 )*(fltAbsReal + 1 ) + fltAbsImg*fltAbsImg);

+		fltS = ssqrts( (fltAbsReal - 1 )*(fltAbsReal - 1 ) + fltAbsImg*fltAbsImg);

+		fltA = 0.5f * ( fltR + fltS );

+		fltB = fltAbsReal / fltA;

+

+

+		/* compute the real part */

+		if(fltB <= sfltBcross)

+			_pfltReal = sacoss(fltB);

+		else if( fltAbsReal <= 1)

+			_pfltReal = satans(ssqrts(0.5f * (fltA + fltAbsReal) * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + (fltS + (1 - fltAbsReal)))) / fltAbsReal);

+		else

+			_pfltReal = satans((fltAbsImg * ssqrts(0.5f * ((fltA + fltAbsReal) / (fltR + (fltAbsReal + 1)) + (fltA + fltAbsReal) / (fltS + (fltAbsReal - 1))))) / fltAbsReal);

+

+		/* compute the imaginary part */

+		if(fltA <= sfltAcross)

+		{

+			float fltImg1 = 0;

+

+			if(fltAbsReal < 1)

+			  /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */

+				fltImg1 = 0.5f * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + fltAbsImg*fltAbsImg / (fltS + (1 - fltAbsReal)));

+			else

+			  /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */

+				fltImg1 = 0.5f * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + (fltS + (fltAbsReal - 1)));

+			/* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */

+			fltTemp = fltImg1 + ssqrts(fltImg1 *( fltA + 1));

+			_pfltImg = slog1ps(fltTemp);

+		}

+		else

+		  /* ai = log(A + sqrt(A**2 - 1.d0)) */

+			_pfltImg = slogs(fltA + ssqrts(fltA*fltA - 1));

+	}

+	else

+	  {/* evaluation in the special regions ... */

+		if(fltAbsImg <= fltEpsm * sabss(fltAbsReal - 1))

+		{

+			if(fltAbsReal < 1)

+			{

+				_pfltReal	= sacoss(fltAbsReal);

+				_pfltImg	= fltAbsImg / ssqrts((1 + fltAbsReal) * (1 - fltAbsReal));

+			}

+			else

+			{

+				_pfltReal = 0;

+				if(fltAbsReal <= fltLsup)

+				{

+					fltTemp		= (fltAbsReal - 1) + ssqrts((fltAbsReal - 1) * (fltAbsReal + 1));

+					_pfltImg	= slog1ps(fltTemp);

+				}

+				else

+					_pfltImg	= sfltLn2 + slogs(fltAbsReal);

+			}

+		}

+		else if(fltAbsImg < fltLinf)

+		{

+			_pfltReal	= ssqrts(fltAbsImg);

+			_pfltImg	= _pfltReal;

+		}

+		else if((fltEpsm * fltAbsImg - 1 >= fltAbsReal))

+		{

+			_pfltReal	= sfltPi_2;

+			_pfltImg	= sfltLn2 + slogs(fltAbsImg);

+		}

+		else if(fltAbsReal > 1)

+		{

+			_pfltReal	= satans(fltAbsImg / fltAbsReal);

+			fltTemp		= (fltAbsReal / fltAbsImg)*(fltAbsReal / fltAbsImg);

+			_pfltImg	= sfltLn2 + slogs(fltAbsImg) + 0.5f * slog1ps(fltTemp);

+		}

+		else

+		{

+			float fltTemp2 = ssqrts(1 + fltAbsImg*fltAbsImg);

+			_pfltReal	= sfltPi_2;

+			fltTemp		= 2 * fltAbsImg * (fltAbsImg + fltTemp2);

+			_pfltImg	= 0.5f * slog1ps(fltTemp);

+		}

+	}

+	if(fltSignReal < 0)

+		_pfltReal = sfltPi - _pfltReal;

+

+	if(fltAbsImg != 0 || fltSignReal < 0)

+		_pfltImg = - fltSignImg * _pfltImg;

+

+	return FloatComplex(_pfltReal, _pfltImg);

+}

diff --git a/src/c/elementaryFunctions/acos/zacoss.c b/src/c/elementaryFunctions/acos/zacoss.c
index 7758a932..de6f3fe9 100644
--- a/src/c/elementaryFunctions/acos/zacoss.c
+++ b/src/c/elementaryFunctions/acos/zacoss.c
@@ -1,146 +1,147 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2007-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
- *
- */
-
-/*
- * This fonction is a translation of fortran wacos write by Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>
- *     REFERENCE
- *        This is a Fortran-77 translation of an algorithm by
- *        T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which
- *        appears in their article :
- *          "Implementing the Complex Arcsine and Arccosine
- *           Functions Using Exception Handling", ACM, TOMS,
- *           Vol 23, No. 3, Sept 1997, p. 299-335
- */
-
-#include "acos.h"
-#include "atan.h"
-#include "log.h"
-#include "log1p.h"
-#include "sqrt.h"
-#include "abs.h"
-#include "lapack.h"
-#include "min.h"
-#include "max.h"
-
-#define localSign(x) (x>0 ? 1 : -1)
-
-doubleComplex		zacoss(doubleComplex z) {
-	static double sdblPi		= 3.1415926535897932384626433;
-	static double sdblPi_2		= 1.5707963267948966192313216;
-	static double sdblLn2		= 0.6931471805599453094172321;
-	static double sdblAcross	= 1.5;
-	static double sdblBcross	= 0.6417;
-
-	double dblLsup = dsqrts(getOverflowThreshold())/8.0;
-	double dblLinf = 4.0 * dsqrts(getUnderflowThreshold());
-	double dblEpsm = dsqrts(getRelativeMachinePrecision());
-
-	double dblAbsReal	= dabss(zreals(z));
-	double dblAbsImg	= dabss(zimags(z));
-	double dblSignReal	= localSign(zreals(z));
-	double dblSignImg	= localSign(zimags(z));
-
-	double dblR = 0, dblS = 0, dblA = 0, dblB = 0;
-
-	double dblTemp = 0;
-
-	double _pdblReal = 0;
-	double _pdblImg = 0;
-
-	if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup)
-	  {/* we are in the safe region */
-		dblR = dsqrts( (dblAbsReal + 1 )*(dblAbsReal + 1 ) + dblAbsImg*dblAbsImg);
-		dblS = dsqrts( (dblAbsReal - 1 )*(dblAbsReal - 1 ) + dblAbsImg*dblAbsImg);
-		dblA = 0.5 * ( dblR + dblS );
-		dblB = dblAbsReal / dblA;
-
-
-		/* compute the real part */
-		if(dblB <= sdblBcross)
-			_pdblReal = dacoss(dblB);
-		else if( dblAbsReal <= 1)
-			_pdblReal = datans(dsqrts(0.5 * (dblA + dblAbsReal) * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal)))) / dblAbsReal);
-		else
-			_pdblReal = datans((dblAbsImg * dsqrts(0.5 * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal - 1))))) / dblAbsReal);
-
-		/* compute the imaginary part */
-		if(dblA <= sdblAcross)
-		{
-			double dblImg1 = 0;
-
-			if(dblAbsReal < 1)
-			  /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */
-				dblImg1 = 0.5 * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg*dblAbsImg / (dblS + (1 - dblAbsReal)));
-			else
-			  /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */
-				dblImg1 = 0.5 * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1)));
-			/* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */
-			dblTemp = dblImg1 + dsqrts(dblImg1 *( dblA + 1));
-			_pdblImg = dlog1ps(dblTemp);
-		}
-		else
-		  /* ai = log(A + sqrt(A**2 - 1.d0)) */
-			_pdblImg = dlogs(dblA + dsqrts(dblA*dblA - 1));
-	}
-	else
-	  {/* evaluation in the special regions ... */
-		if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1))
-		{
-			if(dblAbsReal < 1)
-			{
-				_pdblReal	= dacoss(dblAbsReal);
-				_pdblImg	= dblAbsImg / dsqrts((1 + dblAbsReal) * (1 - dblAbsReal));
-			}
-			else
-			{
-				_pdblReal = 0;
-				if(dblAbsReal <= dblLsup)
-				{
-					dblTemp		= (dblAbsReal - 1) + dsqrts((dblAbsReal - 1) * (dblAbsReal + 1));
-					_pdblImg	= dlog1ps(dblTemp);
-				}
-				else
-					_pdblImg	= sdblLn2 + dlogs(dblAbsReal);
-			}
-		}
-		else if(dblAbsImg < dblLinf)
-		{
-			_pdblReal	= dsqrts(dblAbsImg);
-			_pdblImg	= _pdblReal;
-		}
-		else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal))
-		{
-			_pdblReal	= sdblPi_2;
-			_pdblImg	= sdblLn2 + dlogs(dblAbsImg);
-		}
-		else if(dblAbsReal > 1)
-		{
-			_pdblReal	= datans(dblAbsImg / dblAbsReal);
-			dblTemp		= (dblAbsReal / dblAbsImg)*(dblAbsReal / dblAbsImg);
-			_pdblImg	= sdblLn2 + dlogs(dblAbsImg) + 0.5 * dlog1ps(dblTemp);
-		}
-		else
-		{
-			double dblTemp2 = dsqrts(1 + dblAbsImg*dblAbsImg);
-			_pdblReal	= sdblPi_2;
-			dblTemp		= 2 * dblAbsImg * (dblAbsImg + dblTemp2);
-			_pdblImg	= 0.5 * dlog1ps(dblTemp);
-		}
-	}
-	if(dblSignReal < 0)
-		_pdblReal = sdblPi - _pdblReal;
-
-	if(dblAbsImg != 0 || dblSignReal < 0)
-		_pdblImg = - dblSignImg * _pdblImg;
-
-	return DoubleComplex(_pdblReal, _pdblImg);
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2007-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

+ *

+ */

+

+/*

+ * This fonction is a translation of fortran wacos write by Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>

+ *     REFERENCE

+ *        This is a Fortran-77 translation of an algorithm by

+ *        T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which

+ *        appears in their article :

+ *          "Implementing the Complex Arcsine and Arccosine

+ *           Functions Using Exception Handling", ACM, TOMS,

+ *           Vol 23, No. 3, Sept 1997, p. 299-335

+ */

+

+#include "acos.h"

+#include "atan.h"

+#include "log.h"

+#include "log1p.h"

+#include "sqrt.h"

+#include "abs.h"

+#include "lapack.h"

+#include "min.h"

+#include "max.h"

+

+#define localSign(x) (x>0 ? 1 : -1)

+

+doubleComplex		zacoss(doubleComplex z) {

+	static double sdblPi		= 3.1415926535897932384626433;

+	static double sdblPi_2		= 1.5707963267948966192313216;

+	static double sdblLn2		= 0.6931471805599453094172321;

+	static double sdblAcross	= 1.5;

+	static double sdblBcross	= 0.6417;

+

+	double dblLsup = dsqrts(getOverflowThreshold())/8.0;

+	double dblLinf = 4.0 * dsqrts(getUnderflowThreshold());

+	double dblEpsm = dsqrts(getRelativeMachinePrecision());

+

+	double dblAbsReal	= dabss(zreals(z));

+	double dblAbsImg	= dabss(zimags(z));

+	double dblSignReal	= localSign(zreals(z));

+	double dblSignImg	= localSign(zimags(z));

+

+	double dblR = 0, dblS = 0, dblA = 0, dblB = 0;

+

+	double dblTemp = 0;

+

+	double _pdblReal = 0;

+	double _pdblImg = 0;

+

+	if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup)

+	  {/* we are in the safe region */

+		dblR = dsqrts( (dblAbsReal + 1 )*(dblAbsReal + 1 ) + dblAbsImg*dblAbsImg);

+		dblS = dsqrts( (dblAbsReal - 1 )*(dblAbsReal - 1 ) + dblAbsImg*dblAbsImg);

+		dblA = 0.5 * ( dblR + dblS );

+		dblB = dblAbsReal / dblA;

+

+

+		/* compute the real part */

+		if(dblB <= sdblBcross)

+			_pdblReal = dacoss(dblB);

+		else if( dblAbsReal <= 1)

+			_pdblReal = datans(dsqrts(0.5 * (dblA + dblAbsReal) * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal)))) / dblAbsReal);

+		else

+			_pdblReal = datans((dblAbsImg * dsqrts(0.5 * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal - 1))))) / dblAbsReal);

+

+		/* compute the imaginary part */

+		if(dblA <= sdblAcross)

+		{

+			double dblImg1 = 0;

+

+			if(dblAbsReal < 1)

+			  /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */

+				dblImg1 = 0.5 * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg*dblAbsImg / (dblS + (1 - dblAbsReal)));

+			else

+			  /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */

+				dblImg1 = 0.5 * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1)));

+			/* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */

+			dblTemp = dblImg1 + dsqrts(dblImg1 *( dblA + 1));

+			_pdblImg = dlog1ps(dblTemp);

+		}

+		else

+		  /* ai = log(A + sqrt(A**2 - 1.d0)) */

+			_pdblImg = dlogs(dblA + dsqrts(dblA*dblA - 1));

+	}

+	else

+	  {/* evaluation in the special regions ... */

+		if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1))

+		{

+			if(dblAbsReal < 1)

+			{

+				_pdblReal	= dacoss(dblAbsReal);

+				_pdblImg	= dblAbsImg / dsqrts((1 + dblAbsReal) * (1 - dblAbsReal));

+			}

+			else

+			{

+				_pdblReal = 0;

+				if(dblAbsReal <= dblLsup)

+				{

+					dblTemp		= (dblAbsReal - 1) + dsqrts((dblAbsReal - 1) * (dblAbsReal + 1));

+					_pdblImg	= dlog1ps(dblTemp);

+				}

+				else

+					_pdblImg	= sdblLn2 + dlogs(dblAbsReal);

+			}

+		}

+		else if(dblAbsImg < dblLinf)

+		{

+			_pdblReal	= dsqrts(dblAbsImg);

+			_pdblImg	= _pdblReal;

+		}

+		else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal))

+		{

+			_pdblReal	= sdblPi_2;

+			_pdblImg	= sdblLn2 + dlogs(dblAbsImg);

+		}

+		else if(dblAbsReal > 1)

+		{

+			_pdblReal	= datans(dblAbsImg / dblAbsReal);

+			dblTemp		= (dblAbsReal / dblAbsImg)*(dblAbsReal / dblAbsImg);

+			_pdblImg	= sdblLn2 + dlogs(dblAbsImg) + 0.5 * dlog1ps(dblTemp);

+		}

+		else

+		{

+			double dblTemp2 = dsqrts(1 + dblAbsImg*dblAbsImg);

+			_pdblReal	= sdblPi_2;

+			dblTemp		= 2 * dblAbsImg * (dblAbsImg + dblTemp2);

+			_pdblImg	= 0.5 * dlog1ps(dblTemp);

+		}

+	}

+	if(dblSignReal < 0)

+		_pdblReal = sdblPi - _pdblReal;

+

+	if(dblAbsImg != 0 || dblSignReal < 0)

+		_pdblImg = - dblSignImg * _pdblImg;

+

+	return DoubleComplex(_pdblReal, _pdblImg);

+}

diff --git a/src/c/elementaryFunctions/asin/casins.c b/src/c/elementaryFunctions/asin/casins.c
index 9c15ce53..35a4a8d8 100644
--- a/src/c/elementaryFunctions/asin/casins.c
+++ b/src/c/elementaryFunctions/asin/casins.c
@@ -1,144 +1,146 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2008-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
- *
- */
-
-/*
- *     REFERENCE
- *        This is a Fortran-77 translation of an algorithm by
- *        T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which
- *        appears in their article :
- *          "Implementing the Complex Arcsine and Arccosine
- *           Functions Using Exception Handling", ACM, TOMS,
- *           Vol 23, No. 3, Sept 1997, p. 299-335
- */
-
-#include "lapack.h"
-#include "asin.h"
-#include "atan.h"
-#include "sqrt.h"
-#include "abs.h"
-#include "log.h"
-#include "log1p.h"
-#include "min.h"
-#include "max.h"
-
-floatComplex		casins(floatComplex z) {
-  static float sdblPi_2		= 1.5707963267948966192313216f;
-  static float sdblLn2		= 0.6931471805599453094172321f;
-  static float sdblAcross	= 1.5f;
-  static float sdblBcross	= 0.6417f;
-
-  float dblLsup = ssqrts((float) getOverflowThreshold())/ 8.0f;
-  float dblLinf = 4.0f * ssqrts((float) getUnderflowThreshold());
-  float dblEpsm = ssqrts((float) getRelativeMachinePrecision());
-
-  float _dblReal	= creals(z);
-  float _dblImg		= cimags(z);
-
-  float dblAbsReal	= sabss(_dblReal);
-  float dblAbsImg	= sabss(_dblImg);
-  float iSignReal	= _dblReal < 0 ? -1.0f : 1.0f;
-  float iSignImg	= _dblImg < 0 ? -1.0f : 1.0f;
-
-  float dblR = 0, dblS = 0, dblA = 0, dblB = 0;
-
-  float dblTemp = 0;
-
-  float _pdblReal = 0;
-  float _pdblImg = 0;
-
-  if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup)
-    {
-      /* we are in the safe region */
-      dblR = ssqrts( (dblAbsReal + 1) * (dblAbsReal + 1) + dblAbsImg * dblAbsImg);
-      dblS = ssqrts( (dblAbsReal - 1) * (dblAbsReal - 1) + dblAbsImg * dblAbsImg);
-      dblA = (float) 0.5 * ( dblR + dblS );
-      dblB = dblAbsReal / dblA;
-
-
-      /* compute the real part */
-      if(dblB <= sdblBcross)
-	_pdblReal = sasins(dblB);
-      else if( dblAbsReal <= 1)
-	_pdblReal = satans(dblAbsReal / ssqrts( 0.5f * (dblA + dblAbsReal) * ( (dblAbsImg * dblAbsImg) / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal)))));
-      else
-	_pdblReal = satans(dblAbsReal / (dblAbsImg * ssqrts( 0.5f * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal-1))))));
-
-      /* compute the imaginary part */
-      if(dblA <= sdblAcross)
-	{
-	  float dblImg1 = 0;
-
-	  if(dblAbsReal < 1)
-	    /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */
-	    dblImg1 = 0.5f * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg * dblAbsImg / (dblS + (1 - dblAbsReal)));
-	  else
-	    /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */
-	    dblImg1 = 0.5f * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1)));
-	  /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */
-	  dblTemp = dblImg1 + ssqrts(dblImg1 * (dblA + 1));
-	  _pdblImg = slog1ps(dblTemp);
-	}
-      else
-	/* ai = log(A + sqrt(A**2 - 1.d0)) */
-	_pdblImg = slogs(dblA + ssqrts(dblA * dblA - (float) 1.0));
-    }
-  else
-    {
-      /* evaluation in the special regions ... */
-      if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1))
-	{
-	  if(dblAbsReal < 1)
-	    {
-	      _pdblReal	= sasins(dblAbsReal);
-	      _pdblImg	= dblAbsImg / ssqrts((1 + dblAbsReal) * (1 - dblAbsReal));
-	    }
-	  else
-	    {
-	      _pdblReal = sdblPi_2;
-	      if(dblAbsReal <= dblLsup)
-		{
-		  dblTemp		= (dblAbsReal - 1) + ssqrts((dblAbsReal - 1) * (dblAbsReal + 1));
-		  _pdblImg	= slog1ps(dblTemp);
-		}
-	      else
-		_pdblImg	= sdblLn2 + slogs(dblAbsReal);
-	    }
-	}
-      else if(dblAbsImg < dblLinf)
-	{
-	  _pdblReal	= sdblPi_2 - ssqrts(dblAbsImg);
-	  _pdblImg	= ssqrts(dblAbsImg);
-	}
-      else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal))
-	{
-	  _pdblReal	= dblAbsReal * dblAbsImg;
-	  _pdblImg	= sdblLn2 + slogs(dblAbsReal);
-	}
-      else if(dblAbsReal > 1)
-	{
-	  _pdblReal	= satans(dblAbsReal / dblAbsImg);
-	  dblTemp	= (dblAbsReal / dblAbsImg) * (dblAbsReal / dblAbsImg);
-	  _pdblImg	= sdblLn2 + slogs(dblAbsReal) + 0.5f * slog1ps(dblTemp);
-	}
-      else
-	{
-	  float dblTemp2 = ssqrts(1 + dblAbsImg * dblAbsImg);
-	  _pdblReal	= dblAbsReal / dblTemp2;
-	  dblTemp	= 2.0f * dblAbsImg * (dblAbsImg + dblTemp2);
-	  _pdblImg	= 0.5f * slog1ps(dblTemp);
-	}
-    }
-  _pdblReal *= iSignReal;
-  _pdblImg *= iSignImg;
-
-  return (FloatComplex(_pdblReal, _pdblImg));
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2008-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

+ *

+ */

+

+/*

+ *     REFERENCE

+ *        This is a Fortran-77 translation of an algorithm by

+ *        T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which

+ *        appears in their article :

+ *          "Implementing the Complex Arcsine and Arccosine

+ *           Functions Using Exception Handling", ACM, TOMS,

+ *           Vol 23, No. 3, Sept 1997, p. 299-335

+ *     Thanks to Tom Fairgrieve

+ */

+

+#include "lapack.h"

+#include "asin.h"

+#include "atan.h"

+#include "sqrt.h"

+#include "abs.h"

+#include "log.h"

+#include "log1p.h"

+#include "min.h"

+#include "max.h"

+

+floatComplex		casins(floatComplex z) {

+  static float sdblPi_2		= 1.5707963267948966192313216f;

+  static float sdblLn2		= 0.6931471805599453094172321f;

+  static float sdblAcross	= 1.5f;

+  static float sdblBcross	= 0.6417f;

+

+  float dblLsup = ssqrts((float) getOverflowThreshold())/ 8.0f;

+  float dblLinf = 4.0f * ssqrts((float) getUnderflowThreshold());

+  float dblEpsm = ssqrts((float) getRelativeMachinePrecision());

+

+  float _dblReal	= creals(z);

+  float _dblImg		= cimags(z);

+

+  float dblAbsReal	= sabss(_dblReal);

+  float dblAbsImg	= sabss(_dblImg);

+  float iSignReal	= _dblReal < 0 ? -1.0f : 1.0f;

+  float iSignImg	= _dblImg < 0 ? -1.0f : 1.0f;

+

+  float dblR = 0, dblS = 0, dblA = 0, dblB = 0;

+

+  float dblTemp = 0;

+

+  float _pdblReal = 0;

+  float _pdblImg = 0;

+

+  if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup)

+    {

+      /* we are in the safe region */

+      dblR = ssqrts( (dblAbsReal + 1) * (dblAbsReal + 1) + dblAbsImg * dblAbsImg);

+      dblS = ssqrts( (dblAbsReal - 1) * (dblAbsReal - 1) + dblAbsImg * dblAbsImg);

+      dblA = (float) 0.5 * ( dblR + dblS );

+      dblB = dblAbsReal / dblA;

+

+

+      /* compute the real part */

+      if(dblB <= sdblBcross)

+	_pdblReal = sasins(dblB);

+      else if( dblAbsReal <= 1)

+	_pdblReal = satans(dblAbsReal / ssqrts( 0.5f * (dblA + dblAbsReal) * ( (dblAbsImg * dblAbsImg) / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal)))));

+      else

+	_pdblReal = satans(dblAbsReal / (dblAbsImg * ssqrts( 0.5f * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal-1))))));

+

+      /* compute the imaginary part */

+      if(dblA <= sdblAcross)

+	{

+	  float dblImg1 = 0;

+

+	  if(dblAbsReal < 1)

+	    /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */

+	    dblImg1 = 0.5f * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg * dblAbsImg / (dblS + (1 - dblAbsReal)));

+	  else

+	    /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */

+	    dblImg1 = 0.5f * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1)));

+	  /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */

+	  dblTemp = dblImg1 + ssqrts(dblImg1 * (dblA + 1));

+	  _pdblImg = slog1ps(dblTemp);

+	}

+      else

+	/* ai = log(A + sqrt(A**2 - 1.d0)) */

+	_pdblImg = slogs(dblA + ssqrts(dblA * dblA - (float) 1.0));

+    }

+  else

+    {

+      /* evaluation in the special regions ... */

+      if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1))

+	{

+	  if(dblAbsReal < 1)

+	    {

+	      _pdblReal	= sasins(dblAbsReal);

+	      _pdblImg	= dblAbsImg / ssqrts((1 + dblAbsReal) * (1 - dblAbsReal));

+	    }

+	  else

+	    {

+	      _pdblReal = sdblPi_2;

+	      if(dblAbsReal <= dblLsup)

+		{

+		  dblTemp		= (dblAbsReal - 1) + ssqrts((dblAbsReal - 1) * (dblAbsReal + 1));

+		  _pdblImg	= slog1ps(dblTemp);

+		}

+	      else

+		_pdblImg	= sdblLn2 + slogs(dblAbsReal);

+	    }

+	}

+      else if(dblAbsImg < dblLinf)

+	{

+	  _pdblReal	= sdblPi_2 - ssqrts(dblAbsImg);

+	  _pdblImg	= ssqrts(dblAbsImg);

+	}

+      else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal))

+	{

+	  _pdblReal	= dblAbsReal * dblAbsImg;

+	  _pdblImg	= sdblLn2 + slogs(dblAbsReal);

+	}

+      else if(dblAbsReal > 1)

+	{

+	  _pdblReal	= satans(dblAbsReal / dblAbsImg);

+	  dblTemp	= (dblAbsReal / dblAbsImg) * (dblAbsReal / dblAbsImg);

+	  _pdblImg	= sdblLn2 + slogs(dblAbsReal) + 0.5f * slog1ps(dblTemp);

+	}

+      else

+	{

+	  float dblTemp2 = ssqrts(1 + dblAbsImg * dblAbsImg);

+	  _pdblReal	= dblAbsReal / dblTemp2;

+	  dblTemp	= 2.0f * dblAbsImg * (dblAbsImg + dblTemp2);

+	  _pdblImg	= 0.5f * slog1ps(dblTemp);

+	}

+    }

+  _pdblReal *= iSignReal;

+  _pdblImg *= iSignImg;

+

+  return (FloatComplex(_pdblReal, _pdblImg));

+}

diff --git a/src/c/elementaryFunctions/asin/zasins.c b/src/c/elementaryFunctions/asin/zasins.c
index 5d2110ff..5bd586a8 100644
--- a/src/c/elementaryFunctions/asin/zasins.c
+++ b/src/c/elementaryFunctions/asin/zasins.c
@@ -1,144 +1,146 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2007-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
- *
- */
-
-/*
- *     REFERENCE
- *        This is a Fortran-77 translation of an algorithm by
- *        T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which
- *        appears in their article :
- *          "Implementing the Complex Arcsine and Arccosine
- *           Functions Using Exception Handling", ACM, TOMS,
- *           Vol 23, No. 3, Sept 1997, p. 299-335
- */
-
-#include "lapack.h"
-#include "asin.h"
-#include "atan.h"
-#include "sqrt.h"
-#include "abs.h"
-#include "log.h"
-#include "log1p.h"
-#include "min.h"
-#include "max.h"
-
-doubleComplex		zasins(doubleComplex z) {
-  static double sdblPi_2	= 1.5707963267948966192313216;
-  static double sdblLn2		= 0.6931471805599453094172321;
-  static double sdblAcross	= 1.5;
-  static double sdblBcross	= 0.6417;
-
-  double dblLsup = dsqrts(getOverflowThreshold())/8.0;
-  double dblLinf = 4 * dsqrts(getUnderflowThreshold());
-  double dblEpsm = dsqrts(getRelativeMachinePrecision());
-
-  double _dblReal	= zreals(z);
-  double _dblImg	= zimags(z);
-
-  double dblAbsReal	= dabss(_dblReal);
-  double dblAbsImg	= dabss(_dblImg);
-  int iSignReal		= _dblReal < 0 ? -1 : 1;
-  int iSignImg		= _dblImg < 0 ? -1 : 1;
-
-  double dblR = 0, dblS = 0, dblA = 0, dblB = 0;
-
-  double dblTemp = 0;
-
-  double _pdblReal = 0;
-  double _pdblImg = 0;
-
-  if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup)
-    {
-      /* we are in the safe region */
-      dblR = dsqrts( (dblAbsReal + 1) * (dblAbsReal + 1) + dblAbsImg * dblAbsImg);
-      dblS = dsqrts( (dblAbsReal - 1) * (dblAbsReal - 1) + dblAbsImg * dblAbsImg);
-      dblA = 0.5 * ( dblR + dblS );
-      dblB = dblAbsReal / dblA;
-
-
-      /* compute the real part */
-      if(dblB <= sdblBcross)
-	_pdblReal = dasins(dblB);
-      else if( dblAbsReal <= 1)
-	_pdblReal = datans(dblAbsReal / dsqrts( 0.5 * (dblA + dblAbsReal) * ( (dblAbsImg * dblAbsImg) / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal)))));
-      else
-	_pdblReal = datans(dblAbsReal / (dblAbsImg * dsqrts(0.5 * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal-1))))));
-
-      /* compute the imaginary part */
-      if(dblA <= sdblAcross)
-	{
-	  double dblImg1 = 0;
-
-	  if(dblAbsReal < 1)
-	    /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */
-	    dblImg1 = 0.5 * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg * dblAbsImg / (dblS + (1 - dblAbsReal)));
-	  else
-	    /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */
-	    dblImg1 = 0.5 * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1)));
-	  /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */
-	  dblTemp = dblImg1 + dsqrts(dblImg1 * (dblA + 1));
-	  _pdblImg = dlog1ps(dblTemp);
-	}
-      else
-	/* ai = log(A + sqrt(A**2 - 1.d0)) */
-	_pdblImg = dlogs(dblA + dsqrts(dblA * dblA - 1));
-    }
-  else
-    {
-      /* evaluation in the special regions ... */
-      if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1))
-	{
-	  if(dblAbsReal < 1)
-	    {
-	      _pdblReal	= dasins(dblAbsReal);
-	      _pdblImg	= dblAbsImg / dsqrts((1 + dblAbsReal) * (1 - dblAbsReal));
-	    }
-	  else
-	    {
-	      _pdblReal = sdblPi_2;
-	      if(dblAbsReal <= dblLsup)
-		{
-		  dblTemp		= (dblAbsReal - 1) + dsqrts((dblAbsReal - 1) * (dblAbsReal + 1));
-		  _pdblImg	= dlog1ps(dblTemp);
-		}
-	      else
-		_pdblImg	= sdblLn2 + dlogs(dblAbsReal);
-	    }
-	}
-      else if(dblAbsImg < dblLinf)
-	{
-	  _pdblReal	= sdblPi_2 - dsqrts(dblAbsImg);
-	  _pdblImg	= dsqrts(dblAbsImg);
-	}
-      else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal))
-	{
-	  _pdblReal	= dblAbsReal * dblAbsImg;
-	  _pdblImg	= sdblLn2 + dlogs(dblAbsReal);
-	}
-      else if(dblAbsReal > 1)
-	{
-	  _pdblReal	= datans(dblAbsReal / dblAbsImg);
-	  dblTemp		= (dblAbsReal / dblAbsImg) * (dblAbsReal / dblAbsImg);
-	  _pdblImg	= sdblLn2 + dlogs(dblAbsReal) + 0.5 * dlog1ps(dblTemp);
-	}
-      else
-	{
-	  double dblTemp2 = dsqrts(1 + dblAbsImg * dblAbsImg);
-	  _pdblReal	= dblAbsReal / dblTemp2;
-	  dblTemp		= 2 * dblAbsImg * (dblAbsImg + dblTemp2);
-	  _pdblImg	= 0.5 * dlog1ps(dblTemp);
-	}
-    }
-  _pdblReal *= iSignReal;
-  _pdblImg *= iSignImg;
-
-  return (DoubleComplex(_pdblReal, _pdblImg));
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2007-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

+ *

+ */

+

+/*

+ *     REFERENCE

+ *        This is a Fortran-77 translation of an algorithm by

+ *        T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which

+ *        appears in their article :

+ *          "Implementing the Complex Arcsine and Arccosine

+ *           Functions Using Exception Handling", ACM, TOMS,

+ *           Vol 23, No. 3, Sept 1997, p. 299-335

+ *     Thanks to Tom Fairgrieve

+ */

+

+#include "lapack.h"

+#include "asin.h"

+#include "atan.h"

+#include "sqrt.h"

+#include "abs.h"

+#include "log.h"

+#include "log1p.h"

+#include "min.h"

+#include "max.h"

+

+doubleComplex		zasins(doubleComplex z) {

+  static double sdblPi_2	= 1.5707963267948966192313216;

+  static double sdblLn2		= 0.6931471805599453094172321;

+  static double sdblAcross	= 1.5;

+  static double sdblBcross	= 0.6417;

+

+  double dblLsup = dsqrts(getOverflowThreshold())/8.0;

+  double dblLinf = 4 * dsqrts(getUnderflowThreshold());

+  double dblEpsm = dsqrts(getRelativeMachinePrecision());

+

+  double _dblReal	= zreals(z);

+  double _dblImg	= zimags(z);

+

+  double dblAbsReal	= dabss(_dblReal);

+  double dblAbsImg	= dabss(_dblImg);

+  int iSignReal		= _dblReal < 0 ? -1 : 1;

+  int iSignImg		= _dblImg < 0 ? -1 : 1;

+

+  double dblR = 0, dblS = 0, dblA = 0, dblB = 0;

+

+  double dblTemp = 0;

+

+  double _pdblReal = 0;

+  double _pdblImg = 0;

+

+  if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup)

+    {

+      /* we are in the safe region */

+      dblR = dsqrts( (dblAbsReal + 1) * (dblAbsReal + 1) + dblAbsImg * dblAbsImg);

+      dblS = dsqrts( (dblAbsReal - 1) * (dblAbsReal - 1) + dblAbsImg * dblAbsImg);

+      dblA = 0.5 * ( dblR + dblS );

+      dblB = dblAbsReal / dblA;

+

+

+      /* compute the real part */

+      if(dblB <= sdblBcross)

+	_pdblReal = dasins(dblB);

+      else if( dblAbsReal <= 1)

+	_pdblReal = datans(dblAbsReal / dsqrts( 0.5 * (dblA + dblAbsReal) * ( (dblAbsImg * dblAbsImg) / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal)))));

+      else

+	_pdblReal = datans(dblAbsReal / (dblAbsImg * dsqrts(0.5 * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal-1))))));

+

+      /* compute the imaginary part */

+      if(dblA <= sdblAcross)

+	{

+	  double dblImg1 = 0;

+

+	  if(dblAbsReal < 1)

+	    /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */

+	    dblImg1 = 0.5 * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg * dblAbsImg / (dblS + (1 - dblAbsReal)));

+	  else

+	    /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */

+	    dblImg1 = 0.5 * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1)));

+	  /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */

+	  dblTemp = dblImg1 + dsqrts(dblImg1 * (dblA + 1));

+	  _pdblImg = dlog1ps(dblTemp);

+	}

+      else

+	/* ai = log(A + sqrt(A**2 - 1.d0)) */

+	_pdblImg = dlogs(dblA + dsqrts(dblA * dblA - 1));

+    }

+  else

+    {

+      /* evaluation in the special regions ... */

+      if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1))

+	{

+	  if(dblAbsReal < 1)

+	    {

+	      _pdblReal	= dasins(dblAbsReal);

+	      _pdblImg	= dblAbsImg / dsqrts((1 + dblAbsReal) * (1 - dblAbsReal));

+	    }

+	  else

+	    {

+	      _pdblReal = sdblPi_2;

+	      if(dblAbsReal <= dblLsup)

+		{

+		  dblTemp		= (dblAbsReal - 1) + dsqrts((dblAbsReal - 1) * (dblAbsReal + 1));

+		  _pdblImg	= dlog1ps(dblTemp);

+		}

+	      else

+		_pdblImg	= sdblLn2 + dlogs(dblAbsReal);

+	    }

+	}

+      else if(dblAbsImg < dblLinf)

+	{

+	  _pdblReal	= sdblPi_2 - dsqrts(dblAbsImg);

+	  _pdblImg	= dsqrts(dblAbsImg);

+	}

+      else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal))

+	{

+	  _pdblReal	= dblAbsReal * dblAbsImg;

+	  _pdblImg	= sdblLn2 + dlogs(dblAbsReal);

+	}

+      else if(dblAbsReal > 1)

+	{

+	  _pdblReal	= datans(dblAbsReal / dblAbsImg);

+	  dblTemp		= (dblAbsReal / dblAbsImg) * (dblAbsReal / dblAbsImg);

+	  _pdblImg	= sdblLn2 + dlogs(dblAbsReal) + 0.5 * dlog1ps(dblTemp);

+	}

+      else

+	{

+	  double dblTemp2 = dsqrts(1 + dblAbsImg * dblAbsImg);

+	  _pdblReal	= dblAbsReal / dblTemp2;

+	  dblTemp		= 2 * dblAbsImg * (dblAbsImg + dblTemp2);

+	  _pdblImg	= 0.5 * dlog1ps(dblTemp);

+	}

+    }

+  _pdblReal *= iSignReal;

+  _pdblImg *= iSignImg;

+

+  return (DoubleComplex(_pdblReal, _pdblImg));

+}

diff --git a/src/c/elementaryFunctions/atan/catans.c b/src/c/elementaryFunctions/atan/catans.c
index d3a232fb..d2081181 100644
--- a/src/c/elementaryFunctions/atan/catans.c
+++ b/src/c/elementaryFunctions/atan/catans.c
@@ -1,248 +1,249 @@
-/*
- *  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 "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);
-}
+/*

+ *  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);

+}

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);

+}

diff --git a/src/c/elementaryFunctions/lnp1m1/dlnp1m1s.c b/src/c/elementaryFunctions/lnp1m1/dlnp1m1s.c
index 28cdbe0c..7e1759be 100644
--- a/src/c/elementaryFunctions/lnp1m1/dlnp1m1s.c
+++ b/src/c/elementaryFunctions/lnp1m1/dlnp1m1s.c
@@ -1,76 +1,77 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2008-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
- *
- */
-
-#include "lnp1m1.h"
-#include "abs.h"
-
-/*
-	PURPOSE :  Compute   v = log ( (1 + s)/(1 - s) )
-	for small s, this is for |s| < SLIM = 0.20
-
-	ALGORITHM :
-	1/ if |s| is "very small" we use a truncated
-	taylor dvp (by keeping 3 terms) from :
-	2        4          6
-	t = 2 * s * ( 1 + 1/3 s  + 1/5 s  + [ 1/7 s  + ....] )
-	2        4
-	t = 2 * s * ( 1 + 1/3 s  + 1/5 s  + er)
-
-	The limit E until we use this formula may be simply
-	gotten so that the negliged part er is such that :
-	2        4
-	(#) er <= epsm * ( 1 + 1/3 s  + 1/5 s )   for all |s|<= E
-
-	As  er  = 1/7 s^6 + 1/9 s^8 + ...
-	er <= 1/7 * s^6 ( 1 + s^2 + s^4 + ...) = 1/7  s^6/(1-s^2)
-
-	the inequality (#) is forced if :
-
-	1/7  s^6 / (1-s^2)  <= epsm * ( 1 + 1/3 s^2  + 1/5 s^4 )
-
-	s^6 <= 7 epsm * (1 - 2/3 s^2 - 3/15 s^4 - 1/5 s^6)
-
-	So that E is very near (7 epsm)^(1/6) (approximately 3.032d-3):
-
-	2/ For larger |s| we used a minimax polynome :
-
-	yi = s * (2  + d3 s^3 + d5 s^5 .... + d13 s^13 + d15 s^15)
-
-	This polynome was computed (by some remes algorithm) following
-	(*) the sin(x) example (p 39) of the book :
-
-	"ELEMENTARY FUNCTIONS"
-	"Algorithms and implementation"
-	J.M. Muller (Birkhauser)
-
-	(*) without the additionnal raffinement to get the first coefs
-	very near floating point numbers)
-*/
-double dlnp1m1s(double Var)
-{
-	static double D3	= 0.66666666666672679472;
-	static double D5	= 0.39999999996176889299;
-	static double D7	= 0.28571429392829380980;
-	static double D9	= 0.22222138684562683797;
-	static double D11	= 0.18186349187499222459;
-	static double D13	= 0.15250315884469364710;
-	static double D15	= 0.15367270224757008114;
-	static double E		= 3.032E-3;
-	static double C3	= 2.0/3.0;
-	static double C5	= 2.0/5.0;
-
-	double S2 = Var * Var;
-	if( dabss(Var) <= E)
-		return Var * (2 + S2 * (C3 + C5 * S2));
-	else
-		return Var * (2 + S2 * (D3 + S2 * (D5 + S2 * (D7 + S2 * (D9 + S2 * (D11 + S2 * (D13 + S2 * D15)))))));
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2008-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

+ *

+ */

+

+#include "lnp1m1.h"

+#include "abs.h"

+

+/*

+	PURPOSE :  Compute   v = log ( (1 + s)/(1 - s) )

+	for small s, this is for |s| < SLIM = 0.20

+

+	ALGORITHM :

+	1/ if |s| is "very small" we use a truncated

+	taylor dvp (by keeping 3 terms) from :

+	2        4          6

+	t = 2 * s * ( 1 + 1/3 s  + 1/5 s  + [ 1/7 s  + ....] )

+	2        4

+	t = 2 * s * ( 1 + 1/3 s  + 1/5 s  + er)

+

+	The limit E until we use this formula may be simply

+	gotten so that the negliged part er is such that :

+	2        4

+	(#) er <= epsm * ( 1 + 1/3 s  + 1/5 s )   for all |s|<= E

+

+	As  er  = 1/7 s^6 + 1/9 s^8 + ...

+	er <= 1/7 * s^6 ( 1 + s^2 + s^4 + ...) = 1/7  s^6/(1-s^2)

+

+	the inequality (#) is forced if :

+

+	1/7  s^6 / (1-s^2)  <= epsm * ( 1 + 1/3 s^2  + 1/5 s^4 )

+

+	s^6 <= 7 epsm * (1 - 2/3 s^2 - 3/15 s^4 - 1/5 s^6)

+

+	So that E is very near (7 epsm)^(1/6) (approximately 3.032d-3):

+

+	2/ For larger |s| we used a minimax polynome :

+

+	yi = s * (2  + d3 s^3 + d5 s^5 .... + d13 s^13 + d15 s^15)

+

+	This polynome was computed (by some remes algorithm) following

+	(*) the sin(x) example (p 39) of the book :

+

+	"ELEMENTARY FUNCTIONS"

+	"Algorithms and implementation"

+	J.M. Muller (Birkhauser)

+

+	(*) without the additionnal raffinement to get the first coefs

+	very near floating point numbers)

+*/

+double dlnp1m1s(double Var)

+{

+	static double D3	= 0.66666666666672679472;

+	static double D5	= 0.39999999996176889299;

+	static double D7	= 0.28571429392829380980;

+	static double D9	= 0.22222138684562683797;

+	static double D11	= 0.18186349187499222459;

+	static double D13	= 0.15250315884469364710;

+	static double D15	= 0.15367270224757008114;

+	static double E		= 3.032E-3;

+	static double C3	= 2.0/3.0;

+	static double C5	= 2.0/5.0;

+

+	double S2 = Var * Var;

+	if( dabss(Var) <= E)

+		return Var * (2 + S2 * (C3 + C5 * S2));

+	else

+		return Var * (2 + S2 * (D3 + S2 * (D5 + S2 * (D7 + S2 * (D9 + S2 * (D11 + S2 * (D13 + S2 * D15)))))));

+}

diff --git a/src/c/elementaryFunctions/lnp1m1/slnp1m1s.c b/src/c/elementaryFunctions/lnp1m1/slnp1m1s.c
index 6c991cc0..9940810c 100644
--- a/src/c/elementaryFunctions/lnp1m1/slnp1m1s.c
+++ b/src/c/elementaryFunctions/lnp1m1/slnp1m1s.c
@@ -1,76 +1,77 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2007-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
- *
- */
-
-#include "lnp1m1.h"
-#include "abs.h"
-
-/*
-  PURPOSE :  Compute   v = log ( (1 + s)/(1 - s) )
-  for small s, this is for |s| < SLIM = 0.20
-
-  ALGORITHM :
-  1/ if |s| is "very small" we use a truncated
-  taylor dvp (by keeping 3 terms) from :
-  2        4          6
-  t = 2 * s * ( 1 + 1/3 s  + 1/5 s  + [ 1/7 s  + ....] )
-  2        4
-  t = 2 * s * ( 1 + 1/3 s  + 1/5 s  + er)
-
-  The limit E until we use this formula may be simply
-  gotten so that the negliged part er is such that :
-  2        4
-  (#) er <= epsm * ( 1 + 1/3 s  + 1/5 s )   for all |s|<= E
-
-  As  er  = 1/7 s^6 + 1/9 s^8 + ...
-  er <= 1/7 * s^6 ( 1 + s^2 + s^4 + ...) = 1/7  s^6/(1-s^2)
-
-  the inequality (#) is forced if :
-
-  1/7  s^6 / (1-s^2)  <= epsm * ( 1 + 1/3 s^2  + 1/5 s^4 )
-
-  s^6 <= 7 epsm * (1 - 2/3 s^2 - 3/15 s^4 - 1/5 s^6)
-
-  So that E is very near (7 epsm)^(1/6) (approximately 3.032d-3):
-
-  2/ For larger |s| we used a minimax polynome :
-
-  yi = s * (2  + d3 s^3 + d5 s^5 .... + d13 s^13 + d15 s^15)
-
-  This polynome was computed (by some remes algorithm) following
-  (*) the sin(x) example (p 39) of the book :
-
-  "ELEMENTARY FUNCTIONS"
-  "Algorithms and implementation"
-  J.M. Muller (Birkhauser)
-
-  (*) without the additionnal raffinement to get the first coefs
-  very near floating point numbers)
-*/
-float slnp1m1s(float Var)
-{
-  static float D3	= 0.66666666666672679472f;
-  static float D5	= 0.39999999996176889299f;
-  static float D7	= 0.28571429392829380980f;
-  static float D9	= 0.22222138684562683797f;
-  static float D11	= 0.18186349187499222459f;
-  static float D13	= 0.15250315884469364710f;
-  static float D15	= 0.15367270224757008114f;
-  static float E	= 3.032E-3f;
-  static float C3	= 2.0f/3.0f;
-  static float C5	= 2.0f/5.0f;
-
-  float S2 = Var * Var;
-  if( sabss(Var) <= E)
-    return Var * (2 + S2 * (C3 + C5 * S2));
-  else
-    return Var * (2 + S2 * (D3 + S2 * (D5 + S2 * (D7 + S2 * (D9 + S2 * (D11 + S2 * (D13 + S2 * D15)))))));
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2007-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

+ *

+ */

+

+#include "lnp1m1.h"

+#include "abs.h"

+

+/*

+  PURPOSE :  Compute   v = log ( (1 + s)/(1 - s) )

+  for small s, this is for |s| < SLIM = 0.20

+

+  ALGORITHM :

+  1/ if |s| is "very small" we use a truncated

+  taylor dvp (by keeping 3 terms) from :

+  2        4          6

+  t = 2 * s * ( 1 + 1/3 s  + 1/5 s  + [ 1/7 s  + ....] )

+  2        4

+  t = 2 * s * ( 1 + 1/3 s  + 1/5 s  + er)

+

+  The limit E until we use this formula may be simply

+  gotten so that the negliged part er is such that :

+  2        4

+  (#) er <= epsm * ( 1 + 1/3 s  + 1/5 s )   for all |s|<= E

+

+  As  er  = 1/7 s^6 + 1/9 s^8 + ...

+  er <= 1/7 * s^6 ( 1 + s^2 + s^4 + ...) = 1/7  s^6/(1-s^2)

+

+  the inequality (#) is forced if :

+

+  1/7  s^6 / (1-s^2)  <= epsm * ( 1 + 1/3 s^2  + 1/5 s^4 )

+

+  s^6 <= 7 epsm * (1 - 2/3 s^2 - 3/15 s^4 - 1/5 s^6)

+

+  So that E is very near (7 epsm)^(1/6) (approximately 3.032d-3):

+

+  2/ For larger |s| we used a minimax polynome :

+

+  yi = s * (2  + d3 s^3 + d5 s^5 .... + d13 s^13 + d15 s^15)

+

+  This polynome was computed (by some remes algorithm) following

+  (*) the sin(x) example (p 39) of the book :

+

+  "ELEMENTARY FUNCTIONS"

+  "Algorithms and implementation"

+  J.M. Muller (Birkhauser)

+

+  (*) without the additionnal raffinement to get the first coefs

+  very near floating point numbers)

+*/

+float slnp1m1s(float Var)

+{

+  static float D3	= 0.66666666666672679472f;

+  static float D5	= 0.39999999996176889299f;

+  static float D7	= 0.28571429392829380980f;

+  static float D9	= 0.22222138684562683797f;

+  static float D11	= 0.18186349187499222459f;

+  static float D13	= 0.15250315884469364710f;

+  static float D15	= 0.15367270224757008114f;

+  static float E	= 3.032E-3f;

+  static float C3	= 2.0f/3.0f;

+  static float C5	= 2.0f/5.0f;

+

+  float S2 = Var * Var;

+  if( sabss(Var) <= E)

+    return Var * (2 + S2 * (C3 + C5 * S2));

+  else

+    return Var * (2 + S2 * (D3 + S2 * (D5 + S2 * (D7 + S2 * (D9 + S2 * (D11 + S2 * (D13 + S2 * D15)))))));

+}

diff --git a/src/c/elementaryFunctions/log/clogs.c b/src/c/elementaryFunctions/log/clogs.c
index d402d1d1..3b9f6911 100644
--- a/src/c/elementaryFunctions/log/clogs.c
+++ b/src/c/elementaryFunctions/log/clogs.c
@@ -1,64 +1,65 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2007-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
- *
- */
-
-#include <math.h>
-#include "log.h"
-#include "lapack.h"
-#include "log1p.h"
-#include "pythag.h"
-
-floatComplex	clogs(floatComplex in) {
-  static float sR2	= 1.41421356237309504f;
-
-  float _RealIn = creals(in);
-  float _ImgIn = cimags(in);
-
-  float _RealOut = 0;
-  float _ImgOut = 0;
-
-  float RMax = (float) getOverflowThreshold();
-  float LInf = sqrtf((float) getUnderflowThreshold());
-  float LSup = sqrtf(0.5f * RMax);
-
-  float AbsReal = fabsf(_RealIn);
-  float AbsImg	= fabsf(_ImgIn);
-
-  _ImgOut = atan2f(_ImgIn, _RealIn);
-
-  if(_ImgIn > _RealIn)
-    {/* switch Real part and Imaginary part */
-      float Temp	= AbsReal;
-      AbsReal		= AbsImg;
-      AbsImg		= Temp;
-    }
-
-  if((0.5 <= AbsReal) && (AbsReal <= sR2))
-    _RealOut = 0.5f * slog1ps((AbsReal - 1.0f) * (AbsReal + 1.0f) + AbsImg * AbsImg);
-  else if(LInf < AbsImg && AbsReal < LSup)
-    _RealOut = 0.5f * slogs(AbsReal * AbsReal + AbsImg * AbsImg);
-  else if(AbsReal > RMax)
-    _RealOut = AbsReal;
-  else
-    {
-      float Temp = spythags(AbsReal, AbsImg);
-      if(Temp <= RMax)
-	{
-	  _RealOut = slogs(Temp);
-	}
-      else /* handle rare spurious overflow with : */
-	{
-	  float Temp2 = AbsImg/AbsReal;
-	  _RealOut = slogs(AbsReal) + 0.5f * slog1ps(Temp2 * Temp2);
-	}
-    }
-  return FloatComplex(_RealOut, _ImgOut);
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2007-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

+ *

+ */

+

+#include <math.h>

+#include "log.h"

+#include "lapack.h"

+#include "log1p.h"

+#include "pythag.h"

+

+floatComplex	clogs(floatComplex in) {

+  static float sR2	= 1.41421356237309504f;

+

+  float _RealIn = creals(in);

+  float _ImgIn = cimags(in);

+

+  float _RealOut = 0;

+  float _ImgOut = 0;

+

+  float RMax = (float) getOverflowThreshold();

+  float LInf = sqrtf((float) getUnderflowThreshold());

+  float LSup = sqrtf(0.5f * RMax);

+

+  float AbsReal = fabsf(_RealIn);

+  float AbsImg	= fabsf(_ImgIn);

+

+  _ImgOut = atan2f(_ImgIn, _RealIn);

+

+  if(_ImgIn > _RealIn)

+    {/* switch Real part and Imaginary part */

+      float Temp	= AbsReal;

+      AbsReal		= AbsImg;

+      AbsImg		= Temp;

+    }

+

+  if((0.5 <= AbsReal) && (AbsReal <= sR2))

+    _RealOut = 0.5f * slog1ps((AbsReal - 1.0f) * (AbsReal + 1.0f) + AbsImg * AbsImg);

+  else if(LInf < AbsImg && AbsReal < LSup)

+    _RealOut = 0.5f * slogs(AbsReal * AbsReal + AbsImg * AbsImg);

+  else if(AbsReal > RMax)

+    _RealOut = AbsReal;

+  else

+    {

+      float Temp = spythags(AbsReal, AbsImg);

+      if(Temp <= RMax)

+	{

+	  _RealOut = slogs(Temp);

+	}

+      else /* handle rare spurious overflow with : */

+	{

+	  float Temp2 = AbsImg/AbsReal;

+	  _RealOut = slogs(AbsReal) + 0.5f * slog1ps(Temp2 * Temp2);

+	}

+    }

+  return FloatComplex(_RealOut, _ImgOut);

+}

diff --git a/src/c/elementaryFunctions/log/zlogs.c b/src/c/elementaryFunctions/log/zlogs.c
index dd0792a9..e5e9ded1 100644
--- a/src/c/elementaryFunctions/log/zlogs.c
+++ b/src/c/elementaryFunctions/log/zlogs.c
@@ -1,64 +1,65 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2007-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
- *
- */
-
-#include <math.h>
-#include "log.h"
-#include "lapack.h"
-#include "log1p.h"
-#include "pythag.h"
-
-doubleComplex	zlogs(doubleComplex in) {
-  static double sR2	= 1.41421356237309504;
-
-  double _RealIn = zreals(in);
-  double _ImgIn = zimags(in);
-
-  double _RealOut = 0;
-  double _ImgOut = 0;
-
-  double RMax = getOverflowThreshold();
-  double LInf = sqrt(getUnderflowThreshold());
-  double LSup = sqrt(0.5 * RMax);
-
-  double AbsReal = fabs(_RealIn);
-  double AbsImg	= fabs(_ImgIn);
-
-  _ImgOut = atan2(_ImgIn, _RealIn);
-
-  if(_ImgIn > _RealIn)
-    {/* switch Real part and Imaginary part */
-      double Temp	= AbsReal;
-      AbsReal		= AbsImg;
-      AbsImg		= Temp;
-    }
-
-  if((0.5 <= AbsReal) && (AbsReal <= sR2))
-    _RealOut = 0.5 * dlog1ps((AbsReal - 1) * (AbsReal + 1) + AbsImg * AbsImg);
-  else if(LInf < AbsImg && AbsReal < LSup)
-    _RealOut = 0.5 * dlogs(AbsReal * AbsReal + AbsImg * AbsImg);
-  else if(AbsReal > RMax)
-    _RealOut = AbsReal;
-  else
-    {
-      double Temp = dpythags(AbsReal, AbsImg);
-      if(Temp <= RMax)
-	{
-	  _RealOut = dlogs(Temp);
-	}
-      else /* handle rare spurious overflow with : */
-	{
-	  double Temp2 = AbsImg/AbsReal;
-	  _RealOut = dlogs(AbsReal) + 0.5 * dlog1ps(Temp2 * Temp2);
-	}
-    }
-  return DoubleComplex(_RealOut, _ImgOut);
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2007-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

+ *

+ */

+

+#include <math.h>

+#include "log.h"

+#include "lapack.h"

+#include "log1p.h"

+#include "pythag.h"

+

+doubleComplex	zlogs(doubleComplex in) {

+  static double sR2	= 1.41421356237309504;

+

+  double _RealIn = zreals(in);

+  double _ImgIn = zimags(in);

+

+  double _RealOut = 0;

+  double _ImgOut = 0;

+

+  double RMax = getOverflowThreshold();

+  double LInf = sqrt(getUnderflowThreshold());

+  double LSup = sqrt(0.5 * RMax);

+

+  double AbsReal = fabs(_RealIn);

+  double AbsImg	= fabs(_ImgIn);

+

+  _ImgOut = atan2(_ImgIn, _RealIn);

+

+  if(_ImgIn > _RealIn)

+    {/* switch Real part and Imaginary part */

+      double Temp	= AbsReal;

+      AbsReal		= AbsImg;

+      AbsImg		= Temp;

+    }

+

+  if((0.5 <= AbsReal) && (AbsReal <= sR2))

+    _RealOut = 0.5 * dlog1ps((AbsReal - 1) * (AbsReal + 1) + AbsImg * AbsImg);

+  else if(LInf < AbsImg && AbsReal < LSup)

+    _RealOut = 0.5 * dlogs(AbsReal * AbsReal + AbsImg * AbsImg);

+  else if(AbsReal > RMax)

+    _RealOut = AbsReal;

+  else

+    {

+      double Temp = dpythags(AbsReal, AbsImg);

+      if(Temp <= RMax)

+	{

+	  _RealOut = dlogs(Temp);

+	}

+      else /* handle rare spurious overflow with : */

+	{

+	  double Temp2 = AbsImg/AbsReal;

+	  _RealOut = dlogs(AbsReal) + 0.5 * dlog1ps(Temp2 * Temp2);

+	}

+    }

+  return DoubleComplex(_RealOut, _ImgOut);

+}

diff --git a/src/c/elementaryFunctions/log1p/dlog1ps.c b/src/c/elementaryFunctions/log1p/dlog1ps.c
index 4c055c86..e75a05fd 100644
--- a/src/c/elementaryFunctions/log1p/dlog1ps.c
+++ b/src/c/elementaryFunctions/log1p/dlog1ps.c
@@ -1,33 +1,34 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2008-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
- *
- */
-
-#include "log1p.h"
-#include "log.h"
-#include "lnp1m1.h"
-
-double		dlog1ps(double in) {
-  static double A = -1.0/3.0;
-  static double B = 0.5;
-
-  if(in < -1)
-    {/* got NaN */
-      return (in - in) / (in - in); /* NaN */
-    }
-  else if(A <= in && in <= B)
-    {/* use the function log((1+g)/(1-g)) with g = x/(x + 2) */
-      return dlnp1m1s(in / (in + 2));
-    }
-  else
-    {/* use the standard formula */
-      return dlogs(in + 1);
-    }
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2008-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

+ *

+ */

+

+#include "log1p.h"

+#include "log.h"

+#include "lnp1m1.h"

+

+double		dlog1ps(double in) {

+  static double A = -1.0/3.0;

+  static double B = 0.5;

+

+  if(in < -1)

+    {/* got NaN */

+      return (in - in) / (in - in); /* NaN */

+    }

+  else if(A <= in && in <= B)

+    {/* use the function log((1+g)/(1-g)) with g = x/(x + 2) */

+      return dlnp1m1s(in / (in + 2));

+    }

+  else

+    {/* use the standard formula */

+      return dlogs(in + 1);

+    }

+}

diff --git a/src/c/elementaryFunctions/log1p/slog1ps.c b/src/c/elementaryFunctions/log1p/slog1ps.c
index 501de6b0..04786524 100644
--- a/src/c/elementaryFunctions/log1p/slog1ps.c
+++ b/src/c/elementaryFunctions/log1p/slog1ps.c
@@ -1,33 +1,35 @@
-/*
- *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- *  Copyright (C) 2008-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
- *
- */
-
-#include "log1p.h"
-#include "log.h"
-#include "lnp1m1.h"
-
-float		slog1ps(float in) {
-  static double A = -1.0/3.0;
-  static double B = 0.5;
-
-  if(in < -1)
-    {/* got NaN */
-      return (in - in) / (in - in); /* NaN */
-    }
-  else if(A <= in && in <= B)
-    {/* use the function log((1+g)/(1-g)) with g = x/(x + 2) */
-      return slnp1m1s(in / (in + 2));
-    }
-  else
-    {/* use the standard formula */
-      return slogs(in + 1);
-    }
-}
+/*

+ *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

+ *  Copyright (C) 2008-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

+ *

+ */

+

+#include "log1p.h"

+#include "log.h"

+#include "lnp1m1.h"

+

+float		slog1ps(float in) {

+  static double A = -1.0/3.0;

+  static double B = 0.5;

+

+  if(in < -1)

+    {/* got NaN */

+      return (in - in) / (in - in); /* NaN */

+    }

+  else if(A <= in && in <= B)

+    {/* use the function log((1+g)/(1-g)) with g = x/(x + 2) */

+      return slnp1m1s(in / (in + 2));

+    }

+  else

+    {/* use the standard formula */

+      return slogs(in + 1);

+    }

+}

diff --git a/src/c/elementaryFunctions/sqrt/csqrts.c b/src/c/elementaryFunctions/sqrt/csqrts.c
index a2ed819e..a24f9558 100644
--- a/src/c/elementaryFunctions/sqrt/csqrts.c
+++ b/src/c/elementaryFunctions/sqrt/csqrts.c
@@ -1,6 +1,7 @@
 /*

  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

  *  Copyright (C) 2008-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

diff --git a/src/c/elementaryFunctions/sqrt/zsqrts.c b/src/c/elementaryFunctions/sqrt/zsqrts.c
index f1a2b05c..3637ddd6 100644
--- a/src/c/elementaryFunctions/sqrt/zsqrts.c
+++ b/src/c/elementaryFunctions/sqrt/zsqrts.c
@@ -1,6 +1,7 @@
 /*

  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab

  *  Copyright (C) 2008-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

diff --git a/src/c/elementaryFunctions/tan/ztans.c b/src/c/elementaryFunctions/tan/ztans.c
index ff130dcc..761da36b 100644
--- a/src/c/elementaryFunctions/tan/ztans.c
+++ b/src/c/elementaryFunctions/tan/ztans.c
@@ -1,103 +1,104 @@
-/*
- *  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
- *
- */
-
-/*
-  ALGORITHM
-   based on the formula :
-
-                    0.5 sin(2 xr) +  i 0.5 sinh(2 xi)
-   tan(xr + i xi) = ---------------------------------
-                        cos(xr)^2  + sinh(xi)^2
-
-   noting  d = cos(xr)^2 + sinh(xi)^2, we have :
-
-           yr = 0.5 * sin(2 * xr) / d       (1)
-
-           yi = 0.5 * sinh(2 * xi) / d      (2)
-
-   to avoid spurious overflows in computing yi with
-   formula (2) (which results in NaN for yi)
-   we use also the following formula :
-
-           yi = sign(xi)   when |xi| > LIM  (3)
-
-   Explanations for (3) :
-
-   we have d = sinh(xi)^2 ( 1 + (cos(xr)/sinh(xi))^2 ),
-   so when :
-
-      (cos(xr)/sinh(xi))^2 < epsm   ( epsm = max relative error
-                                     for coding a real in a f.p.
-                                     number set F(b,p,emin,emax)
-                                        epsm = 0.5 b^(1-p) )
-   which is forced  when :
-
-       1/sinh(xi)^2 < epsm   (4)
-   <=> |xi| > asinh(1/sqrt(epsm)) (= 19.06... in ieee 754 double)
-
-   sinh(xi)^2 is a good approximation for d (relative to the f.p.
-   arithmetic used) and then yr may be approximate with :
-
-    yr = cosh(xi)/sinh(xi)
-       = sign(xi) (1 + exp(-2 |xi|))/(1 - exp(-2|xi|))
-       = sign(xi) (1 + 2 u + 2 u^2 + 2 u^3 + ...)
-
-   with u = exp(-2 |xi|)). Now when :
-
-       2 exp(-2|xi|) < epsm  (2)
-   <=> |xi| > 0.5 * log(2/epsm) (= 18.71... in ieee 754 double)
-
-   sign(xi)  is a good approximation for yr.
-
-   Constraint (1) is stronger than (2) and we take finaly
-
-   LIM = 1 + log(2/sqrt(epsm))
-
-   (log(2/sqrt(epsm)) being very near asinh(1/sqrt(epsm))
-
-AUTHOR
-   Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>
-*/
-
-#include <math.h>
-#include "lapack.h"
-#include "tan.h"
-#include "sqrt.h"
-#include "cos.h"
-#include "sinh.h"
-#include "sin.h"
-#include "log.h"
-#include "abs.h"
-
-#define localSign(x) x >= 0 ? 1.0 : -1.0
-
-doubleComplex		ztans(doubleComplex z) {
-  double Temp = 0;
-  double Lim = 1 + dlogs(2.0 / dsqrts( getRelativeMachinePrecision()));
-  double RealIn = zreals(z);
-  double ImagIn = zimags(z);
-  double RealOut = 0;
-  double ImagOut = 0;
-
-  Temp = pow(dcoss(RealIn), 2) + pow(dsinhs(ImagIn), 2);
-  RealOut = 0.5 * dsins(2 * RealIn) / Temp;
-  if(dabss(ImagIn) < Lim)
-    {
-      ImagOut = 0.5 * dsinhs(2 * ImagIn) / Temp;
-    }
-  else
-    {
-      ImagOut = localSign(ImagIn);
-    }
-
-  return DoubleComplex(RealOut, ImagOut);
-}
+/*

+ *  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

+ *

+ */

+

+/*

+  ALGORITHM

+   based on the formula :

+

+                    0.5 sin(2 xr) +  i 0.5 sinh(2 xi)

+   tan(xr + i xi) = ---------------------------------

+                        cos(xr)^2  + sinh(xi)^2

+

+   noting  d = cos(xr)^2 + sinh(xi)^2, we have :

+

+           yr = 0.5 * sin(2 * xr) / d       (1)

+

+           yi = 0.5 * sinh(2 * xi) / d      (2)

+

+   to avoid spurious overflows in computing yi with

+   formula (2) (which results in NaN for yi)

+   we use also the following formula :

+

+           yi = sign(xi)   when |xi| > LIM  (3)

+

+   Explanations for (3) :

+

+   we have d = sinh(xi)^2 ( 1 + (cos(xr)/sinh(xi))^2 ),

+   so when :

+

+      (cos(xr)/sinh(xi))^2 < epsm   ( epsm = max relative error

+                                     for coding a real in a f.p.

+                                     number set F(b,p,emin,emax)

+                                        epsm = 0.5 b^(1-p) )

+   which is forced  when :

+

+       1/sinh(xi)^2 < epsm   (4)

+   <=> |xi| > asinh(1/sqrt(epsm)) (= 19.06... in ieee 754 double)

+

+   sinh(xi)^2 is a good approximation for d (relative to the f.p.

+   arithmetic used) and then yr may be approximate with :

+

+    yr = cosh(xi)/sinh(xi)

+       = sign(xi) (1 + exp(-2 |xi|))/(1 - exp(-2|xi|))

+       = sign(xi) (1 + 2 u + 2 u^2 + 2 u^3 + ...)

+

+   with u = exp(-2 |xi|)). Now when :

+

+       2 exp(-2|xi|) < epsm  (2)

+   <=> |xi| > 0.5 * log(2/epsm) (= 18.71... in ieee 754 double)

+

+   sign(xi)  is a good approximation for yr.

+

+   Constraint (1) is stronger than (2) and we take finaly

+

+   LIM = 1 + log(2/sqrt(epsm))

+

+   (log(2/sqrt(epsm)) being very near asinh(1/sqrt(epsm))

+

+AUTHOR

+   Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>

+*/

+

+#include <math.h>

+#include "lapack.h"

+#include "tan.h"

+#include "sqrt.h"

+#include "cos.h"

+#include "sinh.h"

+#include "sin.h"

+#include "log.h"

+#include "abs.h"

+

+#define localSign(x) x >= 0 ? 1.0 : -1.0

+

+doubleComplex		ztans(doubleComplex z) {

+  double Temp = 0;

+  double Lim = 1 + dlogs(2.0 / dsqrts( getRelativeMachinePrecision()));

+  double RealIn = zreals(z);

+  double ImagIn = zimags(z);

+  double RealOut = 0;

+  double ImagOut = 0;

+

+  Temp = pow(dcoss(RealIn), 2) + pow(dsinhs(ImagIn), 2);

+  RealOut = 0.5 * dsins(2 * RealIn) / Temp;

+  if(dabss(ImagIn) < Lim)

+    {

+      ImagOut = 0.5 * dsinhs(2 * ImagIn) / Temp;

+    }

+  else

+    {

+      ImagOut = localSign(ImagIn);

+    }

+

+  return DoubleComplex(RealOut, ImagOut);

+}