1217309Snwhitehorn/* e_sqrtf.c -- float version of e_sqrt.c. 2251843Sbapt * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3217309Snwhitehorn */ 4217309Snwhitehorn 5217309Snwhitehorn/* 6251843Sbapt * ==================================================== 7217309Snwhitehorn * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8217309Snwhitehorn * 9217309Snwhitehorn * Developed at SunPro, a Sun Microsystems, Inc. business. 10217309Snwhitehorn * Permission to use, copy, modify, and distribute this 11217309Snwhitehorn * software is freely granted, provided that this notice 12217309Snwhitehorn * is preserved. 13217309Snwhitehorn * ==================================================== 14217309Snwhitehorn */ 15217309Snwhitehorn 16217309Snwhitehorn#ifndef lint 17217309Snwhitehornstatic char rcsid[] = "$FreeBSD$"; 18217309Snwhitehorn#endif 19217309Snwhitehorn 20217309Snwhitehorn#include "math.h" 21217309Snwhitehorn#include "math_private.h" 22217309Snwhitehorn 23217309Snwhitehornstatic const float one = 1.0, tiny=1.0e-30; 24217309Snwhitehorn 25217309Snwhitehornfloat 26217309Snwhitehorn__ieee754_sqrtf(float x) 27217309Snwhitehorn{ 28217309Snwhitehorn float z; 29217309Snwhitehorn int32_t sign = (int)0x80000000; 30217309Snwhitehorn int32_t ix,s,q,m,t,i; 31217309Snwhitehorn u_int32_t r; 32217309Snwhitehorn 33217309Snwhitehorn GET_FLOAT_WORD(ix,x); 34217309Snwhitehorn 35217309Snwhitehorn /* take care of Inf and NaN */ 36217309Snwhitehorn if((ix&0x7f800000)==0x7f800000) { 37217309Snwhitehorn return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf 38217309Snwhitehorn sqrt(-inf)=sNaN */ 39217309Snwhitehorn } 40217309Snwhitehorn /* take care of zero */ 41217309Snwhitehorn if(ix<=0) { 42217309Snwhitehorn if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */ 43217309Snwhitehorn else if(ix<0) 44217309Snwhitehorn return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ 45217309Snwhitehorn } 46217309Snwhitehorn /* normalize x */ 47217309Snwhitehorn m = (ix>>23); 48217309Snwhitehorn if(m==0) { /* subnormal x */ 49217309Snwhitehorn for(i=0;(ix&0x00800000)==0;i++) ix<<=1; 50217309Snwhitehorn m -= i-1; 51217309Snwhitehorn } 52217309Snwhitehorn m -= 127; /* unbias exponent */ 53217309Snwhitehorn ix = (ix&0x007fffff)|0x00800000; 54217309Snwhitehorn if(m&1) /* odd m, double x to make it even */ 55217309Snwhitehorn ix += ix; 56217309Snwhitehorn m >>= 1; /* m = [m/2] */ 57217309Snwhitehorn 58217309Snwhitehorn /* generate sqrt(x) bit by bit */ 59217309Snwhitehorn ix += ix; 60217309Snwhitehorn q = s = 0; /* q = sqrt(x) */ 61217309Snwhitehorn r = 0x01000000; /* r = moving bit from right to left */ 62217309Snwhitehorn 63217309Snwhitehorn while(r!=0) { 64217309Snwhitehorn t = s+r; 65217309Snwhitehorn if(t<=ix) { 66217309Snwhitehorn s = t+r; 67217309Snwhitehorn ix -= t; 68217309Snwhitehorn q += r; 69217309Snwhitehorn } 70217309Snwhitehorn ix += ix; 71217309Snwhitehorn r>>=1; 72217309Snwhitehorn } 73217309Snwhitehorn 74217309Snwhitehorn /* use floating add to find out rounding direction */ 75217309Snwhitehorn if(ix!=0) { 76217309Snwhitehorn z = one-tiny; /* trigger inexact flag */ 77217309Snwhitehorn if (z>=one) { 78217309Snwhitehorn z = one+tiny; 79217309Snwhitehorn if (z>one) 80217309Snwhitehorn q += 2; 81217309Snwhitehorn else 82217309Snwhitehorn q += (q&1); 83217309Snwhitehorn } 84217309Snwhitehorn } 85217309Snwhitehorn ix = (q>>1)+0x3f000000; 86217309Snwhitehorn ix += (m <<23); 87217309Snwhitehorn SET_FLOAT_WORD(z,ix); 88217309Snwhitehorn return z; 89217309Snwhitehorn} 90217309Snwhitehorn