e_sqrtf.c revision 8870
1251883Speter/* e_sqrtf.c -- float version of e_sqrt.c. 2251883Speter * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3251883Speter */ 4289166Speter 5251883Speter/* 6289166Speter * ==================================================== 7289166Speter * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8251883Speter * 9251883Speter * Developed at SunPro, a Sun Microsystems, Inc. business. 10251883Speter * Permission to use, copy, modify, and distribute this 11251883Speter * software is freely granted, provided that this notice 12251883Speter * is preserved. 13251883Speter * ==================================================== 14251883Speter */ 15251883Speter 16251883Speter#ifndef lint 17251883Speterstatic char rcsid[] = "$Id: e_sqrtf.c,v 1.1.1.1 1994/08/19 09:39:57 jkh Exp $"; 18251883Speter#endif 19251883Speter 20251883Speter#include "math.h" 21289166Speter#include "math_private.h" 22251883Speter 23251883Speter#ifdef __STDC__ 24251883Speterstatic const float one = 1.0, tiny=1.0e-30; 25251883Speter#else 26251883Speterstatic float one = 1.0, tiny=1.0e-30; 27251883Speter#endif 28251883Speter 29251883Speter#ifdef __STDC__ 30251883Speter float __ieee754_sqrtf(float x) 31251883Speter#else 32251883Speter float __ieee754_sqrtf(x) 33251883Speter float x; 34289166Speter#endif 35289166Speter{ 36251883Speter float z; 37251883Speter int32_t sign = (int)0x80000000; 38251883Speter int32_t ix,s,q,m,t,i; 39251883Speter u_int32_t r; 40251883Speter 41251883Speter GET_FLOAT_WORD(ix,x); 42251883Speter 43251883Speter /* take care of Inf and NaN */ 44251883Speter if((ix&0x7f800000)==0x7f800000) { 45251883Speter return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf 46251883Speter sqrt(-inf)=sNaN */ 47289166Speter } 48251883Speter /* take care of zero */ 49251883Speter if(ix<=0) { 50251883Speter if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */ 51251883Speter else if(ix<0) 52251883Speter return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ 53251883Speter } 54251883Speter /* normalize x */ 55251883Speter m = (ix>>23); 56251883Speter if(m==0) { /* subnormal x */ 57251883Speter for(i=0;(ix&0x00800000)==0;i++) ix<<=1; 58251883Speter m -= i-1; 59251883Speter } 60251883Speter m -= 127; /* unbias exponent */ 61251883Speter ix = (ix&0x007fffff)|0x00800000; 62251883Speter if(m&1) /* odd m, double x to make it even */ 63251883Speter ix += ix; 64251883Speter m >>= 1; /* m = [m/2] */ 65251883Speter 66251883Speter /* generate sqrt(x) bit by bit */ 67251883Speter ix += ix; 68251883Speter q = s = 0; /* q = sqrt(x) */ 69251883Speter r = 0x01000000; /* r = moving bit from right to left */ 70251883Speter 71251883Speter while(r!=0) { 72251883Speter t = s+r; 73251883Speter if(t<=ix) { 74251883Speter s = t+r; 75251883Speter ix -= t; 76251883Speter q += r; 77251883Speter } 78251883Speter ix += ix; 79251883Speter r>>=1; 80289166Speter } 81251883Speter 82251883Speter /* use floating add to find out rounding direction */ 83251883Speter if(ix!=0) { 84251883Speter z = one-tiny; /* trigger inexact flag */ 85251883Speter if (z>=one) { 86251883Speter z = one+tiny; 87251883Speter if (z>one) 88251883Speter q += 2; 89289166Speter else 90289166Speter q += (q&1); 91289166Speter } 92251883Speter } 93251883Speter ix = (q>>1)+0x3f000000; 94251883Speter ix += (m <<23); 95251883Speter SET_FLOAT_WORD(z,ix); 96251883Speter return z; 97251883Speter} 98251883Speter