There was an error while loading. Please reload this page.
1 parent e6dcd37 commit 27de286Copy full SHA for 27de286
1 file changed
Modules/mathmodule.c
@@ -2447,6 +2447,14 @@ precision. When a value at or below 1.0 is correctly rounded, it
2447
never goes above 1.0. And when values at or below 1.0 are squared,
2448
they remain at or below 1.0, thus preserving the summation invariant.
2449
2450
+Another interesting assertion is that csum+lo*lo == csum. In the loop,
2451
+each scaled vector element has a magnitude less than 1.0. After the
2452
+Veltkamp split, *lo* has a maximum value of 2**-27. So the maximum
2453
+value of *lo* squared is 2**-54. The value of ulp(1.0)/2.0 is 2**-53.
2454
+Given that csum >= 1.0, we have:
2455
+ lo**2 <= 2**-54 < 2**-53 == 1/2*ulp(1.0) <= ulp(csum)/2
2456
+Since lo**2 is less than 1/2 ulp(csum), we have csum+lo*lo == csum.
2457
+
2458
The square root differential correction is needed because a
2459
correctly rounded square root of a correctly rounded sum of
2460
squares can still be off by as much as one ulp.
@@ -2519,11 +2527,8 @@ vector_norm(Py_ssize_t n, double *vec, double max, int found_nan)
2519
2527
csum += x;
2520
2528
frac += (oldcsum - csum) + x;
2521
2529
2522
- x = lo * lo;
2523
- assert(fabs(csum) >= fabs(x));
2524
- oldcsum = csum;
2525
- csum += x;
2526
- frac += (oldcsum - csum) + x;
2530
+ assert(csum + lo * lo == csum);
2531
+ frac += lo * lo;
2532
}
2533
h = sqrt(csum - 1.0 + frac);
2534
0 commit comments