IsoSpec  1.95
summator.h
1 /*
2  * Copyright (C) 2015-2018 Mateusz Łącki and Michał Startek.
3  *
4  * This file is part of IsoSpec.
5  *
6  * IsoSpec is free software: you can redistribute it and/or modify
7  * it under the terms of the Simplified ("2-clause") BSD licence.
8  *
9  * IsoSpec is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
12  *
13  * You should have received a copy of the Simplified BSD Licence
14  * along with IsoSpec. If not, see <https://opensource.org/licenses/BSD-2-Clause>.
15  */
16 
17 #pragma once
18 
19 #include <cmath>
20 #include <atomic>
21 
22 namespace IsoSpec
23 {
24 
25 class SSummator
26 {
27  // Shewchuk algorithm
28  std::vector<double> partials;
29  int maxpart;
30 public:
31  inline SSummator()
32  { maxpart = 0; }
33 
34  inline SSummator(SSummator& other)
35  {
36  this->partials = other.partials;
37  this->maxpart = other.maxpart;
38  }
39  inline void add(double x)
40  {
41  unsigned int i=0;
42  for(int pidx=0; pidx<maxpart; pidx++)
43  {
44  double y = partials[pidx];
45  if(std::abs(x) < std::abs(y))
46  std::swap(x, y);
47  double hi = x+y;
48  double lo = y-(hi-x);
49  if(lo != 0.0)
50  {
51  partials[i] = lo;
52  i += 1;
53  }
54  x = hi;
55  }
56  while(partials.size() <= i)
57  partials.push_back(0.0);
58  partials[i] = x;
59  maxpart = i+1;
60  }
61  inline double get()
62  {
63  double ret = 0.0;
64  for(int i=0; i<maxpart; i++)
65  ret += partials[i];
66  return ret;
67  }
68 };
69 
70 
71 
72 
73 
74 
75 
76 class Summator{
77  // Kahan algorithm
78  double sum;
79  double c;
80 
81 public:
82  inline Summator()
83  { sum = 0.0; c = 0.0;}
84 
85  inline void add(double what)
86  {
87  double y = what - c;
88  double t = sum + y;
89  c = (t - sum) - y;
90  sum = t;
91  }
92 
93  inline double get()
94  {
95  return sum;
96  }
97 };
98 
99 class TSummator
100 {
101  // Trivial algorithm, for testing only
102  double sum;
103 public:
104  inline TSummator()
105  { sum = 0.0; }
106 
107  inline void add(double what)
108  {
109  sum += what;
110  }
111  inline double get()
112  {
113  return sum;
114  }
115 };
116 
117 
118 } // namespace IsoSpec
119