[337] | 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 | }
|
---|