math: fix pow signed shift ub
[oweals/musl.git] / src / math / sqrtf.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.c */
2 /*
3  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
4  */
5 /*
6  * ====================================================
7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8  *
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
12  * is preserved.
13  * ====================================================
14  */
15
16 #include "libm.h"
17
18 static const float tiny = 1.0e-30;
19
20 float sqrtf(float x)
21 {
22         float z;
23         int32_t sign = (int)0x80000000;
24         int32_t ix,s,q,m,t,i;
25         uint32_t r;
26
27         GET_FLOAT_WORD(ix, x);
28
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 */
32
33         /* take care of zero */
34         if (ix <= 0) {
35                 if ((ix&~sign) == 0)
36                         return x;  /* sqrt(+-0) = +-0 */
37                 if (ix < 0)
38                         return (x-x)/(x-x);  /* sqrt(-ve) = sNaN */
39         }
40         /* normalize x */
41         m = ix>>23;
42         if (m == 0) {  /* subnormal x */
43                 for (i = 0; (ix&0x00800000) == 0; i++)
44                         ix<<=1;
45                 m -= i - 1;
46         }
47         m -= 127;  /* unbias exponent */
48         ix = (ix&0x007fffff)|0x00800000;
49         if (m&1)  /* odd m, double x to make it even */
50                 ix += ix;
51         m >>= 1;  /* m = [m/2] */
52
53         /* generate sqrt(x) bit by bit */
54         ix += ix;
55         q = s = 0;       /* q = sqrt(x) */
56         r = 0x01000000;  /* r = moving bit from right to left */
57
58         while (r != 0) {
59                 t = s + r;
60                 if (t <= ix) {
61                         s = t+r;
62                         ix -= t;
63                         q += r;
64                 }
65                 ix += ix;
66                 r >>= 1;
67         }
68
69         /* use floating add to find out rounding direction */
70         if (ix != 0) {
71                 z = 1.0f - tiny; /* raise inexact flag */
72                 if (z >= 1.0f) {
73                         z = 1.0f + tiny;
74                         if (z > 1.0f)
75                                 q += 2;
76                         else
77                                 q += q & 1;
78                 }
79         }
80         ix = (q>>1) + 0x3f000000;
81         ix += m << 23;
82         SET_FLOAT_WORD(z, ix);
83         return z;
84 }