* Add Scilab square root algorithm.

* Need some further testing

10 files changed, 254 insertions, 69 deletions

diff --git a/src/elementaryFunctions/sqrt/Makefile.am b/src/elementaryFunctions/sqrt/Makefile.am
index 94e096c3..6bd9e6d5 100644
--- a/src/elementaryFunctions/sqrt/Makefile.am
+++ b/src/elementaryFunctions/sqrt/Makefile.am
@@ -10,8 +10,9 @@
 ##
 ##
 
-libSqrt_la_CFLAGS =	-I ../../type \
-		-I ../includes
+libSqrt_la_CFLAGS =	-I  $(top_builddir)/type \
+			-I $(top_builddir)/elementaryFunctions/includes \
+			-I $(top_builddir)/auxiliaryFunctions/includes
 
 instdir = $(top_builddir)/lib
 
@@ -39,7 +40,10 @@ check_INCLUDES	=	-I $(top_builddir)/elementaryFunctions/includes \
 
 check_LDADD	=	$(top_builddir)/type/libDoubleComplex.la \
 			$(top_builddir)/type/libFloatComplex.la \
+			$(top_builddir)/lib/lapack/libscilapack.la \
 			$(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
+			$(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
+			$(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
 			@LIBMATH@
 
 check_PROGRAMS	= 	testFloatSqrt testDoubleSqrt
diff --git a/src/elementaryFunctions/sqrt/Makefile.in b/src/elementaryFunctions/sqrt/Makefile.in
index 99a95b27..11f6e142 100644
--- a/src/elementaryFunctions/sqrt/Makefile.in
+++ b/src/elementaryFunctions/sqrt/Makefile.in
@@ -66,7 +66,10 @@ am_testDoubleSqrt_OBJECTS = testDoubleSqrt-testDoubleSqrt.$(OBJEXT)
 testDoubleSqrt_OBJECTS = $(am_testDoubleSqrt_OBJECTS)
 am__DEPENDENCIES_1 = $(top_builddir)/type/libDoubleComplex.la \
 	$(top_builddir)/type/libFloatComplex.la \
-	$(top_builddir)/elementaryFunctions/sqrt/libSqrt.la
+	$(top_builddir)/lib/lapack/libscilapack.la \
+	$(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
+	$(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
+	$(top_builddir)/auxiliaryFunctions/pythag/libPythag.la
 testDoubleSqrt_DEPENDENCIES = $(am__DEPENDENCIES_1)
 testDoubleSqrt_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
 	$(LIBTOOLFLAGS) --mode=link $(CCLD) $(testDoubleSqrt_CFLAGS) \
@@ -206,8 +209,9 @@ sysconfdir = @sysconfdir@
 target_alias = @target_alias@
 top_builddir = @top_builddir@
 top_srcdir = @top_srcdir@
-libSqrt_la_CFLAGS = -I ../../type \
-		-I ../includes
+libSqrt_la_CFLAGS = -I  $(top_builddir)/type \
+			-I $(top_builddir)/elementaryFunctions/includes \
+			-I $(top_builddir)/auxiliaryFunctions/includes
 
 instdir = $(top_builddir)/lib
 pkglib_LTLIBRARIES = libSqrt.la
@@ -231,7 +235,10 @@ check_INCLUDES = -I $(top_builddir)/elementaryFunctions/includes \
 
 check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
 			$(top_builddir)/type/libFloatComplex.la \
+			$(top_builddir)/lib/lapack/libscilapack.la \
 			$(top_builddir)/elementaryFunctions/sqrt/libSqrt.la \
+			$(top_builddir)/auxiliaryFunctions/abs/libAbs.la \
+			$(top_builddir)/auxiliaryFunctions/pythag/libPythag.la \
 			@LIBMATH@
 
 
diff --git a/src/elementaryFunctions/sqrt/csqrta.c b/src/elementaryFunctions/sqrt/csqrta.c
index 62851ff1..fc9d36cd 100644
--- a/src/elementaryFunctions/sqrt/csqrta.c
+++ b/src/elementaryFunctions/sqrt/csqrta.c
@@ -1,14 +1,14 @@
 /*
-**  -*- C -*-
-**
-** csqrta.c
-** Made by  Bruno JOFRET <bruno.jofret@inria.fr>
-**
-** Started on  Fri Sep  7 11:40:55 2007 bruno
-** Last update Mon Oct 22 09:49:44 2007 bruno
-**
-** Copyright INRIA 2007
-*/
+ *  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 "sqrt.h"
 
diff --git a/src/elementaryFunctions/sqrt/csqrts.c b/src/elementaryFunctions/sqrt/csqrts.c
index f39b67f5..da7000e8 100644
--- a/src/elementaryFunctions/sqrt/csqrts.c
+++ b/src/elementaryFunctions/sqrt/csqrts.c
@@ -10,9 +10,96 @@
  *
  */
 
+#include <math.h>
 #include "sqrt.h"
+#include "lapack.h"
+#include "abs.h"
+#include "sign.h"
+#include "pythag.h"
+
+#define _sign(a, b) b >=0 ? a : -a
 
 floatComplex	csqrts(floatComplex in) {
-  /* FIXME : Dummy ... */
-  return in;
+  float RMax =  (float) getOverflowThreshold();
+  float BRMin = 2.0f * (float) getUnderflowThreshold();
+
+  float RealIn = creals(in);
+  float ImgIn = cimags(in);
+
+  float RealOut = 0;
+  float ImgOut = 0;
+
+  if(RealIn == 0)
+    {/* pure imaginary case */
+      if(dabss(ImgIn >= BRMin))
+	RealOut = ssqrts(0.5f * sabss(ImgIn));
+      else
+	RealOut = ssqrts(sabss(ImgIn)) * ssqrts(0.5);
+
+      ImgOut = _sign(1, ImgIn) * RealOut;
+    }
+  else if( sabss(RealIn) <= RMax && sabss(ImgIn) <= RMax)
+    {/* standard case : a (not zero) and b are finite */
+      float Temp = ssqrts(2.0f * (sabss(RealIn) + spythags(RealIn, ImgIn)));
+      /* overflow test */
+      if(Temp > RMax)
+	{/*     handle (spurious) overflow by scaling a and b */
+	  float RealTemp = RealIn / 16.0f;
+	  float ImgTemp = ImgIn / 16.0f;
+	  Temp = ssqrts(2.0f * (sabss(RealIn) + spythags(RealIn, ImgTemp)));
+	  if(RealTemp >= 0)
+	    {
+	      RealOut	= 2 * Temp;
+	      ImgOut	= 4 * ImgTemp / Temp;
+	    }
+	  else
+	    {
+	      RealOut	= 4 * sabss(ImgIn) / Temp;
+	      ImgOut	= _sign(2, ImgIn) * Temp;
+	    }
+	}
+      else if(RealIn >= 0) /* classic switch to get the stable formulas */
+	{
+	  RealOut	= 0.5f * Temp;
+	  ImgOut	= ImgIn / Temp;
+	}
+      else
+	{
+	  RealOut	= sabss(ImgIn) / Temp;
+	  ImgOut	= _sign(0.5f, ImgIn) * Temp;
+	}
+    }
+  else
+    {
+      /*
+      //Here we treat the special cases where a and b are +- 00 or NaN.
+      //The following is the treatment recommended by the C99 standard
+      //with the simplification of returning NaN + i NaN if the
+      //the real part or the imaginary part is NaN (C99 recommends
+      //something more complicated)
+      */
+
+      if(isnan(RealIn) == 1 || isnan(ImgIn) == 1)
+	{/* got NaN + i NaN */
+	  RealOut	= RealIn + ImgIn;
+	  ImgOut	= RealOut;
+	}
+      else if( dabss(ImgIn) > RMax)
+	{/* case a +- i oo -> result must be +oo +- i oo  for all a (finite or not) */
+	  RealOut	= sabss(ImgIn);
+	  ImgOut	= ImgIn;
+	}
+      else if(RealIn < -RMax)
+	{/* here a is -Inf and b is finite */
+	  RealOut	= 0;
+	  ImgOut	= _sign(1, ImgIn) * sabss(RealIn);
+	}
+      else
+	{/* here a is +Inf and b is finite */
+	  RealOut	= RealIn;
+	  ImgOut	= 0;
+	}
+    }
+
+  return FloatComplex(RealOut, ImgOut);
 }
diff --git a/src/elementaryFunctions/sqrt/dsqrta.c b/src/elementaryFunctions/sqrt/dsqrta.c
index 56d6ea1f..a948f35b 100644
--- a/src/elementaryFunctions/sqrt/dsqrta.c
+++ b/src/elementaryFunctions/sqrt/dsqrta.c
@@ -1,14 +1,14 @@
 /*
-**  -*- C -*-
-**
-** dsqrta.c
-** Made by  Bruno JOFRET <bruno.jofret@inria.fr>
-**
-** Started on  Fri Sep  7 11:40:21 2007 bruno
-** Last update Mon Oct 22 09:49:34 2007 bruno
-**
-** Copyright INRIA 2007
-*/
+ *  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 "sqrt.h"
 
diff --git a/src/elementaryFunctions/sqrt/dsqrts.c b/src/elementaryFunctions/sqrt/dsqrts.c
index a8159e38..350cbbd0 100644
--- a/src/elementaryFunctions/sqrt/dsqrts.c
+++ b/src/elementaryFunctions/sqrt/dsqrts.c
@@ -1,14 +1,14 @@
 /*
-**  -*- C -*-
-**
-** dsqrts.c
-** Made by  Bruno JOFRET <bruno.jofret@inria.fr>
-**
-** Started on  Fri Apr 20 10:59:31 2007 jofret
-** Last update Mon Apr 23 17:37:02 2007 jofret
-**
-** Copyright INRIA 2007
-*/
+ *  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 "sqrt.h"
diff --git a/src/elementaryFunctions/sqrt/ssqrta.c b/src/elementaryFunctions/sqrt/ssqrta.c
index 6ffa5c26..9c379751 100644
--- a/src/elementaryFunctions/sqrt/ssqrta.c
+++ b/src/elementaryFunctions/sqrt/ssqrta.c
@@ -1,14 +1,14 @@
 /*
-**  -*- C -*-
-**
-** ssqrta.c
-** Made by  Bruno JOFRET <bruno.jofret@inria.fr>
-**
-** Started on  Fri Sep  7 11:06:03 2007 bruno
-** Last update Mon Oct 22 09:49:22 2007 bruno
-**
-** Copyright INRIA 2007
-*/
+ *  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 "sqrt.h"
 
diff --git a/src/elementaryFunctions/sqrt/testSqrt.h b/src/elementaryFunctions/sqrt/testSqrt.h
index 5761d293..2b9c3de2 100644
--- a/src/elementaryFunctions/sqrt/testSqrt.h
+++ b/src/elementaryFunctions/sqrt/testSqrt.h
@@ -1,14 +1,14 @@
 /*
-**  -*- C -*-
-**
-** testSqrt.h
-** Made by  Bruno JOFRET <bruno.jofret@inria.fr>
-**
-** Started on  Fri Mar 30 12:01:28 2007 jofret
-** Last update Fri Apr 20 10:52:56 2007 jofret
-**
-** Copyright INRIA 2007
-*/
+ *  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
+ *
+ */
 
 #ifndef _TESTSQRT_H_
 #define _TESTSQRT_H_
diff --git a/src/elementaryFunctions/sqrt/zsqrta.c b/src/elementaryFunctions/sqrt/zsqrta.c
index 92e6fb1c..b6472951 100644
--- a/src/elementaryFunctions/sqrt/zsqrta.c
+++ b/src/elementaryFunctions/sqrt/zsqrta.c
@@ -1,14 +1,14 @@
 /*
-**  -*- C -*-
-**
-** zsqrta.c
-** Made by  Bruno JOFRET <bruno.jofret@inria.fr>
-**
-** Started on  Fri Sep  7 11:41:40 2007 bruno
-** Last update Mon Oct 22 09:49:12 2007 bruno
-**
-** Copyright INRIA 2007
-*/
+ *  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 "sqrt.h"
 
diff --git a/src/elementaryFunctions/sqrt/zsqrts.c b/src/elementaryFunctions/sqrt/zsqrts.c
index 05d9192e..357a752e 100644
--- a/src/elementaryFunctions/sqrt/zsqrts.c
+++ b/src/elementaryFunctions/sqrt/zsqrts.c
@@ -10,9 +10,96 @@
  *
  */
 
+#include <math.h>
 #include "sqrt.h"
+#include "lapack.h"
+#include "abs.h"
+#include "sign.h"
+#include "pythag.h"
+
+#define _sign(a, b) b >=0 ? a : -a
 
 doubleComplex	zsqrts(doubleComplex in) {
-  /* FIXME : Dummy ... */
-  return in;
+  double RMax =  getOverflowThreshold();
+  double BRMin = 2 * getUnderflowThreshold();
+
+  double RealIn = zreals(in);
+  double ImgIn = zimags(in);
+
+  double RealOut = 0;
+  double ImgOut = 0;
+
+  if(RealIn == 0)
+    {/* pure imaginary case */
+      if(dabss(ImgIn >= BRMin))
+	RealOut = dsqrts(0.5 * dabss(ImgIn));
+      else
+	RealOut = dsqrts(dabss(ImgIn)) * dsqrts(0.5);
+
+      ImgOut = _sign(1, ImgIn) * RealOut;
+    }
+  else if( dabss(RealIn) <= RMax && dabss(ImgIn) <= RMax)
+    {/* standard case : a (not zero) and b are finite */
+      double Temp = dsqrts(2 * (dabss(RealIn) + dpythags(RealIn, ImgIn)));
+      /* overflow test */
+      if(Temp > RMax)
+	{/*     handle (spurious) overflow by scaling a and b */
+	  double RealTemp = RealIn / 16;
+	  double ImgTemp = ImgIn / 16;
+	  Temp = dsqrts(2 * (dabss(RealIn) + dpythags(RealIn, ImgTemp)));
+	  if(RealTemp >= 0)
+	    {
+	      RealOut	= 2 * Temp;
+	      ImgOut	= 4 * ImgTemp / Temp;
+	    }
+	  else
+	    {
+	      RealOut	= 4 * dabss(ImgIn) / Temp;
+	      ImgOut	= _sign(2, ImgIn) * Temp;
+	    }
+	}
+      else if(RealIn >= 0) /* classic switch to get the stable formulas */
+	{
+	  RealOut	= 0.5 * Temp;
+	  ImgOut	= ImgIn / Temp;
+	}
+      else
+	{
+	  RealOut	= dabss(ImgIn) / Temp;
+	  ImgOut	= _sign(0.5, ImgIn) * Temp;
+	}
+    }
+  else
+    {
+      /*
+      //Here we treat the special cases where a and b are +- 00 or NaN.
+      //The following is the treatment recommended by the C99 standard
+      //with the simplification of returning NaN + i NaN if the
+      //the real part or the imaginary part is NaN (C99 recommends
+      //something more complicated)
+      */
+
+      if(isnan(RealIn) == 1 || isnan(ImgIn) == 1)
+	{/* got NaN + i NaN */
+	  RealOut	= RealIn + ImgIn;
+	  ImgOut	= RealOut;
+	}
+      else if( dabss(ImgIn) > RMax)
+	{/* case a +- i oo -> result must be +oo +- i oo  for all a (finite or not) */
+	  RealOut	= dabss(ImgIn);
+	  ImgOut	= ImgIn;
+	}
+      else if(RealIn < -RMax)
+	{/* here a is -Inf and b is finite */
+	  RealOut	= 0;
+	  ImgOut	= _sign(1, ImgIn) * dabss(RealIn);
+	}
+      else
+	{/* here a is +Inf and b is finite */
+	  RealOut	= RealIn;
+	  ImgOut	= 0;
+	}
+    }
+
+  return DoubleComplex(RealOut, ImgOut);
 }