68 lines
1.6 KiB
C++
68 lines
1.6 KiB
C++
#include <array>
|
|
#include <cmath>
|
|
#include <iomanip>
|
|
#include <iostream>
|
|
|
|
double ByVaccaSeries(int numTerms)
|
|
{
|
|
// this method is simple but converges slowly
|
|
// calculate gamma by:
|
|
// 1 * (1/2 - 1/3) +
|
|
// 2 * (1/4 - 1/5 + 1/6 - 1/7) +
|
|
// 3 * (1/8 - 1/9 + 1/10 - 1/11 + 1/12 - 1/13 + 1/14 - 1/15) +
|
|
// 4 * ( . . . ) +
|
|
// . . .
|
|
double gamma = 0;
|
|
size_t next = 4;
|
|
|
|
for(double numerator = 1; numerator < numTerms; ++numerator)
|
|
{
|
|
double delta = 0;
|
|
for(size_t denominator = next/2; denominator < next; denominator+=2)
|
|
{
|
|
// calculate terms two at a time
|
|
delta += 1.0/denominator - 1.0/(denominator + 1);
|
|
}
|
|
|
|
gamma += numerator * delta;
|
|
next *= 2;
|
|
}
|
|
return gamma;
|
|
}
|
|
|
|
// based on the C entry
|
|
double ByEulersMethod()
|
|
{
|
|
//Bernoulli numbers with even indices
|
|
const std::array<double, 8> B2 {1.0, 1.0/6, -1.0/30, 1.0/42, -1.0/30,
|
|
5.0/66, -691.0/2730, 7.0/6};
|
|
|
|
const int n = 10;
|
|
|
|
//n-th harmonic number
|
|
const double h = [] // immediately invoked lambda
|
|
{
|
|
double sum = 1;
|
|
for (int k = 2; k <= n; k++) { sum += 1.0 / k; }
|
|
return sum - log(n);
|
|
}();
|
|
|
|
//expansion C = -digamma(1)
|
|
double a = -1.0 / (2*n);
|
|
double r = 1;
|
|
for (int k = 1; k < ssize(B2); k++)
|
|
{
|
|
r *= n * n;
|
|
a += B2[k] / (2*k * r);
|
|
}
|
|
|
|
return h + a;
|
|
}
|
|
|
|
|
|
int main()
|
|
{
|
|
std::cout << std::setprecision(16) << "Vacca series: " << ByVaccaSeries(32);
|
|
std::cout << std::setprecision(16) << "\nEulers method: " << ByEulersMethod();
|
|
}
|