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