11dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 21dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* @(#)e_sqrt.c 1.3 95/01/18 */ 31dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 41dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ==================================================== 51dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 61dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 71dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Developed at SunSoft, a Sun Microsystems, Inc. business. 81dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Permission to use, copy, modify, and distribute this 91dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * software is freely granted, provided that this notice 101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * is preserved. 111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ==================================================== 121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 14a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes#include <sys/cdefs.h> 15a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes__FBSDID("$FreeBSD$"); 161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* __ieee754_sqrt(x) 181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Return correctly rounded sqrt. 191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ------------------------------------------ 201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * | Use the hardware sqrt if you have one | 211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ------------------------------------------ 221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Method: 231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Bit by bit method using integer arithmetic. (Slow, but portable) 241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 1. Normalization 251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Scale x to y in [1,4) with even powers of 2: 261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then 271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * sqrt(x) = 2^k * sqrt(y) 281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 2. Bit by bit computation 291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Let q = sqrt(y) truncated to i bit after binary point (q = 1), 301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i 0 311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i+1 2 321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * s = 2*q , and y = 2 * ( y - q ). (1) 331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i i i i 341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * To compute q from q , one checks whether 361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i+1 i 371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * -(i+1) 2 391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * (q + 2 ) <= y. (2) 401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i 411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * -(i+1) 421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * If (2) is false, then q = q ; otherwise q = q + 2 . 431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i+1 i i+1 i 441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * With some algebric manipulation, it is not difficult to see 461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * that (2) is equivalent to 471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * -(i+1) 481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * s + 2 <= y (3) 491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i i 501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * The advantage of (3) is that s and y can be computed by 521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i i 531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * the following recurrence formula: 541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * if (3) is false 551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * s = s , y = y ; (4) 571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i+1 i i+1 i 581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * otherwise, 601dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * -i -(i+1) 611dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * s = s + 2 , y = y - s - 2 (5) 621dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * i+1 i i+1 i i 631dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 641dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * One may easily use induction to prove (4) and (5). 651dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Note. Since the left hand side of (3) contain only i+2 bits, 661dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * it does not necessary to do a full (53-bit) comparison 671dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * in (3). 681dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 3. Final rounding 691dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * After generating the 53 bits result, we compute one more bit. 701dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Together with the remainder, we can decide whether the 711dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * result is exact, bigger than 1/2ulp, or less than 1/2ulp 721dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * (it will never equal to 1/2ulp). 731dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * The rounding mode can be detected by checking whether 741dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * huge + tiny is equal to huge, and whether huge - tiny is 751dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * equal to huge for some floating point number "huge" and "tiny". 761dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 771dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Special cases: 781dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * sqrt(+-0) = +-0 ... exact 791dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * sqrt(inf) = inf 801dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * sqrt(-ve) = NaN ... with invalid signal 811dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * sqrt(NaN) = NaN ... with invalid signal for signaling NaN 821dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 831dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Other methods : see the appended file at the end of the program below. 841dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project *--------------- 851dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 861dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 87a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes#include <float.h> 88a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes 891dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#include "math.h" 901dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#include "math_private.h" 911dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 921dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic const double one = 1.0, tiny=1.0e-300; 931dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 941dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectdouble 951dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project__ieee754_sqrt(double x) 961dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project{ 971dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double z; 981dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project int32_t sign = (int)0x80000000; 991dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project int32_t ix0,s0,q,m,t,i; 1001dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u_int32_t r,t1,s1,ix1,q1; 1011dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1021dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project EXTRACT_WORDS(ix0,ix1,x); 1031dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1041dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* take care of Inf and NaN */ 1051dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if((ix0&0x7ff00000)==0x7ff00000) { 1061dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf 1071dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project sqrt(-inf)=sNaN */ 1081dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1091dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* take care of zero */ 1101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(ix0<=0) { 1111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */ 1121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else if(ix0<0) 1131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ 1141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* normalize x */ 1161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project m = (ix0>>20); 1171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(m==0) { /* subnormal x */ 1181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project while(ix0==0) { 1191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project m -= 21; 1201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 |= (ix1>>11); ix1 <<= 21; 1211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1; 1231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project m -= i-1; 1241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 |= (ix1>>(32-i)); 1251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix1 <<= i; 1261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project m -= 1023; /* unbias exponent */ 1281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 = (ix0&0x000fffff)|0x00100000; 1291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(m&1){ /* odd m, double x to make it even */ 1301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 += ix0 + ((ix1&sign)>>31); 1311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix1 += ix1; 1321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project m >>= 1; /* m = [m/2] */ 1341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* generate sqrt(x) bit by bit */ 1361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 += ix0 + ((ix1&sign)>>31); 1371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix1 += ix1; 1381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */ 1391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r = 0x00200000; /* r = moving bit from right to left */ 1401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project while(r!=0) { 1421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t = s0+r; 1431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(t<=ix0) { 1441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project s0 = t+r; 1451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 -= t; 1461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q += r; 1471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 += ix0 + ((ix1&sign)>>31); 1491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix1 += ix1; 1501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r>>=1; 1511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r = sign; 1541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project while(r!=0) { 1551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t1 = s1+r; 1561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t = s0; 1571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 1581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project s1 = t1+r; 1591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1; 1601dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 -= t; 1611dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (ix1 < t1) ix0 -= 1; 1621dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix1 -= t1; 1631dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q1 += r; 1641dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1651dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 += ix0 + ((ix1&sign)>>31); 1661dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix1 += ix1; 1671dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r>>=1; 1681dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1691dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1701dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* use floating add to find out rounding direction */ 1711dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if((ix0|ix1)!=0) { 1721dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z = one-tiny; /* trigger inexact flag */ 1731dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (z>=one) { 1741dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z = one+tiny; 1751dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;} 1761dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else if (z>one) { 1771dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (q1==(u_int32_t)0xfffffffe) q+=1; 1781dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q1+=2; 1791dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else 1801dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q1 += (q1&1); 1811dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1821dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1831dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 = (q>>1)+0x3fe00000; 1841dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix1 = q1>>1; 1851dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if ((q&1)==1) ix1 |= sign; 1861dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ix0 += (m <<20); 1871dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project INSERT_WORDS(z,ix0,ix1); 1881dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return z; 1891dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project} 1901dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 191a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes#if (LDBL_MANT_DIG == 53) 192a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes__weak_reference(sqrt, sqrtl); 193a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes#endif 194a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes 1951dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 1961dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source ProjectOther methods (use floating-point arithmetic) 1971dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project------------- 1981dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project(This is a copy of a drafted paper by Prof W. Kahan 1991dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectand K.C. Ng, written in May, 1986) 2001dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2011dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Two algorithms are given here to implement sqrt(x) 2021dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (IEEE double precision arithmetic) in software. 2031dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Both supply sqrt(x) correctly rounded. The first algorithm (in 2041dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Section A) uses newton iterations and involves four divisions. 2051dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project The second one uses reciproot iterations to avoid division, but 2061dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project requires more multiplications. Both algorithms need the ability 2071dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project to chop results of arithmetic operations instead of round them, 2081dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project and the INEXACT flag to indicate when an arithmetic operation 2091dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project is executed exactly with no roundoff error, all part of the 2101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project standard (IEEE 754-1985). The ability to perform shift, add, 2111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project subtract and logical AND operations upon 32-bit words is needed 2121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project too, though not part of the standard. 2131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source ProjectA. sqrt(x) by Newton Iteration 2151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (1) Initial approximation 2171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Let x0 and x1 be the leading and the trailing 32-bit words of 2191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project a floating point number x (in IEEE double format) respectively 2201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1 11 52 ...widths 2221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ------------------------------------------------------ 2231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project x: |s| e | f | 2241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ------------------------------------------------------ 2251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project msb lsb msb lsb ...order 2261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ------------------------ ------------------------ 2291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project x0: |s| e | f1 | x1: | f2 | 2301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ------------------------ ------------------------ 2311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project By performing shifts and subtracts on x0 and x1 (both regarded 2331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project as integers), we obtain an 8-bit approximation of sqrt(x) as 2341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project follows. 2351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project k := (x0>>1) + 0x1ff80000; 2371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits 2381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Here k is a 32-bit integer and T1[] is an integer array containing 2391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project correction terms. Now magically the floating value of y (y's 2401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project leading 32-bit word is y0, the value of its trailing word is 0) 2411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project approximates sqrt(x) to almost 8-bit. 2421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Value of T1: 2441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project static int T1[32]= { 2451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 0, 1024, 3062, 5746, 9193, 13348, 18162, 23592, 2461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215, 2471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581, 2481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 16499, 12183, 8588, 5674, 3403, 1742, 661, 130,}; 2491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (2) Iterative refinement 2511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Apply Heron's rule three times to y, we have y approximates 2531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project sqrt(x) to within 1 ulp (Unit in the Last Place): 2541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y := (y+x/y)/2 ... almost 17 sig. bits 2561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y := (y+x/y)/2 ... almost 35 sig. bits 2571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y := y-(y-x/y)/2 ... within 1 ulp 2581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2601dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Remark 1. 2611dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Another way to improve y to within 1 ulp is: 2621dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2631dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x) 2641dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x) 2651dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2661dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2 2671dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (x-y )*y 2681dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y := y + 2* ---------- ...within 1 ulp 2691dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2 2701dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3y + x 2711dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2721dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2731dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project This formula has one division fewer than the one above; however, 2741dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project it requires more multiplications and additions. Also x must be 2751dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project scaled in advance to avoid spurious overflow in evaluating the 2761dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project expression 3y*y+x. Hence it is not recommended uless division 2771dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project is slow. If division is very slow, then one should use the 2781dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project reciproot algorithm given in section B. 2791dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2801dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (3) Final adjustment 2811dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2821dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project By twiddling y's last bit it is possible to force y to be 2831dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project correctly rounded according to the prevailing rounding mode 2841dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project as follows. Let r and i be copies of the rounding mode and 2851dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project inexact flag before entering the square root program. Also we 2861dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project use the expression y+-ulp for the next representable floating 2871dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project numbers (up and down) of y. Note that y+-ulp = either fixed 2881dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project point y+-1, or multiply y by nextafter(1,+-inf) in chopped 2891dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project mode. 2901dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2911dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project I := FALSE; ... reset INEXACT flag I 2921dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project R := RZ; ... set rounding mode to round-toward-zero 2931dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z := x/y; ... chopped quotient, possibly inexact 2941dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project If(not I) then { ... if the quotient is exact 2951dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(z=y) { 2961dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project I := i; ... restore inexact flag 2971dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project R := r; ... restore rounded mode 2981dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return sqrt(x):=y. 2991dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { 3001dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z := z - ulp; ... special rounding 3011dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 3021dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 3031dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project i := TRUE; ... sqrt(x) is inexact 3041dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project If (r=RN) then z=z+ulp ... rounded-to-nearest 3051dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project If (r=RP) then { ... round-toward-+inf 3061dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = y+ulp; z=z+ulp; 3071dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 3081dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y := y+z; ... chopped sum 3091dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y0:=y0-0x00100000; ... y := y/2 is correctly rounded. 3101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project I := i; ... restore inexact flag 3111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project R := r; ... restore rounded mode 3121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return sqrt(x):=y. 3131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (4) Special cases 3151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Square root of +inf, +-0, or NaN is itself; 3171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Square root of a negative number is NaN with invalid signal. 3181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source ProjectB. sqrt(x) by Reciproot Iteration 3211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (1) Initial approximation 3231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Let x0 and x1 be the leading and the trailing 32-bit words of 3251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project a floating point number x (in IEEE double format) respectively 3261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (see section A). By performing shifs and subtracts on x0 and y0, 3271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project we obtain a 7.8-bit approximation of 1/sqrt(x) as follows. 3281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project k := 0x5fe80000 - (x0>>1); 3301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits 3311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Here k is a 32-bit integer and T2[] is an integer array 3331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project containing correction terms. Now magically the floating 3341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project value of y (y's leading 32-bit word is y0, the value of 3351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project its trailing word y1 is set to zero) approximates 1/sqrt(x) 3361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project to almost 7.8-bit. 3371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Value of T2: 3391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project static int T2[64]= { 3401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866, 3411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f, 3421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d, 3431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0, 3441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989, 3451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd, 3461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e, 3471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,}; 3481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (2) Iterative refinement 3501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Apply Reciproot iteration three times to y and multiply the 3521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project result by x to get an approximation z that matches sqrt(x) 3531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project to about 1 ulp. To be exact, we will have 3541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project -1ulp < sqrt(x)-z<1.0625ulp. 3551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ... set rounding mode to Round-to-nearest 3571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x) 3581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x) 3591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ... special arrangement for better accuracy 3601dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z := x*y ... 29 bits to sqrt(x), with z*y<1 3611dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x) 3621dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3631dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that 3641dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (a) the term z*y in the final iteration is always less than 1; 3651dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (b) the error in the final result is biased upward so that 3661dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project -1 ulp < sqrt(x) - z < 1.0625 ulp 3671dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project instead of |sqrt(x)-z|<1.03125ulp. 3681dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3691dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (3) Final adjustment 3701dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3711dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project By twiddling y's last bit it is possible to force y to be 3721dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project correctly rounded according to the prevailing rounding mode 3731dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project as follows. Let r and i be copies of the rounding mode and 3741dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project inexact flag before entering the square root program. Also we 3751dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project use the expression y+-ulp for the next representable floating 3761dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project numbers (up and down) of y. Note that y+-ulp = either fixed 3771dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project point y+-1, or multiply y by nextafter(1,+-inf) in chopped 3781dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project mode. 3791dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3801dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project R := RZ; ... set rounding mode to round-toward-zero 3811dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project switch(r) { 3821dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project case RN: ... round-to-nearest 3831dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(x<= z*(z-ulp)...chopped) z = z - ulp; else 3841dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp; 3851dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project break; 3861dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project case RZ:case RM: ... round-to-zero or round-to--inf 3871dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project R:=RP; ... reset rounding mod to round-to-+inf 3881dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(x<z*z ... rounded up) z = z - ulp; else 3891dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp; 3901dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project break; 3911dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project case RP: ... round-to-+inf 3921dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else 3931dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(x>z*z ...chopped) z = z+ulp; 3941dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project break; 3951dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 3961dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 3971dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Remark 3. The above comparisons can be done in fixed point. For 3981dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project example, to compare x and w=z*z chopped, it suffices to compare 3991dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project x1 and w1 (the trailing parts of x and w), regarding them as 4001dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project two's complement integers. 4011dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4021dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ...Is z an exact square root? 4031dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project To determine whether z is an exact square root of x, let z1 be the 4041dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project trailing part of z, and also let x0 and x1 be the leading and 4051dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project trailing parts of x. 4061dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4071dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0 4081dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project I := 1; ... Raise Inexact flag: z is not exact 4091dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else { 4101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2 4111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project k := z1 >> 26; ... get z's 25-th and 26-th 4121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project fraction bits 4131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project I := i or (k&j) or ((k&(j+j+1))!=(x1&3)); 4141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 4151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project R:= r ... restore rounded mode 4161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return sqrt(x):=z. 4171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project If multiplication is cheaper then the foregoing red tape, the 4191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Inexact flag can be evaluated by 4201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project I := i; 4221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project I := (z*z!=x) or I. 4231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Note that z*z can overwrite I; this value must be sensed if it is 4251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project True. 4261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be 4281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project zero. 4291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project -------------------- 4311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z1: | f2 | 4321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project -------------------- 4331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project bit 31 bit 0 4341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd 4361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project or even of logb(x) have the following relations: 4371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ------------------------------------------------- 4391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project bit 27,26 of z1 bit 1,0 of x1 logb(x) 4401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ------------------------------------------------- 4411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 00 00 odd and even 4421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 01 01 even 4431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 10 10 odd 4441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 10 00 even 4451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 11 01 even 4461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ------------------------------------------------- 4471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project (4) Special cases (see (4) of Section A). 4491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 4501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 4511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 452