diff options
author | Sunil Shetye | 2019-05-16 12:18:48 +0530 |
---|---|---|
committer | Sunil Shetye | 2019-05-29 11:08:01 +0530 |
commit | 26b77d7593b5ee0792b6b556f5569ea4227c2b02 (patch) | |
tree | 8f92052234b01bf39b9c3a6e3cb12b3962d96b1b /src/c/elementaryFunctions/sqrt/csqrts.c | |
parent | 5a73e6bec4a12db7afae9de300e39256f754d8d3 (diff) | |
download | scilab2c-26b77d7593b5ee0792b6b556f5569ea4227c2b02.tar.gz scilab2c-26b77d7593b5ee0792b6b556f5569ea4227c2b02.tar.bz2 scilab2c-26b77d7593b5ee0792b6b556f5569ea4227c2b02.zip |
convert to unix format
Diffstat (limited to 'src/c/elementaryFunctions/sqrt/csqrts.c')
-rw-r--r-- | src/c/elementaryFunctions/sqrt/csqrts.c | 222 |
1 files changed, 111 insertions, 111 deletions
diff --git a/src/c/elementaryFunctions/sqrt/csqrts.c b/src/c/elementaryFunctions/sqrt/csqrts.c index a24f9558..07ef3bb0 100644 --- a/src/c/elementaryFunctions/sqrt/csqrts.c +++ b/src/c/elementaryFunctions/sqrt/csqrts.c @@ -1,111 +1,111 @@ -/*
- * 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 <math.h>
-#include "sqrt.h"
-#include "lapack.h"
-#include "abs.h"
-#include "sign.h"
-#include "pythag.h"
-
-#ifdef _MSC_VER
-#include <float.h>
-#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);
-}
+/* + * 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 <math.h> +#include "sqrt.h" +#include "lapack.h" +#include "abs.h" +#include "sign.h" +#include "pythag.h" + +#ifdef _MSC_VER +#include <float.h> +#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); +} |