1b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project/* @(#)e_sqrt.c 1.3 95/01/18 */ 2b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project/* 3b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * ==================================================== 4b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 6b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Developed at SunSoft, a Sun Microsystems, Inc. business. 7b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Permission to use, copy, modify, and distribute this 8b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * software is freely granted, provided that this notice 9b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * is preserved. 10b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * ==================================================== 11b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project */ 12b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 13b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project/* __ieee754_sqrt(x) 14b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Return correctly rounded sqrt. 15b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * ------------------------------------------ 16b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * | Use the hardware sqrt if you have one | 17b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * ------------------------------------------ 18b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Method: 19b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Bit by bit method using integer arithmetic. (Slow, but portable) 20b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 1. Normalization 21b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Scale x to y in [1,4) with even powers of 2: 22b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then 23b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * sqrt(x) = 2^k * ieee_sqrt(y) 24b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 2. Bit by bit computation 25b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Let q = ieee_sqrt(y) truncated to i bit after binary point (q = 1), 26b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i 0 27b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i+1 2 28b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * s = 2*q , and y = 2 * ( y - q ). (1) 29b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i i i i 30b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 31b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * To compute q from q , one checks whether 32b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i+1 i 33b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 34b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * -(i+1) 2 35b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * (q + 2 ) <= y. (2) 36b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i 37b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * -(i+1) 38b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * If (2) is false, then q = q ; otherwise q = q + 2 . 39b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i+1 i i+1 i 40b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 41b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * With some algebric manipulation, it is not difficult to see 42b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * that (2) is equivalent to 43b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * -(i+1) 44b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * s + 2 <= y (3) 45b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i i 46b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 47b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * The advantage of (3) is that s and y can be computed by 48b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i i 49b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * the following recurrence formula: 50b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * if (3) is false 51b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 52b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * s = s , y = y ; (4) 53b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i+1 i i+1 i 54b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 55b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * otherwise, 56b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * -i -(i+1) 57b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * s = s + 2 , y = y - s - 2 (5) 58b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * i+1 i i+1 i i 59b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 60b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * One may easily use induction to prove (4) and (5). 61b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Note. Since the left hand side of (3) contain only i+2 bits, 62b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * it does not necessary to do a full (53-bit) comparison 63b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * in (3). 64b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 3. Final rounding 65b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * After generating the 53 bits result, we compute one more bit. 66b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Together with the remainder, we can decide whether the 67b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * result is exact, bigger than 1/2ulp, or less than 1/2ulp 68b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * (it will never equal to 1/2ulp). 69b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * The rounding mode can be detected by checking whether 70b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * huge + tiny is equal to huge, and whether huge - tiny is 71b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * equal to huge for some floating point number "huge" and "tiny". 72b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 73b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Special cases: 74b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * sqrt(+-0) = +-0 ... exact 75b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * sqrt(inf) = inf 76b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * sqrt(-ve) = NaN ... with invalid signal 77b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * sqrt(NaN) = NaN ... with invalid signal for signaling NaN 78b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * 79b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project * Other methods : see the appended file at the end of the program below. 80b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project *--------------- 81b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project */ 82b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 83b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project#include "fdlibm.h" 84b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 85b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project#ifdef __STDC__ 86b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Projectstatic const double one = 1.0, tiny=1.0e-300; 87b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project#else 88b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Projectstatic double one = 1.0, tiny=1.0e-300; 89b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project#endif 90b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 91b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project#ifdef __STDC__ 92b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project double __ieee754_sqrt(double x) 93b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project#else 94b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project double __ieee754_sqrt(x) 95b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project double x; 96b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project#endif 97b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project{ 98b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project double z; 99b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project int sign = (int)0x80000000; 100b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project unsigned r,t1,s1,ix1,q1; 101b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project int ix0,s0,q,m,t,i; 102b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 103b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 = __HI(x); /* high word of x */ 104b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix1 = __LO(x); /* low word of x */ 105b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 106b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project /* take care of Inf and NaN */ 107b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if((ix0&0x7ff00000)==0x7ff00000) { 108b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project return x*x+x; /* ieee_sqrt(NaN)=NaN, ieee_sqrt(+inf)=+inf 109b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ieee_sqrt(-inf)=sNaN */ 110b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 111b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project /* take care of zero */ 112b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(ix0<=0) { 113b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(((ix0&(~sign))|ix1)==0) return x;/* ieee_sqrt(+-0) = +-0 */ 114b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project else if(ix0<0) 115b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project return (x-x)/(x-x); /* ieee_sqrt(-ve) = sNaN */ 116b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 117b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project /* normalize x */ 118b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project m = (ix0>>20); 119b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(m==0) { /* subnormal x */ 120b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project while(ix0==0) { 121b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project m -= 21; 122b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 |= (ix1>>11); ix1 <<= 21; 123b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 124b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1; 125b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project m -= i-1; 126b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 |= (ix1>>(32-i)); 127b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix1 <<= i; 128b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 129b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project m -= 1023; /* unbias exponent */ 130b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 = (ix0&0x000fffff)|0x00100000; 131b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(m&1){ /* odd m, double x to make it even */ 132b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 += ix0 + ((ix1&sign)>>31); 133b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix1 += ix1; 134b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 135b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project m >>= 1; /* m = [m/2] */ 136b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 137b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project /* generate ieee_sqrt(x) bit by bit */ 138b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 += ix0 + ((ix1&sign)>>31); 139b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix1 += ix1; 140b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project q = q1 = s0 = s1 = 0; /* [q,q1] = ieee_sqrt(x) */ 141b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project r = 0x00200000; /* r = moving bit from right to left */ 142b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 143b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project while(r!=0) { 144b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project t = s0+r; 145b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(t<=ix0) { 146b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project s0 = t+r; 147b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 -= t; 148b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project q += r; 149b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 150b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 += ix0 + ((ix1&sign)>>31); 151b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix1 += ix1; 152b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project r>>=1; 153b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 154b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 155b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project r = sign; 156b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project while(r!=0) { 157b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project t1 = s1+r; 158b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project t = s0; 159b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 160b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project s1 = t1+r; 161b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1; 162b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 -= t; 163b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if (ix1 < t1) ix0 -= 1; 164b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix1 -= t1; 165b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project q1 += r; 166b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 167b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 += ix0 + ((ix1&sign)>>31); 168b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix1 += ix1; 169b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project r>>=1; 170b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 171b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 172b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project /* use floating add to find out rounding direction */ 173b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if((ix0|ix1)!=0) { 174b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project z = one-tiny; /* trigger inexact flag */ 175b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if (z>=one) { 176b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project z = one+tiny; 177b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if (q1==(unsigned)0xffffffff) { q1=0; q += 1;} 178b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project else if (z>one) { 179b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if (q1==(unsigned)0xfffffffe) q+=1; 180b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project q1+=2; 181b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } else 182b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project q1 += (q1&1); 183b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 184b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 185b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 = (q>>1)+0x3fe00000; 186b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix1 = q1>>1; 187b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if ((q&1)==1) ix1 |= sign; 188b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ix0 += (m <<20); 189b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project __HI(z) = ix0; 190b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project __LO(z) = ix1; 191b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project return z; 192b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project} 193b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 194b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project/* 195b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source ProjectOther methods (use floating-point arithmetic) 196b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project------------- 197b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project(This is a copy of a drafted paper by Prof W. Kahan 198b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Projectand K.C. Ng, written in May, 1986) 199b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 200b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Two algorithms are given here to implement ieee_sqrt(x) 201b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (IEEE double precision arithmetic) in software. 202b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Both supply ieee_sqrt(x) correctly rounded. The first algorithm (in 203b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Section A) uses newton iterations and involves four divisions. 204b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project The second one uses reciproot iterations to avoid division, but 205b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project requires more multiplications. Both algorithms need the ability 206b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project to chop results of arithmetic operations instead of round them, 207b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project and the INEXACT flag to indicate when an arithmetic operation 208b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project is executed exactly with no roundoff error, all part of the 209b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project standard (IEEE 754-1985). The ability to perform shift, add, 210b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project subtract and logical AND operations upon 32-bit words is needed 211b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project too, though not part of the standard. 212b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 213b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source ProjectA. ieee_sqrt(x) by Newton Iteration 214b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 215b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (1) Initial approximation 216b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 217b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Let x0 and x1 be the leading and the trailing 32-bit words of 218b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project a floating point number x (in IEEE double format) respectively 219b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 220b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 1 11 52 ...widths 221b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ------------------------------------------------------ 222b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project x: |s| e | f | 223b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ------------------------------------------------------ 224b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project msb lsb msb lsb ...order 225b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 226b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 227b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ------------------------ ------------------------ 228b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project x0: |s| e | f1 | x1: | f2 | 229b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ------------------------ ------------------------ 230b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 231b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project By performing shifts and subtracts on x0 and x1 (both regarded 232b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project as integers), we obtain an 8-bit approximation of ieee_sqrt(x) as 233b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project follows. 234b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 235b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project k := (x0>>1) + 0x1ff80000; 236b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y0 := k - T1[31&(k>>15)]. ... y ~ ieee_sqrt(x) to 8 bits 237b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Here k is a 32-bit integer and T1[] is an integer array containing 238b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project correction terms. Now magically the floating value of y (y's 239b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project leading 32-bit word is y0, the value of its trailing word is 0) 240b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project approximates ieee_sqrt(x) to almost 8-bit. 241b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 242b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Value of T1: 243b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project static int T1[32]= { 244b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 0, 1024, 3062, 5746, 9193, 13348, 18162, 23592, 245b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215, 246b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581, 247b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 16499, 12183, 8588, 5674, 3403, 1742, 661, 130,}; 248b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 249b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (2) Iterative refinement 250b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 251b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Apply Heron's rule three times to y, we have y approximates 252b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project sqrt(x) to within 1 ulp (Unit in the Last Place): 253b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 254b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y := (y+x/y)/2 ... almost 17 sig. bits 255b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y := (y+x/y)/2 ... almost 35 sig. bits 256b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y := y-(y-x/y)/2 ... within 1 ulp 257b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 258b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 259b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Remark 1. 260b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Another way to improve y to within 1 ulp is: 261b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 262b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y := (y+x/y) ... almost 17 sig. bits to 2*ieee_sqrt(x) 263b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y := y - 0x00100006 ... almost 18 sig. bits to ieee_sqrt(x) 264b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 265b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 2 266b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (x-y )*y 267b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y := y + 2* ---------- ...within 1 ulp 268b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 2 269b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 3y + x 270b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 271b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 272b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project This formula has one division fewer than the one above; however, 273b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project it requires more multiplications and additions. Also x must be 274b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project scaled in advance to avoid spurious overflow in evaluating the 275b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project expression 3y*y+x. Hence it is not recommended uless division 276b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project is slow. If division is very slow, then one should use the 277b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project reciproot algorithm given in section B. 278b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 279b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (3) Final adjustment 280b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 281b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project By twiddling y's last bit it is possible to force y to be 282b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project correctly rounded according to the prevailing rounding mode 283b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project as follows. Let r and i be copies of the rounding mode and 284b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project inexact flag before entering the square root program. Also we 285b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project use the expression y+-ulp for the next representable floating 286b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project numbers (up and down) of y. Note that y+-ulp = either fixed 287b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project point y+-1, or multiply y by ieee_nextafter(1,+-inf) in chopped 288b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project mode. 289b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 290b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project I := FALSE; ... reset INEXACT flag I 291b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project R := RZ; ... set rounding mode to round-toward-zero 292b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project z := x/y; ... chopped quotient, possibly inexact 293b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project If(not I) then { ... if the quotient is exact 294b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(z=y) { 295b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project I := i; ... restore inexact flag 296b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project R := r; ... restore rounded mode 297b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project return ieee_sqrt(x):=y. 298b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } else { 299b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project z := z - ulp; ... special rounding 300b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 301b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 302b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project i := TRUE; ... ieee_sqrt(x) is inexact 303b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project If (r=RN) then z=z+ulp ... rounded-to-nearest 304b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project If (r=RP) then { ... round-toward-+inf 305b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y = y+ulp; z=z+ulp; 306b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 307b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y := y+z; ... chopped sum 308b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y0:=y0-0x00100000; ... y := y/2 is correctly rounded. 309b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project I := i; ... restore inexact flag 310b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project R := r; ... restore rounded mode 311b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project return ieee_sqrt(x):=y. 312b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 313b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (4) Special cases 314b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 315b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Square root of +inf, +-0, or NaN is itself; 316b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Square root of a negative number is NaN with invalid signal. 317b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 318b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 319b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source ProjectB. ieee_sqrt(x) by Reciproot Iteration 320b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 321b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (1) Initial approximation 322b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 323b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Let x0 and x1 be the leading and the trailing 32-bit words of 324b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project a floating point number x (in IEEE double format) respectively 325b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (see section A). By performing shifs and subtracts on x0 and y0, 326b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project we obtain a 7.8-bit approximation of 1/ieee_sqrt(x) as follows. 327b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 328b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project k := 0x5fe80000 - (x0>>1); 329b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y0:= k - T2[63&(k>>14)]. ... y ~ 1/ieee_sqrt(x) to 7.8 bits 330b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 331b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Here k is a 32-bit integer and T2[] is an integer array 332b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project containing correction terms. Now magically the floating 333b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project value of y (y's leading 32-bit word is y0, the value of 334b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project its trailing word y1 is set to zero) approximates 1/ieee_sqrt(x) 335b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project to almost 7.8-bit. 336b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 337b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Value of T2: 338b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project static int T2[64]= { 339b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866, 340b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f, 341b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d, 342b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0, 343b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989, 344b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd, 345b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e, 346b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,}; 347b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 348b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (2) Iterative refinement 349b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 350b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Apply Reciproot iteration three times to y and multiply the 351b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project result by x to get an approximation z that matches ieee_sqrt(x) 352b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project to about 1 ulp. To be exact, we will have 353b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project -1ulp < ieee_sqrt(x)-z<1.0625ulp. 354b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 355b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ... set rounding mode to Round-to-nearest 356b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/ieee_sqrt(x) 357b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/ieee_sqrt(x) 358b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ... special arrangement for better accuracy 359b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project z := x*y ... 29 bits to ieee_sqrt(x), with z*y<1 360b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project z := z + 0.5*z*(1-z*y) ... about 1 ulp to ieee_sqrt(x) 361b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 362b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that 363b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (a) the term z*y in the final iteration is always less than 1; 364b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (b) the error in the final result is biased upward so that 365b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project -1 ulp < ieee_sqrt(x) - z < 1.0625 ulp 366b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project instead of |ieee_sqrt(x)-z|<1.03125ulp. 367b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 368b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (3) Final adjustment 369b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 370b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project By twiddling y's last bit it is possible to force y to be 371b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project correctly rounded according to the prevailing rounding mode 372b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project as follows. Let r and i be copies of the rounding mode and 373b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project inexact flag before entering the square root program. Also we 374b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project use the expression y+-ulp for the next representable floating 375b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project numbers (up and down) of y. Note that y+-ulp = either fixed 376b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project point y+-1, or multiply y by ieee_nextafter(1,+-inf) in chopped 377b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project mode. 378b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 379b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project R := RZ; ... set rounding mode to round-toward-zero 380b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project switch(r) { 381b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project case RN: ... round-to-nearest 382b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(x<= z*(z-ulp)...chopped) z = z - ulp; else 383b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp; 384b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project break; 385b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project case RZ:case RM: ... round-to-zero or round-to--inf 386b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project R:=RP; ... reset rounding mod to round-to-+inf 387b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(x<z*z ... rounded up) z = z - ulp; else 388b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp; 389b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project break; 390b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project case RP: ... round-to-+inf 391b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else 392b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project if(x>z*z ...chopped) z = z+ulp; 393b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project break; 394b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 395b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 396b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Remark 3. The above comparisons can be done in fixed point. For 397b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project example, to compare x and w=z*z chopped, it suffices to compare 398b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project x1 and w1 (the trailing parts of x and w), regarding them as 399b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project two's complement integers. 400b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 401b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ...Is z an exact square root? 402b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project To determine whether z is an exact square root of x, let z1 be the 403b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project trailing part of z, and also let x0 and x1 be the leading and 404b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project trailing parts of x. 405b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 406b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0 407b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project I := 1; ... Raise Inexact flag: z is not exact 408b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project else { 409b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project j := 1 - [(x0>>20)&1] ... j = ieee_logb(x) mod 2 410b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project k := z1 >> 26; ... get z's 25-th and 26-th 411b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project fraction bits 412b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project I := i or (k&j) or ((k&(j+j+1))!=(x1&3)); 413b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project } 414b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project R:= r ... restore rounded mode 415b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project return ieee_sqrt(x):=z. 416b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 417b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project If multiplication is cheaper then the foregoing red tape, the 418b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Inexact flag can be evaluated by 419b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 420b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project I := i; 421b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project I := (z*z!=x) or I. 422b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 423b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Note that z*z can overwrite I; this value must be sensed if it is 424b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project True. 425b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 426b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be 427b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project zero. 428b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 429b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project -------------------- 430b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project z1: | f2 | 431b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project -------------------- 432b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project bit 31 bit 0 433b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 434b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd 435b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project or even of ieee_logb(x) have the following relations: 436b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 437b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ------------------------------------------------- 438b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project bit 27,26 of z1 bit 1,0 of x1 logb(x) 439b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ------------------------------------------------- 440b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 00 00 odd and even 441b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 01 01 even 442b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 10 10 odd 443b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 10 00 even 444b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 11 01 even 445b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project ------------------------------------------------- 446b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 447b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project (4) Special cases (see (4) of Section A). 448b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 449b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project */ 450b07e1d9fd8d9e4e03698e0bd9bf77154c5390326The Android Open Source Project 451