# Calculating the Riemann Zeta Function for x>2 real numbers

Here's the code:

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string.hpp>
#include <iostream>
#include <cmath>
#include <iomanip>
#include <limits>
#include <exception>

const int PRECISION = 100;
typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<PRECISION> > arbFloat;

// Prototypes
bool isStringValid(const std::string & str);
bool isNumberValid(const arbFloat & x);
inline std::string resizeArbtoString(const arbFloat & x);

int main(){
arbFloat zeta = 0.0;

std::cout << "\u03b6(s), s = ";
std::string inputStr;
std::getline(std::cin, inputStr);

if(!isStringValid(inputStr)) return 2;

arbFloat input = static_cast<arbFloat>(inputStr);
if(!isNumberValid(input)) return 3;

std::cout << std::setprecision(PRECISION);

int i = 1;
if(input == 0) {
zeta = -0.5;
} else {
std::string preComp, postComp;
std::cout << "Convergence:\n";
do {
preComp = resizeArbtoString(zeta);
zeta += pow(i, -input);
postComp = resizeArbtoString(zeta);
std::cout << '\t' + preComp + '\n';
i++;
} while(preComp != postComp);
}

std::cout << "\n\u03b6(" + inputStr + ") = " << zeta << "\nAfter " << i << " iterations";

return 0;
}

// Check input
bool isStringValid(const std::string & str){
// Check if string contains spaces
if(std::count(str.begin(), str.end(), ' ') > 0){
std::cout << "\nError: Input contains spaces\n";
return false;
}

// Check if string contains multiple .
if(std::count(str.begin(), str.end(), '.') > 1){
std::cout << "\nError: Input contains multiple decimal marks\n";
return false;
}

// Check if NaN or Out of bounds (due to parsing failure)
try{
boost::lexical_cast<arbFloat>(str);
} catch(std::runtime_error){
std::cout << "\nError: Unable to parse (value too large or incorrect number type)\n";
return false;
} catch(...){
std::cout << "\nError: Input is NaN\n";
return false;
}

// Check if intentional NaN
if(boost::icontains(str, "nan")){
std::cout << "\nError: Intentional NaN\n";
return false;
}

return true;
}

// Check number
bool isNumberValid(const arbFloat & x){
// Range check
if(x == std::numeric_limits<arbFloat>::infinity()){
std::cout << "\nError: Out of bounds\n";
return false;
}

if(x == 1){
std::cout << "\nError: Complex Infinity\n";
return false;
}

if(x < 2){
std::cout << "\nError: Does not converge (unsupported)\n";
return false;
}

return true;
}

inline std::string resizeArbtoString(const arbFloat & x){
std::string resizedStr = static_cast<std::string>(x);
resizedStr.resize(PRECISION + 2);
return resizedStr;
}


# Explanation

This is used to calculate the Riemann Zeta Function for real numbers that are greater than 2 (other than zero). Smaller values will take much longer. This is a result of the algorithm I use for approximating the Riemann Zeta function, not due to implementation. Likewise, larger values are much faster. If a larger value (around 500+) is entered, the precision of 100 digits is exceeded (it returns 1, but there are still digits past all those zeroes).

I am using three Boost libraries: multi-precision, lexical cast, and string algorithms. First I initialize an arbitrary float with the precision of 100 (the value in the code will be manually changed occasionally). I am using function prototypes.

First I define the variable zeta to be zero. Then I get the user input (example output is below). I then check if it is a valid string that will convert easily to an arbFloat by using the isStringValid function. I use a string at first to catch things like multiple periods or spaces (and so on), which result in differing behavior.

I then convert it to an arbFloat, and check if that number is valid using the isNumberValid function. Next I set the output precision.

Inside of the if-else statement, I do the calculations. Inside of the else portion, I first convert the initial value of zeta to a string using the resizeArbtoString function. The string is being resized to PRECISION + 2 to account for the number before the decimal mark, and the decimal mark itself.

I then do the actual calculations. Next, the resulting value of zeta is converted to a string using the same technique as before. The do-while loop checks if the strings are the same before and after calculations (if they are, it would result in an endless loop anyways).

