Last active
October 12, 2023 11:07
-
-
Save jhale/3b8ddb2eacf89174b19ab562f8240f44 to your computer and use it in GitHub Desktop.
Issue with summation MAMATH
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| import scipy as sp | |
| np.set_printoptions(precision=15) | |
| def maclaurin_series(n): | |
| m2pi = -2 * np.pi | |
| ns = np.arange(0, n) | |
| xs = np.power(m2pi, ns) / sp.special.factorial(ns) | |
| return xs | |
| def rel_error(computed, exact): | |
| return np.abs(computed - exact) / np.abs(exact) | |
| series = np.float32(maclaurin_series(64)) | |
| # Compute sum in much more precise arithmetic | |
| float64_sum = np.sum(series, dtype=np.float64) | |
| single_sum_orig = np.float32(0.0) | |
| for x in series: | |
| single_sum_orig = np.float32(single_sum_orig + x) | |
| order = np.argsort(np.abs(series)) | |
| single_sum_dec = np.float32(0.0) | |
| for x in reversed(series[order]): | |
| single_sum_dec = np.float32(single_sum_dec + x) | |
| print(f"Relative error (original order): {rel_error(single_sum_orig, float64_sum):.6e}") | |
| print(f"Relative error (decreasing order): {rel_error(single_sum_dec, float64_sum):.6e}") |
Author
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I looked into this issue with my colleague Michal Habera. Here is a working solution.
The issue with the overflow when computing the terms of the sum in float32 is separate from the issue of round error when computing the sum.
The solution is to compute the terms of the sum in float64, perform the sum in float64 precision and then cast the result to float32, which is essentially what was proposed in the original paper by Higham.