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