|
| 1 | +/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */ |
| 2 | +/* |
| 3 | + * ==================================================== |
| 4 | + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 5 | + * |
| 6 | + * Developed at SunSoft, a Sun Microsystems, Inc. business. |
| 7 | + * Permission to use, copy, modify, and distribute this |
| 8 | + * software is freely granted, provided that this notice |
| 9 | + * is preserved. |
| 10 | + * ==================================================== |
| 11 | + */ |
| 12 | +/* sqrt(x) |
| 13 | + * Return correctly rounded sqrt. |
| 14 | + * ------------------------------------------ |
| 15 | + * | Use the hardware sqrt if you have one | |
| 16 | + * ------------------------------------------ |
| 17 | + * Method: |
| 18 | + * Bit by bit method using integer arithmetic. (Slow, but portable) |
| 19 | + * 1. Normalization |
| 20 | + * Scale x to y in [1,4) with even powers of 2: |
| 21 | + * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then |
| 22 | + * sqrt(x) = 2^k * sqrt(y) |
| 23 | + * 2. Bit by bit computation |
| 24 | + * Let q = sqrt(y) truncated to i bit after binary point (q = 1), |
| 25 | + * i 0 |
| 26 | + * i+1 2 |
| 27 | + * s = 2*q , and y = 2 * ( y - q ). (1) |
| 28 | + * i i i i |
| 29 | + * |
| 30 | + * To compute q from q , one checks whether |
| 31 | + * i+1 i |
| 32 | + * |
| 33 | + * -(i+1) 2 |
| 34 | + * (q + 2 ) <= y. (2) |
| 35 | + * i |
| 36 | + * -(i+1) |
| 37 | + * If (2) is false, then q = q ; otherwise q = q + 2 . |
| 38 | + * i+1 i i+1 i |
| 39 | + * |
| 40 | + * With some algebric manipulation, it is not difficult to see |
| 41 | + * that (2) is equivalent to |
| 42 | + * -(i+1) |
| 43 | + * s + 2 <= y (3) |
| 44 | + * i i |
| 45 | + * |
| 46 | + * The advantage of (3) is that s and y can be computed by |
| 47 | + * i i |
| 48 | + * the following recurrence formula: |
| 49 | + * if (3) is false |
| 50 | + * |
| 51 | + * s = s , y = y ; (4) |
| 52 | + * i+1 i i+1 i |
| 53 | + * |
| 54 | + * otherwise, |
| 55 | + * -i -(i+1) |
| 56 | + * s = s + 2 , y = y - s - 2 (5) |
| 57 | + * i+1 i i+1 i i |
| 58 | + * |
| 59 | + * One may easily use induction to prove (4) and (5). |
| 60 | + * Note. Since the left hand side of (3) contain only i+2 bits, |
| 61 | + * it does not necessary to do a full (53-bit) comparison |
| 62 | + * in (3). |
| 63 | + * 3. Final rounding |
| 64 | + * After generating the 53 bits result, we compute one more bit. |
| 65 | + * Together with the remainder, we can decide whether the |
| 66 | + * result is exact, bigger than 1/2ulp, or less than 1/2ulp |
| 67 | + * (it will never equal to 1/2ulp). |
| 68 | + * The rounding mode can be detected by checking whether |
| 69 | + * huge + tiny is equal to huge, and whether huge - tiny is |
| 70 | + * equal to huge for some floating point number "huge" and "tiny". |
| 71 | + * |
| 72 | + * Special cases: |
| 73 | + * sqrt(+-0) = +-0 ... exact |
| 74 | + * sqrt(inf) = inf |
| 75 | + * sqrt(-ve) = NaN ... with invalid signal |
| 76 | + * sqrt(NaN) = NaN ... with invalid signal for signaling NaN |
| 77 | + */ |
| 78 | + |
1 | 79 | use core::f64;
|
2 | 80 |
|
3 | 81 | const TINY: f64 = 1.0e-300;
|
|
0 commit comments