1 /*-
2  * SPDX-License-Identifier: BSD-2-Clause
35 /* Return (x + ulp) for normal positive x. Assumes no overflow. */
42 	if (++u.bits.manl == 0) {  in inc()
43 		if (++u.bits.manh == 0) {  in inc()
51 /* Return (x - ulp) for normal positive x. Assumes no underflow. */
58 	if (u.bits.manl-- == 0) {  in dec()
59 		if (u.bits.manh-- == LDBL_NBIT) {  in dec()
60 			u.bits.exp--;  in dec()
86 	/* If x = -Inf, then sqrt(x) = NaN. */  in sqrtl()
87 	if (u.bits.exp == LDBL_MAX_EXP * 2 - 1)  in sqrtl()
90 	/* If x = +-0, then sqrt(x) = +-0. */  in sqrtl()
91 	if ((u.bits.manh | u.bits.manl | u.bits.exp) == 0)  in sqrtl()
94 	/* If x < 0, then raise invalid and return NaN */  in sqrtl()
96 		return ((x - x) / (x - x));  in sqrtl()
100 	if (u.bits.exp == 0) {  in sqrtl()
102 		u.e *= 0x1.0p514;  in sqrtl()
103 		k = -514;  in sqrtl()
105 		k = 0;  in sqrtl()
111 	if ((u.bits.exp - 0x3ffe) & 1) {	/* n is odd.     */  in sqrtl()
112 		k += u.bits.exp - 0x3fff;	/* 2k = n - 1.   */  in sqrtl()
113 		u.bits.exp = 0x3fff;		/* u.e in [1,2). */  in sqrtl()
115 		k += u.bits.exp - 0x4000;	/* 2k = n - 2.   */  in sqrtl()
116 		u.bits.exp = 0x4000;		/* u.e in [2,4). */  in sqrtl()
123 	xn = sqrt(u.e);			/* 53-bit estimate of sqrtl(x). */  in sqrtl()
125 	xn = (xn + (u.e / xn)) * 0.5;	/* 106-bit estimate. */  in sqrtl()
128 	u.bits.manl = 0;		/* Zero out lower bits. */  in sqrtl()
129 	lo = (lo - u.e) / xn;		/* Low bits divided by xn. */  in sqrtl()
132 	u.bits.exp += (k >> 1) - 1;  in sqrtl()
136 	fesetround(FE_TOWARDZERO);	/* Set to round-toward-zero. */  in sqrtl()
144 		/* Round correctly for inputs like x = y**2 - ulp. */  in sqrtl()
145 		xn = dec(xn);		/* xn = xn - ulp. */  in sqrtl()
149 		xn = inc(xn);		/* xn = xn + ulp. */  in sqrtl()
151 		u.e = inc(u.e);		/* u.e = u.e + ulp. */  in sqrtl()
152 		xn = inc(xn);		/* xn  = xn + ulp. */  in sqrtl()
156 	u.bits.exp--;  in sqrtl()