github.com/afumu/libc@v0.0.6/musl/src/math/sqrtf.c (about)

     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  	SET_FLOAT_WORD(z, ix + ((uint32_t)m << 23));
    82  	return z;
    83  }