Programming Praxis

5 Responses to “Reciprocal Palindromes”

1. matthew said

   November 30, 2021 at 10:02 AM

   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

   November 30, 2021 at 10:20 AM

   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

   December 1, 2021 at 5:28 PM

   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

   December 2, 2021 at 8:55 AM

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

   April 18, 2022 at 9:29 PM

   https://pastebin.com/fpniQL4u