]>
git.proxmox.com Git - wasi-libc.git/blob - libc-top-half/musl/src/math/sqrtf.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.c */
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
13 * ====================================================
18 static const float tiny
= 1.0e-30;
23 int32_t sign
= (int)0x80000000;
27 GET_FLOAT_WORD(ix
, x
);
29 /* take care of Inf and NaN */
30 if ((ix
&0x7f800000) == 0x7f800000)
31 return x
*x
+ x
; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
33 /* take care of zero */
36 return x
; /* sqrt(+-0) = +-0 */
38 return (x
-x
)/(x
-x
); /* sqrt(-ve) = sNaN */
42 if (m
== 0) { /* subnormal x */
43 for (i
= 0; (ix
&0x00800000) == 0; i
++)
47 m
-= 127; /* unbias exponent */
48 ix
= (ix
&0x007fffff)|0x00800000;
49 if (m
&1) /* odd m, double x to make it even */
51 m
>>= 1; /* m = [m/2] */
53 /* generate sqrt(x) bit by bit */
55 q
= s
= 0; /* q = sqrt(x) */
56 r
= 0x01000000; /* r = moving bit from right to left */
69 /* use floating add to find out rounding direction */
71 z
= 1.0f
- tiny
; /* raise inexact flag */
80 ix
= (q
>>1) + 0x3f000000;
82 SET_FLOAT_WORD(z
, ix
);