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 | }
|
---|