12116Sjkh/* e_sqrtf.c -- float version of e_sqrt.c. 22116Sjkh * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 32116Sjkh */ 42116Sjkh 52116Sjkh/* 62116Sjkh * ==================================================== 72116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 82116Sjkh * 92116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 102116Sjkh * Permission to use, copy, modify, and distribute this 118870Srgrimes * software is freely granted, provided that this notice 122116Sjkh * is preserved. 132116Sjkh * ==================================================== 142116Sjkh */ 152116Sjkh 162116Sjkh#ifndef lint 1750476Speterstatic char rcsid[] = "$FreeBSD$"; 182116Sjkh#endif 192116Sjkh 202116Sjkh#include "math.h" 212116Sjkh#include "math_private.h" 222116Sjkh 232116Sjkhstatic const float one = 1.0, tiny=1.0e-30; 242116Sjkh 2597413Salfredfloat 2697413Salfred__ieee754_sqrtf(float x) 272116Sjkh{ 282116Sjkh float z; 298870Srgrimes int32_t sign = (int)0x80000000; 302116Sjkh int32_t ix,s,q,m,t,i; 312116Sjkh u_int32_t r; 322116Sjkh 332116Sjkh GET_FLOAT_WORD(ix,x); 342116Sjkh 352116Sjkh /* take care of Inf and NaN */ 368870Srgrimes if((ix&0x7f800000)==0x7f800000) { 372116Sjkh return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf 382116Sjkh sqrt(-inf)=sNaN */ 398870Srgrimes } 402116Sjkh /* take care of zero */ 412116Sjkh if(ix<=0) { 422116Sjkh if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */ 432116Sjkh else if(ix<0) 442116Sjkh return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ 452116Sjkh } 462116Sjkh /* normalize x */ 472116Sjkh m = (ix>>23); 482116Sjkh if(m==0) { /* subnormal x */ 492116Sjkh for(i=0;(ix&0x00800000)==0;i++) ix<<=1; 502116Sjkh m -= i-1; 512116Sjkh } 522116Sjkh m -= 127; /* unbias exponent */ 532116Sjkh ix = (ix&0x007fffff)|0x00800000; 542116Sjkh if(m&1) /* odd m, double x to make it even */ 552116Sjkh ix += ix; 562116Sjkh m >>= 1; /* m = [m/2] */ 572116Sjkh 582116Sjkh /* generate sqrt(x) bit by bit */ 592116Sjkh ix += ix; 602116Sjkh q = s = 0; /* q = sqrt(x) */ 612116Sjkh r = 0x01000000; /* r = moving bit from right to left */ 622116Sjkh 632116Sjkh while(r!=0) { 648870Srgrimes t = s+r; 658870Srgrimes if(t<=ix) { 668870Srgrimes s = t+r; 678870Srgrimes ix -= t; 688870Srgrimes q += r; 698870Srgrimes } 702116Sjkh ix += ix; 712116Sjkh r>>=1; 722116Sjkh } 732116Sjkh 742116Sjkh /* use floating add to find out rounding direction */ 752116Sjkh if(ix!=0) { 762116Sjkh z = one-tiny; /* trigger inexact flag */ 772116Sjkh if (z>=one) { 782116Sjkh z = one+tiny; 792116Sjkh if (z>one) 802116Sjkh q += 2; 812116Sjkh else 822116Sjkh q += (q&1); 832116Sjkh } 842116Sjkh } 852116Sjkh ix = (q>>1)+0x3f000000; 862116Sjkh ix += (m <<23); 872116Sjkh SET_FLOAT_WORD(z,ix); 882116Sjkh return z; 892116Sjkh} 90