8000 Sumprod(): Update citation. Reorder functions. Add final twosum() call. Improve comments. by rhettinger · Pull Request #101249 · python/cpython · GitHub
[go: up one dir, main page]
Skip to content

Sumprod(): Update citation. Reorder functions. Add final twosum() call. Improve comments. #101249

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 8 commits into from
Jan 22, 2023
Merged

Sumprod(): Update citation. Reorder functions. Add final twosum() call. Improve comments. #101249

Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
2aa9a0e
Update citation. Reorder functions. Add final twosum() call.
rhettinger Jan 22, 2023
6ec01c2
The ASCII art was more distracting than helpful.
rhettinger Jan 22, 2023
12018cd
Segregate double length API from triple length API.
rhettinger Jan 22, 2023
7b5cf0c
Make comment match presentation in the paper.
rhettinger Jan 22, 2023
d271af1
Just two sources cover all the code.
rhettinger Jan 22, 2023
52c940b
Fix typo
rhettinger Jan 22, 2023
1974f59
Remove superfluous character
rhettinger Jan 22, 2023
5dc9099
Remove superfluous character
rhettinger Jan 22, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 21 additions & 38 deletions Modules/mathmodule.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2833,34 +2833,20 @@ long_add_would_overflow(long a, long b)

/*
Double and triple length extended precision floating point arithmetic
based on ideas from three sources:

Improved Kahan–Babuška algorithm by Arnold Neumaier
https://www.mat.univie.ac.at/~neum/scan/01.pdf
based on:

A Floating-Point Technique for Extending the Available Precision
by T. J. Dekker
https://csclub.uwaterloo.ca/~pbarfuss/dekker1971.pdf

Ultimately Fast Accurate Summation by Siegfried M. Rump
https://www.tuhh.de/ti3/paper/rump/Ru08b.pdf

Double length functions:
* dl_split() exact split of a C double into two half precision components.
* dl_mul() exact multiplication of two C doubles.

Triple length functions and constant:
* tl_zero is a triple length zero for starting or resetting an accumulation.
* tl_add() compensated addition of a C double to a triple length number.
* tl_fma() performs a triple length fused-multiply-add.
* tl_to_d() converts from triple length number back to a C double.
Accurate Sum and Dot Product
by Takeshi Ogita, Siegfried M. Rump, and Shin’Ichi Oishi
https://doi.org/10.1137/030601818
https://www.tuhh.de/ti3/paper/rump/OgRuOi05.pdf

*/

typedef struct{ double hi; double lo; } DoubleLength;
typedef struct{ double hi; double lo; double tiny; } TripleLength;

static const TripleLength tl_zero = {0.0, 0.0, 0.0};

static inline DoubleLength
twosum(double a, double b)
Expand All @@ -2874,25 +2860,9 @@ twosum(double a, double b)
return (DoubleLength) {s, t};
}

static inline TripleLength
tl_add(TripleLength total, double x)
{
/* Input: x total.hi total.lo total.tiny
|--- twosum ---|
s.hi s.lo
|--- twosum ---|
t.hi t.lo
|--- single sum ---|
Output: s.hi t.hi tiny
*/
DoubleLength s = twosum(x, total.hi);
DoubleLength t = twosum(s.lo, total.lo);
return (TripleLength) {s.hi, t.hi, t.lo + total.tiny};
}

static inline DoubleLength
dl_split(double x) {
double t = x * 134217729.0; /* Veltkamp constant = float(0x8000001) */
double t = x * 134217729.0; // Veltkamp constant = 2.0 ** 27 + 1
double hi = t - (t - x);
double lo = x - hi;
return (DoubleLength) {hi, lo};
Expand All @@ -2911,6 +2881,18 @@ dl_mul(double x, double y)
return (DoubleLength) {z, zz};
}

typedef struct{ double hi; double lo; double tiny; } TripleLength;

static const TripleLength tl_zero = {0.0, 0.0, 0.0};

static inline TripleLength
tl_add(TripleLength total, double x)
{
DoubleLength s = twosum(x, total.hi);
DoubleLength t = twosum(s.lo, total.lo);
return (TripleLength) {s.hi, t.hi, t.lo + total.tiny};
}

static inline TripleLength
tl_fma(TripleLength total, double p, double q)
{
Expand All @@ -2922,7 +2904,8 @@ tl_fma(TripleLength total, double p, double q)
static inline double
tl_to_d(TripleLength total)
{
return total.tiny + total.lo + total.hi;
DoubleLength last = twosum(total.lo, total.hi);
return total.tiny + last.lo + last.hi;
}

/*[clinic input]
Expand Down Expand Up @@ -3039,7 +3022,7 @@ math_sumprod_impl(PyObject *module, PyObject *p, PyObject *q)
}

finalize_int_path:
// # We're finished, overflowed, or have a non-int
// We're finished, overflowed, or have a non-int
int_path_enabled = false;
if (int_total_in_use) {
term_i = PyLong_FromLong(int_total);
Expand Down
0