Finally the value of zeta is printed along with the number of iterations.

# Example output

ζ(s), s = 100
Convergence:
0
1
1.0000000000000000000000000000007888609052210118054117285652827862296732064351090230047702789306640624
1.0000000000000000000000000000007888609052210118073520537827654190672617124639136733859843956171621636
1.0000000000000000000000000000007888609052210118073520537827660413687894985780843877923897736295862226
1.0000000000000000000000000000007888609052210118073520537827660413687896253431444106153299232999067602
1.0000000000000000000000000000007888609052210118073520537827660413687896253431459412620374098062482047
1.0000000000000000000000000000007888609052210118073520537827660413687896253431459412623465788470572268
1.0000000000000000000000000000007888609052210118073520537827660413687896253431459412623465793379665733
1.0000000000000000000000000000007888609052210118073520537827660413687896253431459412623465793379703381
1.0000000000000000000000000000007888609052210118073520537827660413687896253431459412623465793379703382

ζ(100) = 1.000000000000000000000000000000788860905221011807352053782766041368789625343145941262346579337970338
After 12 iterations
RUN SUCCESSFUL (total time: 3s)


Am I doing something inefficiently? How can I improve this code?

• You might be interested in sciencedirect.com/science/article/pii/S0377042700003368 – vnp Oct 28 '16 at 7:22
• (My no-brainer: (doc)comment your code.) – greybeard Nov 2 '16 at 22:14
• @vnp I understand that there are some optimizations to the algorithm I am using to calculate the Riemann Zeta function, however these use higher-level calculus that is difficult to implement in C++. I suppose I am looking for optimizations to the actual program, not the mathematical equation. Thanks though! – esote Nov 3 '16 at 20:50
• (Per)Using boost/math/special_functions/zeta.hpp is out of this question? A pity I can't read math: P.Borwein's Efficient Algorithms for the Riemann Zeta Function – greybeard Nov 5 '16 at 21:53
• @greybeard Haha, no. Even though I am using Boost, I want to calculate it without using the Boost Zeta function. And P. Borwein's document doesn't provide anything usable in C++ (it uses integrals and Bernoulli numbers and such). – esote Nov 6 '16 at 0:32

Even though you noted that boost includes them in their headers, it's better to include them in your code too. That way, you know what features you're going to use. Include guards will make sure that there's no compile-time overhead.

Next, it's a little bit difficult to differentiate between the type arbFloat and your functions. Everything you've written uses camelCase. An IDE makes this simple, of course, but since all your names are available in a single file, it's a little bit overwhelming at first.

# Magic numbers

const int PRECISION = 100;


Here we see a great use of const to prevent a magic number. We don't have to guess what 100 means at another point in your code. We only see the identifier PRECISION.

if(!isStringValid(inputStr)) return 2;

// ...

if(!isNumberValid(input)) return 3;


I suggest you to use a similar approach here, or an enum:

enum errorCode {
INVALID_STRING = 2,
INVALID_NUMBER = 3
};


# Keep the errors on the error stream

When a program has a lot of output, you usually pipe the output into a file:

echo "25" | ./zeta-program > zeta.dat


However, since all your error messages are also on std::cout, one will not see them. Use std::cerr instead of std::cout for errors.

There are some small things of with your algorithm. First of all, you display a string in every iteration. That's slow:

$echo "25" | time ./zeta-program ... lots of output 0.00user 0.00system 1:06.21elapsed$ echo "25" | time ./zeta-program > /dev/null
0.00user 0.00system 0:11.85elapsed


So maybe tone it down a little bit. Display the string every nth iteration. However, even the use of a string to compare the numbers is slightly weird. Transforming a number to base10 is usually slow compared to other artihmetic operations, and it doesn't cover your original intend.

Instead, write a function that does exactly what you want: check whether your value is close to your previous:

bool difference_is_small(const arbFloat & x, const arbFloat & y){
using boost::multiprecision::pow;
using boost::multiprecision::abs;

static const arbFloat margin = pow(10, static_cast<arbFloat>(-PRECISION));
return abs(x - y) <= margin;
}


With this function, we end up with

std::cout << "Convergence:\n";

