Prime Palindromes

October 26, 2021

Today’s exercise is from a beginner programming web site:

A prime number, such as 127, has no factors other than itself and one. A palindrome, such as 121, is the same number when its digits are reversed. A prime palindrome, such as 131, is both prime and a palindrome. Find the 100th prime palindrome.

Your task is to write a program to find the 100th prime palindrome; answer the question as if you were explaining the task to a beginner programmer. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

16 Responses to “Prime Palindromes”

  1. kernelbob said

    No solution, I just feel compelled to mention Belphegor’s Prime, 1000000000000066600000000000001. A few years ago I played with numbers like this, and if I remember correctly, it’s not hard to construct primes of the form 1000…0something000…01 in any base. It’s a lucky coincidence that we attach numerological significance to 666, 13 (the number of consecutive zeros), and 31 (the total number of digits).

    https://en.wikipedia.org/wiki/Belphegor%27s_prime

  2. programmingpraxis said

    We discussed Belphegor’s Prime a long time ago.

  3. Richard A. O'Keefe said

    It is always worth trying the most obvious approach first. If it does not work, you may be able to tweak it to discover why it doesn’t, maybe get some insight into the problem space. For this problem, there are three issues: “prime”, “palindrome”, and “100th”. There is no obvious way of enumerating palindromic integers, but my personal library for Project Euler does include an #isPalindrome method for testing whether a natural number is a palindrome, and of course I have a method for enumerating primes. So the code goes like this:
    k := 100.
    Integer primesUpTo: 100000000 do: [:p |
    (p isPalindrome and: [(k := k – 1) isZero]) ifTrue: [
    Transcript print: p; cr.
    ^true]].
    Transcript nextPutAll: ‘Failed.’; cr.
    ^false

    The output is
    94049
    which is certainly a palindromic prime.
    In Haskell this would look much prettier:
    (filter isPalindrome primes) !! (100-1)
    where filter and !! are built in, primes can be found in many books, and only isPalindrome is tricky.

    Again, the obvious way to test whether an integer is a palindrome is to reverse its digits and see if you get the same number again. You could stop early, but let’s stick with an obvious approach.
    Integer
    methods for: ‘digits’
    reverseDigits
    |n r|
    n := self abs.
    r := 0.
    [0 < n] whileTrue: [
    r := r * 10 + (n \ 10).
    n := n // 10].
    ^self negative ifTrue: [r negated] ifFalse: [r]

    isPalindrome
      ^self reverseDigits = self
    

    Squeak Smalltalk has #primesUpTo:do: but not #isPalindrome. The same is true of Pharo.

  4. Richard A. O'Keefe said

    Drat. What do I do to preserve indentation?

  5. Ernest Gore said

    Note that we need only test palindromes with an odd number of digits since all palindromes with an even number of digits are divisible by 11

  6. matthew said

    Here’s an adaptation of some code from an earlier problem that directly generates palindromic sequences of digits 0-9 & checks for primality after a few syntactic checks (doesn’t end in even digit or 5, length is odd). To avoid special cases, start the enumeration from 101 and subtract 5 from target number (there are 5 palindromic primes less than 101). A rather naive prime checker – I tried the modular exponentiation check suggested by Praxis, using essentially the algorithm here, https://en.wikipedia.org/wiki/Modular_exponentiation#Pseudocode, but the squaring step overflowed long before the naive checker lost steam – I’m sure this can be worked around (“Shrage’s algorithm” maybe).

    #include <iostream>
    #include <vector>
    #include <cstdlib>
    
    #if !defined INT
    #define INT int
    #endif
    
    void nextpal(std::vector<int> &s) {
      int len = s.size();
      for (int 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);
    }
    
    INT toint(std::vector<int> s) {
      INT n = 0, r = 1;
      for (int i = 0; i != (int)s.size(); i++) {
        n += s[i]*r;
        r *= 10;
      }
      return n;
    }
    
    bool isprime(INT n) {
      for (INT i = 3; i*i <= n; i += 2) {
        if (n%i == 0) {
          return false;
        }
      }
      return true;
    }
    
    INT check(const std::vector<int> &s) {
      if (s.size()%2 == 0) return -1;
      if (s[0]%2 == 0 || s[0] == 5) return -1;
      INT n = toint(s);
      if (!isprime(n)) return -1;
      return n;
    }
    
    int main(int argc, char *argv[]) {
      int count = std::atoi(argv[1])-5;
      std::vector<int> s = { 1, 0, 1 };
      INT p = -1;
      while (count > 0) {
        p = check(s);
        if (p >= 0) count--;
        nextpal(s);
      }
      std::cout << p << "\n";
    }
    

    Using 32 bit ints constrains how far we can get (it’s nice that C compilers can now actually check for overflow etc), but 64 bit longs save the day:

    $ g++ -fsanitize=undefined -Wall -O3 palindromes.cpp -o palindromes
    $ ./palindromes 100
    94049
    $ ./palindromes 5000
    920787029
    $ ./palindromes 6000
    palindromes.cpp:28:7: runtime error: signed integer overflow: 1000000000 * 10 cannot be represented in type 'int'
    1448448409
    $ g++ -DINT=long -fsanitize=undefined -Wall -O3 palindromes.cpp -o palindromes
    $ ./palindromes 6000
    10050505001
    $ time ./palindromes 10000
    13649694631
    
    real	0m3.882s
    user	0m3.873s
    sys	0m0.008s
    
  7. matthew said

    BTW, Schrage’s method for doing modular exponentiation without overflow is covered here: https://programmingpraxis.com/2014/01/14/

  8. […] Came across this on the “Programming Praxis” blog […]

  9. A Mathematica version:

        Part[Select[Prime /@ Range[10000], PalindromeQ], 100] 
    

    Evaluates to:

    94049

  10. matthew said

    Actually, we can implement a suitable modular exponent function very easily, using the (non-standard) 128 bit integer types supported by most compilers & architectures:

    int64_t modmul(int64_t p, int64_t q, int64_t m) {
      return __int128(p)*q%m;
    }
    
    INT expm(INT p, INT q, INT m) {
      INT res = 1;
      while (true) {
        if (q&1) res = modmul(res,p,m);
        q >>= 1;
        if (q == 0) break;
        p = modmul(p,p,m);
      }
      return res;
    }
    
    bool isprime(INT n) {
      if (expm(2,n-1,n) != 1) return false;
      for (INT i = 3; i*i <= n; i += 2) {
        if (n%i == 0) return false;
      }
      return true;
    }
    

    This doesn’t speed things up much though, as we spend the majority of the time checking that the numbers that pass the expm test are indeed prime. However, a glance at https://oeis.org/A068445, “Palindromic pseudoprimes (base 2)” shows there aren’t going to be many false positives so we could skip the confirmation (1837381 is the only one in fact, since it doesn’t seem likely our program will get to 127665878878566721 or 1037998220228997301).

  11. matthew said

    This is the assembler output for that modmul function though (gcc 11.2, -O3, on http://godbolt.org):

    modmul(long, long, long):
            mov     rax, rdi
            mov     r8, rdx
            sub     rsp, 8
            imul    rsi
            mov     rcx, r8
            sar     rcx, 63
            mov     rsi, rdx
            mov     rdi, rax
            mov     rdx, r8
            call    __modti3
            add     rsp, 8
            ret
    

    It’s a little long winded & in fact calls out to a function for the actual mod operation. Let’s see if inline assembler can do better:

    inline int64_t mod128to64(int64_t high, int64_t low, int64_t divisor) {
      int64_t quotient, remainder;
      // 128 bit signed division, with 64 bit quotient and remainder
      __asm__ ("idivq %4" // Mnemonic, eg. "idivq rcx"
               : "=a" (quotient), "=d" (remainder) // rax, rdx are outputs
               : "a" (low), "d" (high), "rm" (divisor) // rax, rdx are dividend
               );
      return remainder;
    }
    
    int64_t modmul(int64_t p, int64_t q, int64_t m) {
      __int128 t = p;
      t *= q;
      return mod128to64(t>>64,t,m);
    }
    

    A Godbolt shows a more satisfactory result:

    modmul(long, long, long):
            mov     rax, rsi
            mov     r8, rdx
            imul    rdi
            idivq r8
            mov     rax, rdx
            ret
    

    but again the effect on runtime is fairly small as it turns out that nearly all the time is in that idiv instruction and the function call overhead is comparatively small.

  12. matthew said

    See https://stackoverflow.com/questions/32540740/intrinsics-for-128-multiplication-and-division for more details, and why there aren’t 128 division intrinsics (which would be a better choice if available). Also, my inline assembler should probably indicate the condition code register may be clobbered:

    #define IDIVQ(low,high,divisor,quot,rem) \
      __asm__ ("idivq %4" \
               : "=a" (quot), "=d" (rem) \
               : "a" (low), "d" (high), "rm" (divisor) \
               : "cc");
    
  13. matthew said

    And it would seem, if anyone cares, that 9223372002002733229 is the largest palindromic prime less than 2^63.

  14. Egil H said

    λ> isPrime n = not . any (== 0) . map (mod n) $ [2..(n-1)]
    λ> isPalindrome n = let n’ = show n in n’ == (reverse n’)
    λ> last . take 100 . filter isPrime . filter isPalindrome $ [2..]
    94049

  15. Aazad said

    Answer is 94949

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: