@@ -2419,9 +2419,9 @@ To avoid overflow/underflow and to achieve high accuracy giving results |
2419 | 2419 | that are almost always correctly rounded, four techniques are used: |
2420 | 2420 |
|
2421 | 2421 | * lossless scaling using a power-of-two scaling factor |
2422 | | -* accurate squaring using Veltkamp-Dekker splitting |
2423 | | -* compensated summation using a variant of the Neumaier algorithm |
2424 | | -* differential correction of the square root |
| 2422 | +* accurate squaring using Veltkamp-Dekker splitting [1] |
| 2423 | +* compensated summation using a variant of the Neumaier algorithm [2] |
| 2424 | +* differential correction of the square root [3] |
2425 | 2425 |
|
2426 | 2426 | The usual presentation of the Neumaier summation algorithm has an |
2427 | 2427 | expensive branch depending on which operand has the larger |
@@ -2456,7 +2456,11 @@ Given that csum >= 1.0, we have: |
2456 | 2456 | Since lo**2 is less than 1/2 ulp(csum), we have csum+lo*lo == csum. |
2457 | 2457 |
|
2458 | 2458 | To minimize loss of information during the accumulation of fractional |
2459 | | -values, each term has a separate accumulator. |
| 2459 | +values, each term has a separate accumulator. This also breaks up |
| 2460 | +sequential dependencies in the inner loop so the CPU can maximize |
| 2461 | +floating point throughput. [5] On a 2.6 GHz Haswell, adding one |
| 2462 | +dimension has an incremental cost of only 5ns -- for example when |
| 2463 | +moving from hypot(x,y) to hypot(x,y,z). |
2460 | 2464 |
|
2461 | 2465 | The square root differential correction is needed because a |
2462 | 2466 | correctly rounded square root of a correctly rounded sum of |
@@ -2466,20 +2470,32 @@ The differential correction starts with a value *x* that is |
2466 | 2470 | the difference between the square of *h*, the possibly inaccurately |
2467 | 2471 | rounded square root, and the accurately computed sum of squares. |
2468 | 2472 | The correction is the first order term of the Maclaurin series |
2469 | | -expansion of sqrt(h**2 + x) == h + x/(2*h) + O(x**2). |
| 2473 | +expansion of sqrt(h**2 + x) == h + x/(2*h) + O(x**2). [4] |
2470 | 2474 |
|
2471 | 2475 | Essentially, this differential correction is equivalent to one |
2472 | 2476 | refinement step in Newton's divide-and-average square root |
2473 | 2477 | algorithm, effectively doubling the number of accurate bits. |
2474 | 2478 | This technique is used in Dekker's SQRT2 algorithm and again in |
2475 | 2479 | Borges' ALGORITHM 4 and 5. |
2476 | 2480 |
|
| 2481 | +Without proof for all cases, hypot() cannot claim to be always |
| 2482 | +correctly rounded. However for n <= 1000, prior to the final addition |
| 2483 | +that rounds the overall result, the internal accuracy of "h" together |
| 2484 | +with its correction of "x / (2.0 * h)" is at least 100 bits. [6] |
| 2485 | +Also, hypot() was tested against a Decimal implementation with |
| 2486 | +prec=300. After 100 million trials, no incorrectly rounded examples |
| 2487 | +were found. In addition, perfect commutativity (all permutations are |
| 2488 | +exactly equal) was verified for 1 billion random inputs with n=5. [7] |
| 2489 | +
|
2477 | 2490 | References: |
2478 | 2491 |
|
2479 | 2492 | 1. Veltkamp-Dekker splitting: http://csclub.uwaterloo.ca/~pbarfuss/dekker1971.pdf |
2480 | 2493 | 2. Compensated summation: http://www.ti3.tu-harburg.de/paper/rump/Ru08b.pdf |
2481 | 2494 | 3. Square root differential correction: https://arxiv.org/pdf/1904.09481.pdf |
2482 | 2495 | 4. https://www.wolframalpha.com/input/?i=Maclaurin+series+sqrt%28h**2+%2B+x%29+at+x%3D0 |
| 2496 | +5. https://bugs.python.org/file49439/hypot.png |
| 2497 | +6. https://bugs.python.org/file49435/best_frac.py |
| 2498 | +7. https://bugs.python.org/file49448/test_hypot_commutativity.py |
2483 | 2499 |
|
2484 | 2500 | */ |
2485 | 2501 |
|
|