Reciprocal Palindromes
November 30, 2021
Over on YouTube, Michael Penn proves that the sum of reciprocals of the palindromes converges:
But he doesn’t compute the value to which the sum converges.
Your task is to compute the infinite sum. When you are finished, you are welcome to read a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
You are adding up the reciprocals in reverse order of size, so you will converge on a sum that is unchanged when the next value is added in (seems to be about 3.370283225788489823 using doubles). We could generate the palindromes in reverse order (starting from somewhere suitably small) but a possibly nicer solution is to use the Kahan summation algorithm (https://en.wikipedia.org/wiki/Kahan_summation_algorithm). Alternatively, you can peek at https://oeis.org/A118031 to see the answer.
Some code using Kahan summation:
#include <vector> #include <cstdlib> #include <cstdio> void nextpal(std::vector<int> &s) { auto len = s.size(); for (auto n = (len+1)/2; n > 0; n--) { if (s[n-1] < 9) { s[len-n] = s[n-1] = s[n-1]+1; return; } else { s[len-n] = s[n-1] = 0; } } // wrap around to all zero, so extend s[0] = 1; s.push_back(1); } double todouble(std::vector<int> s) { double x = 0.0, r = 1.0; // Going from smaller to larger should be accurate // here, so no fancy summation needed. for (auto i = s.begin(); i != s.end(); i++) { x += (*i)*r; r *= 10.0; } return x; } int main(int argc, char *argv[]) { long count = std::strtol(argv[1],0,0); std::vector<int> s = { 1 }; double sum = 0, sum0 = 0; double c = 0; // Compensation for (long i = 1; ; i++) { if (i%count == 0) { std::printf("%20.18f %20.18f %ld\n", sum, sum0, i); } double y = 1/todouble(s); sum0 += y; // Add in to naive sum y -= c; // Subtract correction double t = sum + y; // Add to sum c = (t - sum) - y; // New compensation sum = t; // New sum nextpal(s); } }And some output, first column is compensated sum, second is “naive” sum, third is iterations. Convergence is slow so will have to run for a while to compare with OEIS value:
After running for nearly 24 hours and adding up the first trillion reciprocals, we’ve only got as far as the 11th decimal place, but looks like the summation algorithm is doing its job (https://oeis.org/A118031 gives 3.370283259497..):
Each decimal place takes around 10 times longer than the last, so I’ve stopped it there rather than running until Christmas. It would be fun to add something to the code that uses extrapolation to estimate where we are currently converging to.
It goes a little faster if the todouble() parameter is made a const reference.
https://pastebin.com/fpniQL4u