arbFloat previous = 1;

while(!difference_is_small(previous,zeta)) {
previous = zeta;
zeta += pow(i, -input);
i++;
}


which clearly states what we're checking: the difference between the previous computation, and the following one. Note that you can still use your string-based comparison in difference_is_small, if you want to, but at that point it makes your code easier to read (and the check takes only 80% of the time your string-based check took).

Since posting this question, I have found ways to optimize it on my own:

Regarding the #includes: It seems that there are excess includes, instead of having all of these:

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string.hpp>
#include <iostream>
#include <cmath>
#include <iomanip>
#include <limits>
#include <exception>


I really only need these:

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/algorithm/string.hpp>
#include <iostream>


I assume this is because the Boost headers have subsequent includes that make the rest redundant. The functionality is the same.

For the do-while loop, instead of having:

do {
preComp = resizeArbtoString(zeta);
zeta += pow(i, -input);
postComp = resizeArbtoString(zeta);
std::cout << '\t' + preComp + '\n';
i++;
} while(preComp != postComp);


This can be done with a for-loop, like so:

for(i = 1; i == 1 || preComp != postComp; i++){
preComp = resizeArbtoString(zeta);
zeta += pow(i, -input);
postComp = resizeArbtoString(zeta);
std::cout << '\t' + preComp + '\n';
}


It removes the extra line for i++, and keeps the preComp != postComp inherent in the loop. This does not use the techniques described in the other answers, but nonetheless.

PS: For negative even integers, the result is zero. I have since added this to the program. It isn't an optimization, but I figured since it could help future visitors, it should be mentioned.

• In the convergence do ... while, why you are computing twice the same value between two iterations ?

Call the first preComp = resizeArbtoString(zeta); before the do { and store the postComp in the next preComp.

• Don't forget to compute the } while (...); condition before changing the preComp :

   bool bContinue = true;

preComp = resizeArbtoString(zeta);
do {
zeta += pow(i, -input);
postComp = resizeArbtoString(zeta);
std::cout << '\t' + preComp + '\n';
// compute the while condition before
bContinue = (preComp != postComp);
// store the current postComp for the next iteration
preComp = postComp;
i++;
} while(bContinue);

• why [compute the same value in two successive] iterations? why compare arbFloats using their truncated string representations in the first place? (For positive values,) comparing the (cheapest) logarithms of the current approximation and the absolute difference to its predecessor should do. – greybeard Nov 4 '16 at 7:02
• I wouldn't say that this answer improves the code, if anything it needlessly complicates it. Also, @greybeard I truncate the strings because the arbFloat has more precision (almost always 9 more digits) than my variable PRECISION. By truncating the string, it maintains the set precision and also decreases run time. – esote Nov 5 '16 at 13:45
• @Anonymous: I'm fully aware that you use the number of matching digits as a stand-in for logarithms - intentionally or not. I'm not ready to assume this as fast as using, say, logarithm base 2. I doubt it works in edge cases - 1.000000002(0²³²) vs. 1.000000001(9²³²) (meaning digit repeated many times). – greybeard Nov 5 '16 at 17:54
• @greybeard Feel free to submit an answer, if you think it'll improve the code – esote Nov 5 '16 at 19:16

(Not touching direct use of the Dirichlet series or stopping summation on the first difference between successive partial sums that is "too small in relation".)
The difference any given partial sum and the preceding one is the term just added.
Given that the range is from about 1.6 to 1.0, I suggest to just define a constant lower bound for terms to consider: const arbFloat IGNORE("0.015e-"+PRECISION); (syntax?)

    bool trace = true;
// don't bother computing pow(1, -input)
arbFloat zeta = 1.0;
int i = 2;
if(input == 0) {
zeta = -0.5;
} else {
// compute zeta from its series definition
if (trace) std::cout << "Convergence:\n\t0\n\t1\n";
do {
const arbFloat term = pow(i, -input);
zeta += term;
if (trace) std::cout << '\t' + zeta + '\n';
i += 1;
} while (IGNORE < term);
}


(If this worked (don't have a C++-environment I dare to try boost in), this would "trace" partial sums to the precision explicitly set, presumably.)