Lines Matching +full:ulp +full:- +full:0
27 * ------------------------------------------
29 * ------------------------------------------
38 * i 0
40 * s = 2*q , and y = 2 * ( y - q ). (1)
46 * -(i+1) 2
49 * -(i+1)
55 * -(i+1)
68 * -i -(i+1)
69 * s = s + 2 , y = y - s - 2 (5)
74 * it does not necessary to do a full (53-bit) comparison
79 * result is exact, bigger than 1/2ulp, or less than 1/2ulp
80 * (it will never equal to 1/2ulp).
82 * huge + tiny is equal to huge, and whether huge - tiny is
86 * sqrt(+-0) = +-0 ... exact
88 * sqrt(-ve) = NaN ... with invalid signal
92 *---------------
95 static const double one = 1.0, tiny=1.0e-300;
101 int32_t sign = (int)0x80000000; in sqrt()
108 if((ix0&0x7ff00000)==0x7ff00000) { in sqrt()
110 sqrt(-inf)=sNaN */ in sqrt()
113 if(ix0<=0) { in sqrt()
114 if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */ in sqrt()
115 else if(ix0<0) in sqrt()
116 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ in sqrt()
120 if(m==0) { /* subnormal x */ in sqrt()
121 while(ix0==0) { in sqrt()
122 m -= 21; in sqrt()
125 for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1; in sqrt()
126 m -= i-1; in sqrt()
127 ix0 |= (ix1>>(32-i)); in sqrt()
130 m -= 1023; /* unbias exponent */ in sqrt()
131 ix0 = (ix0&0x000fffff)|0x00100000; in sqrt()
141 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */ in sqrt()
142 r = 0x00200000; /* r = moving bit from right to left */ in sqrt()
144 while(r!=0) { in sqrt()
148 ix0 -= t; in sqrt()
157 while(r!=0) { in sqrt()
162 if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1; in sqrt()
163 ix0 -= t; in sqrt()
164 if (ix1 < t1) ix0 -= 1; in sqrt()
165 ix1 -= t1; in sqrt()
174 if((ix0|ix1)!=0) { in sqrt()
175 z = one-tiny; /* trigger inexact flag */ in sqrt()
178 if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;} in sqrt()
180 if (q1==(u_int32_t)0xfffffffe) q+=1; in sqrt()
186 ix0 = (q>>1)+0x3fe00000; in sqrt()
200 Other methods (use floating-point arithmetic)
201 -------------
214 standard (IEEE 754-1985). The ability to perform shift, add,
215 subtract and logical AND operations upon 32-bit words is needed
222 Let x0 and x1 be the leading and the trailing 32-bit words of
226 ------------------------------------------------------
228 ------------------------------------------------------
232 ------------------------ ------------------------
234 ------------------------ ------------------------
237 as integers), we obtain an 8-bit approximation of sqrt(x) as
240 k := (x0>>1) + 0x1ff80000;
241 y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits
242 Here k is a 32-bit integer and T1[] is an integer array containing
244 leading 32-bit word is y0, the value of its trailing word is 0)
245 approximates sqrt(x) to almost 8-bit.
249 0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,
257 sqrt(x) to within 1 ulp (Unit in the Last Place):
261 y := y-(y-x/y)/2 ... within 1 ulp
265 Another way to improve y to within 1 ulp is:
268 y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
271 (x-y )*y
272 y := y + 2* ---------- ...within 1 ulp
290 use the expression y+-ulp for the next representable floating
291 numbers (up and down) of y. Note that y+-ulp = either fixed
292 point y+-1, or multiply y by nextafter(1,+-inf) in chopped
296 R := RZ; ... set rounding mode to round-toward-zero
304 z := z - ulp; ... special rounding
308 If (r=RN) then z=z+ulp ... rounded-to-nearest
309 If (r=RP) then { ... round-toward-+inf
310 y = y+ulp; z=z+ulp;
313 y0:=y0-0x00100000; ... y := y/2 is correctly rounded.
320 Square root of +inf, +-0, or NaN is itself;
328 Let x0 and x1 be the leading and the trailing 32-bit words of
331 we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
333 k := 0x5fe80000 - (x0>>1);
334 y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits
336 Here k is a 32-bit integer and T2[] is an integer array
338 value of y (y's leading 32-bit word is y0, the value of
340 to almost 7.8-bit.
344 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
345 0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
346 0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
347 0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
348 0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
349 0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
350 0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
351 0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
357 to about 1 ulp. To be exact, we will have
358 -1ulp < sqrt(x)-z<1.0625ulp.
360 ... set rounding mode to Round-to-nearest
361 y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)
362 y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
365 z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)
367 Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
370 -1 ulp < sqrt(x) - z < 1.0625 ulp
371 instead of |sqrt(x)-z|<1.03125ulp.
379 use the expression y+-ulp for the next representable floating
380 numbers (up and down) of y. Note that y+-ulp = either fixed
381 point y+-1, or multiply y by nextafter(1,+-inf) in chopped
384 R := RZ; ... set rounding mode to round-toward-zero
386 case RN: ... round-to-nearest
387 if(x<= z*(z-ulp)...chopped) z = z - ulp; else
388 if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
390 case RZ:case RM: ... round-to-zero or round-to--inf
391 R:=RP; ... reset rounding mod to round-to-+inf
392 if(x<z*z ... rounded up) z = z - ulp; else
393 if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
395 case RP: ... round-to-+inf
396 if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
397 if(x>z*z ...chopped) z = z+ulp;
411 If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
414 j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2
415 k := z1 >> 26; ... get z's 25-th and 26-th
431 Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
434 --------------------
436 --------------------
437 bit 31 bit 0
439 Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
442 -------------------------------------------------
443 bit 27,26 of z1 bit 1,0 of x1 logb(x)
444 -------------------------------------------------
450 -------------------------------------------------