Mercurial > hg > CbC > CbC_gcc
diff libquadmath/math/sqrtq.c @ 68:561a7518be6b
update gcc-4.6
author | Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> |
---|---|
date | Sun, 21 Aug 2011 07:07:55 +0900 |
parents | |
children | 04ced10e8804 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libquadmath/math/sqrtq.c Sun Aug 21 07:07:55 2011 +0900 @@ -0,0 +1,57 @@ +#include "quadmath-imp.h" +#include <math.h> +#include <float.h> + +__float128 +sqrtq (const __float128 x) +{ + __float128 y; + int exp; + + if (x == 0) + return x; + + if (isnanq (x)) + return x; + + if (x < 0) + return nanq (""); + + if (x <= DBL_MAX && x >= DBL_MIN) + { + /* Use double result as starting point. */ + y = sqrt ((double) x); + + /* Two Newton iterations. */ + y -= 0.5q * (y - x / y); + y -= 0.5q * (y - x / y); + return y; + } + +#ifdef HAVE_SQRTL + if (x <= LDBL_MAX && x >= LDBL_MIN) + { + /* Use long double result as starting point. */ + y = sqrtl ((long double) x); + + /* One Newton iteration. */ + y -= 0.5q * (y - x / y); + return y; + } +#endif + + /* If we're outside of the range of C types, we have to compute + the initial guess the hard way. */ + y = frexpq (x, &exp); + if (exp % 2) + y *= 2, exp--; + + y = sqrt (y); + y = scalbnq (y, exp / 2); + + /* Two Newton iterations. */ + y -= 0.5q * (y - x / y); + y -= 0.5q * (y - x / y); + return y; +} +