737df6e617
2012-10-31 Tobias Burnus <burnus@net-b.de>
Joseph Myers <joseph@codesourcery.com>
David S. Miller <davem@davemloft.net>
Ulrich Drepper <drepper@redhat.com>
Marek Polacek <polacek@redhat.com>:
Petr Baudis <pasky@suse.cz>
* math/complex.c (csqrtq): NaN and INF fixes.
* math/sqrtq.c (sqrt): NaN, INF and < 0 fixes.
* math/expm1q.c (expm1q): Changes from GLIBC. Use expq for
large parameters. Fix errno for boundary conditions.
* math/finiteq.c (finiteq): Add comment.
* math/fmaq.c (fmaq): Changes from GLIBC. Fix missing underflows
and bad results for some subnormal results. Fix sign of inexact
zero return. Fix sign of exact zero return.
Ensure additions are not scheduled after fetestexcept.
* math/jnq.c (jnq): Changes from GLIBC. Set up errno properly
for ynq. Fix jnq precision.
* math/nearbyintq.c (nearbyintq): Changes from GLIBC. Do not
manipulate bits before adding and subtracting TWO112[sx].
* math/rintq.c (rintq): Ditto.
* math/scalbnq.c (scalbnq): Changes from GLIBC. Fix integer
overflow.
Co-Authored-By: David S. Miller <davem@davemloft.net>
Co-Authored-By: Joseph Myers <joseph@codesourcery.com>
Co-Authored-By: Ulrich Drepper <drepper@redhat.com>
From-SVN: r193037
61 lines
1.1 KiB
C
61 lines
1.1 KiB
C
#include "quadmath-imp.h"
|
|
#include <math.h>
|
|
#include <float.h>
|
|
|
|
__float128
|
|
sqrtq (const __float128 x)
|
|
{
|
|
__float128 y;
|
|
int exp;
|
|
|
|
if (isnanq (x) || (isinfq (x) && x > 0))
|
|
return x;
|
|
|
|
if (x == 0)
|
|
return x;
|
|
|
|
if (x < 0)
|
|
{
|
|
/* Return NaN with invalid signal. */
|
|
return (x - x) / (x - x);
|
|
}
|
|
|
|
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;
|
|
}
|
|
|