Skip to content

Commit b6b5a24

Browse files
RavenwaterRavenwater
committed
adding harmonic series convergence application
1 parent 5a869eb commit b6b5a24

File tree

2 files changed

+165
-1
lines changed

2 files changed

+165
-1
lines changed

‎applications/accuracy/mathematics/harmonic_series.cpp

+162
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
// harmonic_series.cpp: experiments with mixed-precision representations of the Harmonic Series
2+
//
3+
// Copyright (C) 2017 Stillwater Supercomputing, Inc.
4+
// SPDX-License-Identifier: MIT
5+
//
6+
// This file is part of the UNIVERSAL project, which is released under an MIT Open Source license.
7+
#include <universal/utility/directives.hpp>
8+
#include <universal/traits/arithmetic_traits.hpp>
9+
#include <universal/number/cfloat/cfloat.hpp>
10+
#include <universal/number/posit/posit.hpp>
11+
#include <universal/number/dd/dd.hpp>
12+
#include <numbers> // since C++20
13+
14+
namespace sw {
15+
namespace universal {
16+
17+
template<typename Scalar>
18+
Scalar ForwardHarmonicSeries(unsigned long long terms) {
19+
Scalar one{ 1 }, sum{ 0 };
20+
for (unsigned long long i = 1; i <= terms; ++i) {
21+
sum += one / Scalar(i);
22+
}
23+
return sum;
24+
}
25+
26+
template<typename Scalar>
27+
Scalar ReverseHarmonicSeries(unsigned long long terms) {
28+
Scalar one{ 1 }, sum{ 0 };
29+
for (unsigned long long i = terms; i > 0; --i) {
30+
sum += one / Scalar(i);
31+
}
32+
return sum;
33+
}
34+
35+
36+
template<typename Scalar>
37+
Scalar CompensatedHarmonicSeries(unsigned long long terms) {
38+
Scalar one{ 1 }, sum{ 0 };
39+
for (unsigned long long i = terms; i > 0; --i) {
40+
sum += one / Scalar(i);
41+
}
42+
return sum;
43+
}
44+
45+
template<typename Scalar>
46+
void HarmonicSeriesConvergence(unsigned orderOfMagnitude = 6) {
47+
using std::abs;
48+
49+
std::cout << "Harmonic Series Convergence for " << type_tag<Scalar>() << '\n';
50+
51+
std::vector<uint64_t> term;
52+
uint64_t scale = 100;
53+
term.push_back(scale);
54+
for (unsigned i = 3; i < orderOfMagnitude; ++i) {
55+
scale *= 10;
56+
term.push_back(scale);
57+
}
58+
59+
constexpr int max_digits10 = std::numeric_limits<Scalar>::max_digits10;
60+
constexpr unsigned WIDTH = max_digits10 + 8;
61+
auto oldPrec = std::cout.precision();
62+
std::cout << std::setprecision(max_digits10) << std::scientific;
63+
64+
std::cout
65+
<< std::setw(WIDTH) << "terms"
66+
<< std::setw(WIDTH) << "forward"
67+
<< std::setw(WIDTH) << "reverse"
68+
<< std::setw(WIDTH) << "difference"
69+
// << std::setw(WIDTH) << "compensated"
70+
// << std::setw(WIDTH) << "diff"
71+
<< '\n';
72+
for (uint64_t terms : term) {
73+
Scalar forwardSum = ForwardHarmonicSeries<Scalar>(terms);
74+
Scalar reverseSum = ReverseHarmonicSeries<Scalar>(terms);
75+
// Scalar compensatedSum = CompensatedHarmonicSeries<Scalar>(terms);
76+
Scalar diff_fr = abs(forwardSum - reverseSum);
77+
// Scalar diff_rc = abs(reverseSum - compensatedSum);
78+
79+
std::cout
80+
<< std::setw(WIDTH) << terms
81+
<< std::setw(WIDTH) << forwardSum
82+
<< std::setw(WIDTH) << reverseSum
83+
<< std::setw(WIDTH) << diff_fr
84+
// << std::setw(WIDTH) << compensatedSum
85+
// << std::setw(WIDTH) << diff_rc
86+
<< '\n';
87+
}
88+
89+
std::cout << std::setprecision(oldPrec) << std::defaultfloat;
90+
}
91+
92+
}
93+
}
94+
95+
int main()
96+
try {
97+
using namespace sw::universal;
98+
99+
HarmonicSeriesConvergence<float>();
100+
HarmonicSeriesConvergence<double>();
101+
HarmonicSeriesConvergence<dd>();
102+
103+
return EXIT_SUCCESS;
104+
}
105+
catch (char const* msg) {
106+
std::cerr << "Caught ad-hoc exception: " << msg << std::endl;
107+
return EXIT_FAILURE;
108+
}
109+
catch (const sw::universal::universal_arithmetic_exception& err) {
110+
std::cerr << "Caught unexpected universal arithmetic exception: " << err.what() << std::endl;
111+
return EXIT_FAILURE;
112+
}
113+
catch (const sw::universal::universal_internal_exception& err) {
114+
std::cerr << "Caught unexpected universal internal exception: " << err.what() << std::endl;
115+
return EXIT_FAILURE;
116+
}
117+
catch (std::runtime_error& err) {
118+
std::cerr << "Caught unexpected runtime error: " << err.what() << std::endl;
119+
return EXIT_FAILURE;
120+
}
121+
catch (...) {
122+
std::cerr << "Caught unknown exception" << std::endl;
123+
return EXIT_FAILURE;
124+
}
125+
126+
127+
/*
128+
Harmonic Series Convergence for float
129+
terms forward reverse difference
130+
100 5.187377930e+00 5.187376976e+00 9.536743164e-07
131+
1000 7.485478401e+00 7.485471725e+00 6.675720215e-06
132+
10000 9.787612915e+00 9.787604332e+00 8.583068848e-06
133+
100000 1.209085083e+01 1.209015274e+01 6.980895996e-04
134+
1000000 1.435735798e+01 1.439265156e+01 3.529357910e-02
135+
10000000 1.540368271e+01 1.668603134e+01 1.282348633e+00
136+
100000000 1.540368271e+01 1.880791855e+01 3.404235840e+00
137+
1000000000 1.540368271e+01 1.880791855e+01 3.404235840e+00
138+
10000000000 1.540368271e+01 1.880791855e+01 3.404235840e+00
139+
100000000000 1.540368271e+01 1.880791855e+01 3.404235840e+00
140+
Harmonic Series Convergence for double
141+
terms forward reverse difference
142+
100 5.18737751763962063e+00 5.18737751763962152e+00 8.88178419700125232e-16
143+
1000 7.48547086055034328e+00 7.48547086055034061e+00 2.66453525910037570e-15
144+
10000 9.78760603604434820e+00 9.78760603604438550e+00 3.73034936274052598e-14
145+
100000 1.20901461298633350e+01 1.20901461298634079e+01 7.28306304154102691e-14
146+
1000000 1.43927267228649889e+01 1.43927267228657723e+01 7.83373366175510455e-13
147+
10000000 1.66953113658572718e+01 1.66953113658599648e+01 2.69295696853077970e-12
148+
100000000 1.89978964138525548e+01 1.89978964138534465e+01 8.91731133378925733e-13
149+
1000000000 2.13004815023485499e+01 2.13004815023461482e+01 2.40163444686913863e-12
150+
10000000000 2.36030665949975003e+01 2.36030665948882685e+01 1.09231734768400202e-10
151+
100000000000 2.59056516865364301e+01 2.59056516878463441e+01 1.30991395508317510e-09
152+
Harmonic Series Convergence for double-double
153+
terms forward reverse difference
154+
100 5.1873775176396202608051176756582e+00 5.1873775176396202608051176756583e+00 9.8607613152626475676466070660348e-32
155+
1000 7.4854708605503449126565182043340e+00 7.4854708605503449126565182043338e+00 1.7256332301709633243381562365561e-31
156+
10000 9.7876060360443822641784779048557e+00 9.7876060360443822641784779048520e+00 3.7470892997998060757057106850932e-30
157+
100000 1.2090146129863427947363219363515e+01 1.2090146129863427947363219363505e+01 9.7621537021100210919701409953745e-30
158+
1000000 1.4392726722865723631381127493198e+01 1.4392726722865723631381127493196e+01 1.9721522630525295135293214132070e-30
159+
10000000 1.6695311365859851815399118939716e+01 1.6695311365859851815399118939541e+01 1.7562015902482775317978607184608e-28
160+
100000000 1.8997896413853898324417110394294e+01 1.8997896413853898324417110394212e+01 8.1745711303527348335790372577429e-29
161+
1000000000 2.1300481502347944016685101850059e+01 2.1300481502347944016685101848909e+01 1.1499619845859299593389473160410e-27
162+
*/

‎include/universal/number/dd/manipulators.hpp

+3-1
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,9 @@
1414
namespace sw { namespace universal {
1515

1616
// Generate a type tag for a doubledouble
17-
inline std::string type_tag(const dd& = {}) {
17+
template<typename DoubleDoubleType,
18+
std::enable_if_t< is_dd<DoubleDoubleType>, bool> = true >
19+
std::string type_tag(DoubleDoubleType = {}) {
1820
return std::string("double-double");
1921
}
2022

0 commit comments

Comments
 (0)