IsoSpec 2.2.1
Loading...
Searching...
No Matches
summator.h
1/*
2 * Copyright (C) 2015-2020 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 <vector>
21#include <utility>
22
23namespace IsoSpec
24{
25
27{
28 // Shewchuk algorithm
29 std::vector<double> partials;
30 int maxpart;
31 public:
32 inline SSummator()
33 { maxpart = 0; }
34
35 inline SSummator(const SSummator& other) :
36 partials(other.partials),
37 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
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
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