1#include "quadmath-imp.h" 2#include <math.h> 3#include <float.h> 4 5__float128 6sqrtq (const __float128 x) 7{ 8 __float128 y; 9 int exp; 10 11 if (isnanq (x) || (isinfq (x) && x > 0)) 12 return x; 13 14 if (x == 0) 15 return x; 16 17 if (x < 0) 18 { 19 /* Return NaN with invalid signal. */ 20 return (x - x) / (x - x); 21 } 22 23 if (x <= DBL_MAX && x >= DBL_MIN) 24 { 25 /* Use double result as starting point. */ 26 y = sqrt ((double) x); 27 28 /* Two Newton iterations. */ 29 y -= 0.5q * (y - x / y); 30 y -= 0.5q * (y - x / y); 31 return y; 32 } 33 34#ifdef HAVE_SQRTL 35 if (x <= LDBL_MAX && x >= LDBL_MIN) 36 { 37 /* Use long double result as starting point. */ 38 y = sqrtl ((long double) x); 39 40 /* One Newton iteration. */ 41 y -= 0.5q * (y - x / y); 42 return y; 43 } 44#endif 45 46 /* If we're outside of the range of C types, we have to compute 47 the initial guess the hard way. */ 48 y = frexpq (x, &exp); 49 if (exp % 2) 50 y *= 2, exp--; 51 52 y = sqrt (y); 53 y = scalbnq (y, exp / 2); 54 55 /* Two Newton iterations. */ 56 y -= 0.5q * (y - x / y); 57 y -= 0.5q * (y - x / y); 58 return y; 59} 60 61