-
Notifications
You must be signed in to change notification settings - Fork 0
/
pairwisesum.cpp
129 lines (124 loc) · 2.76 KB
/
pairwisesum.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
/******************************************************/
/* */
/* pairwisesum.cpp - add many numbers */
/* */
/******************************************************/
/* Copyright 2015,2016,2019,2021 Pierre Abbat.
* Licensed under the Apache License, Version 2.0.
* This file is part of AGM.
*/
#include <cfloat>
#include <cstring>
#include <iostream>
#include <cmath>
#include <vector>
#include "pairwisesum.h"
using namespace std;
double pairwisesum(double *a,unsigned n)
/* Adds up n numbers at a pairwise, for less roundoff error.
* i^(i+1) is 1 if i is even, 3 if i%4==1, 7 if i%8==3, etc.
* i^(i+8) is 8 if i%16==0, 24 if i%32==16, etc.
*/
{
unsigned i,j,b;
double sums[32],sum=0;
for (i=0;i+7<n;i+=8) // Add up groups of 8, leaving 0 to 7 numbers at the end
{
b=i^(i+8);
if (b==8)
sums[3]=(((a[i]+a[i+1])+(a[i+2]+a[i+3]))+((a[i+4]+a[i+5])+(a[i+6]+a[i+7])));
else
{
sums[3]+=(((a[i]+a[i+1])+(a[i+2]+a[i+3]))+((a[i+4]+a[i+5])+(a[i+6]+a[i+7])));
for (j=4;b>>(j+1);j++)
sums[j]+=sums[j-1];
sums[j]=sums[j-1];
}
}
for (;i<n;i++) // Add up the last 0 to 7 numbers
{
b=i^(i+1);
if (b==1)
sums[0]=a[i];
else
{
sums[0]+=a[i];
for (j=1;b>>(j+1);j++)
sums[j]+=sums[j-1];
sums[j]=sums[j-1];
}
}
for (i=0;i<32;i++)
if ((n>>i)&1)
sum+=sums[i];
return sum;
}
long double pairwisesum(long double *a,unsigned n)
{
unsigned i,j,b;
long double sums[32],sum=0;
for (i=0;i+7<n;i+=8)
{
b=i^(i+8);
if (b==8)
sums[3]=(((a[i]+a[i+1])+(a[i+2]+a[i+3]))+((a[i+4]+a[i+5])+(a[i+6]+a[i+7])));
else
{
sums[3]+=(((a[i]+a[i+1])+(a[i+2]+a[i+3]))+((a[i+4]+a[i+5])+(a[i+6]+a[i+7])));
for (j=4;b>>(j+1);j++)
sums[j]+=sums[j-1];
sums[j]=sums[j-1];
}
}
for (;i<n;i++)
{
b=i^(i+1);
if (b==1)
sums[0]=a[i];
else
{
sums[0]+=a[i];
for (j=1;b>>(j+1);j++)
sums[j]+=sums[j-1];
sums[j]=sums[j-1];
}
}
for (i=0;i<32;i++)
if ((n>>i)&1)
sum+=sums[i];
return sum;
}
/* This is the original version of pairwisesum.
* This code is left here to show what the optimized version is trying
* to accomplish and as a reference for unit tests.
*/
#if 0
double pairwisesum(double *a,unsigned n)
// a is clobbered.
{
unsigned i,j;
if (n)
{
for (i=1;i<n;i*=2)
for (j=0;j+i<n;j+=2*i)
a[j]+=a[j+i];
return a[0];
}
else
return 0;
}
#endif
double pairwisesum(vector<double> &a)
{
if (a.size())
return pairwisesum(&a[0],a.size());
else
return 0;
}
long double pairwisesum(vector<long double> &a)
{
if (a.size())
return pairwisesum(&a[0],a.size());
else
return 0;
}