From 879e2ac0a540fa1b199e20d47830aa5eea484a4c Mon Sep 17 00:00:00 2001 From: cornet Date: Wed, 22 Apr 2009 05:54:43 +0000 Subject: build with visual studio (dynamic libraries) --- src/elementaryFunctions/sqrt/csqrts.c | 215 +++++++++++++++++----------------- 1 file changed, 110 insertions(+), 105 deletions(-) (limited to 'src/elementaryFunctions/sqrt/csqrts.c') diff --git a/src/elementaryFunctions/sqrt/csqrts.c b/src/elementaryFunctions/sqrt/csqrts.c index f06d2dd1..a2ed819e 100644 --- a/src/elementaryFunctions/sqrt/csqrts.c +++ b/src/elementaryFunctions/sqrt/csqrts.c @@ -1,105 +1,110 @@ -/* - * 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 -#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) { - 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); -} +/* + * 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 +#include "sqrt.h" +#include "lapack.h" +#include "abs.h" +#include "sign.h" +#include "pythag.h" + +#ifdef _MSC_VER +#include +#define isnan(x) _isnan((double)x) +#endif + +#define _sign(a, b) b >=0 ? a : -a + +floatComplex csqrts(floatComplex 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); +} -- cgit