e_sqrtf.c revision 8870
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 178870Srgrimesstatic char rcsid[] = "$Id: e_sqrtf.c,v 1.1.1.1 1994/08/19 09:39:57 jkh Exp $"; 182116Sjkh#endif 192116Sjkh 202116Sjkh#include "math.h" 212116Sjkh#include "math_private.h" 222116Sjkh 232116Sjkh#ifdef __STDC__ 242116Sjkhstatic const float one = 1.0, tiny=1.0e-30; 252116Sjkh#else 262116Sjkhstatic float one = 1.0, tiny=1.0e-30; 272116Sjkh#endif 282116Sjkh 292116Sjkh#ifdef __STDC__ 302116Sjkh float __ieee754_sqrtf(float x) 312116Sjkh#else 322116Sjkh float __ieee754_sqrtf(x) 332116Sjkh float x; 342116Sjkh#endif 352116Sjkh{ 362116Sjkh float z; 378870Srgrimes int32_t sign = (int)0x80000000; 382116Sjkh int32_t ix,s,q,m,t,i; 392116Sjkh u_int32_t r; 402116Sjkh 412116Sjkh GET_FLOAT_WORD(ix,x); 422116Sjkh 432116Sjkh /* take care of Inf and NaN */ 448870Srgrimes if((ix&0x7f800000)==0x7f800000) { 452116Sjkh return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf 462116Sjkh sqrt(-inf)=sNaN */ 478870Srgrimes } 482116Sjkh /* take care of zero */ 492116Sjkh if(ix<=0) { 502116Sjkh if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */ 512116Sjkh else if(ix<0) 522116Sjkh return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ 532116Sjkh } 542116Sjkh /* normalize x */ 552116Sjkh m = (ix>>23); 562116Sjkh if(m==0) { /* subnormal x */ 572116Sjkh for(i=0;(ix&0x00800000)==0;i++) ix<<=1; 582116Sjkh m -= i-1; 592116Sjkh } 602116Sjkh m -= 127; /* unbias exponent */ 612116Sjkh ix = (ix&0x007fffff)|0x00800000; 622116Sjkh if(m&1) /* odd m, double x to make it even */ 632116Sjkh ix += ix; 642116Sjkh m >>= 1; /* m = [m/2] */ 652116Sjkh 662116Sjkh /* generate sqrt(x) bit by bit */ 672116Sjkh ix += ix; 682116Sjkh q = s = 0; /* q = sqrt(x) */ 692116Sjkh r = 0x01000000; /* r = moving bit from right to left */ 702116Sjkh 712116Sjkh while(r!=0) { 728870Srgrimes t = s+r; 738870Srgrimes if(t<=ix) { 748870Srgrimes s = t+r; 758870Srgrimes ix -= t; 768870Srgrimes q += r; 778870Srgrimes } 782116Sjkh ix += ix; 792116Sjkh r>>=1; 802116Sjkh } 812116Sjkh 822116Sjkh /* use floating add to find out rounding direction */ 832116Sjkh if(ix!=0) { 842116Sjkh z = one-tiny; /* trigger inexact flag */ 852116Sjkh if (z>=one) { 862116Sjkh z = one+tiny; 872116Sjkh if (z>one) 882116Sjkh q += 2; 892116Sjkh else 902116Sjkh q += (q&1); 912116Sjkh } 922116Sjkh } 932116Sjkh ix = (q>>1)+0x3f000000; 942116Sjkh ix += (m <<23); 952116Sjkh SET_FLOAT_WORD(z,ix); 962116Sjkh return z; 972116Sjkh} 98