## Reciprocal Palindromes

### November 30, 2021

We wrote a palindrome generator in a previous exercise, which makes this easy:

(define-generator (palindromes)
(do ((k 0 (+ k 1))) ((= k 10))
(yield k))
(do ((i 1 (* i 10))) (#f)
(do ((j i (+ j 1))) ((= j (* 10 i)))
(let ((ds (digits j)))
(yield (undigits (append ds (reverse ds))))))
(do ((j i (+ j 1))) ((= j (* 10 i)))
(let ((ds (digits j)))
(do ((k 0 (+ k 1))) ((= k 10))
(yield (undigits (append ds (list k) (reverse ds)))))))))

(define (nth-reciprocal-palindrome-sum n)
(let ((pal (palindromes)))
(pal) ; discard 0
(do ((p (pal) (pal))
(sum 0 (+ sum (/ 1.0 p))))
((< n p) sum))))

> (nth-reciprocal-palindrome-sum 1000)
3.3190020006509102
> (nth-reciprocal-palindrome-sum 1000000)
3.3674689702500697
> (nth-reciprocal-palindrome-sum 1000000000)
3.3702320909393553
> (nth-reciprocal-palindrome-sum 1000000000000)
3.370280445226788
> (nth-reciprocal-palindrome-sum 100000000000000)
3.3702829780846546
> (nth-reciprocal-palindrome-sum 100000000000000000)
3.37028322578849

Okay, this is getting ridiculous.

Pages: 1 2

### 5 Responses to “Reciprocal Palindromes”

1. matthew said

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.

2. matthew said

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:

$g++ -Wall -O3 rpalindromes.cpp -o rpalindromes$ ./rpalindromes 100000000
3.370283197792765417 3.370283199089854076 100000000
3.370283231354666675 3.370283225788489823 200000000
3.370283238286138516 3.370283225788489823 300000000
3.370283242340789354 3.370283225788489823 400000000
3.370283245217610357 3.370283225788489823 500000000
3.370283247449045838 3.370283225788489823 600000000
3.370283249272261195 3.370283225788489823 700000000
3.370283250813768117 3.370283225788489823 800000000
3.370283252149082198 3.370283225788489823 900000000
3.370283253326912476 3.370283225788489823 1000000000
...

3. matthew said

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..):

3.370283259491198447 3.370283225788489823 999600000000
3.370283259491199335 3.370283225788489823 999700000000
3.370283259491200667 3.370283225788489823 999800000000
3.370283259491201555 3.370283225788489823 999900000000
3.370283259491202887 3.370283225788489823 1000000000000
^C


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.

4. matthew said

It goes a little faster if the todouble() parameter is made a const reference.