Prime Palindromes
October 26, 2021
I’m pretty sure the instructor of that beginner programmer had something like this in mind:
(define (palindrome? n)
(define (rev-digits n)
(let loop ((n n) (sum 0))
(if (zero? n) sum
(loop (quotient n 10)
(+ (* 10 sum) (modulo n 10))))))
(= n (rev-digits n)))
The digits of the number are stripped, from right to left, using the quotient and modulo operators, computing the sum of the digits multiplied by their respective powers of ten. At the end, the sum is compared to the original to determine if the number is a palindrome.
(define (prime? n)
(if (even? n) (= n 2)
(let loop ((d 3))
(if (< n (* d d)) #t
(if (zero? (modulo n d)) #f
(loop (+ d 2)))))))
Primality is determined by trial division, with a simple first check for divisibility by 2 followed by trial division by the odd numbers up to the square root of n.
(define (pp k)
(let loop ((k k) (n 2) (ns (list)))
(if (zero? k) (reverse ns)
(if (and (prime? n) (palindrome? n))
(loop (- k 1) (+ n 1) (cons n ns))
(loop k (+ n 1) ns)))))
Here we simply accumulate the list of prime palindromes, starting the test from 2.
> (pp 100) (2 3 5 7 11 101 131 151 181 191 313 353 373 383 727 757 787 797 919 929 10301 10501 10601 11311 11411 12421 12721 12821 13331 13831 13931 14341 14741 15451 15551 16061 16361 16561 16661 17471 17971 18181 18481 19391 19891 19991 30103 30203 30403 30703 30803 31013 31513 32323 32423 33533 34543 34843 35053 35153 35353 35753 36263 36563 37273 37573 38083 38183 38783 39293 70207 70507 70607 71317 71917 72227 72727 73037 73237 73637 74047 74747 75557 76367 76667 77377 77477 77977 78487 78787 78887 79397 79697 79997 90709 91019 93139 93239 93739 94049)
The program as written is very inefficient, primarily because almost all palindromes are composite. Any palindrome that ends in 2, 4, 5, 6, 8 or 0 is certainly not prime, except the trivial cases of 2 and 5. Further, any palindrome with an even number of digits is divisible by 11, so it is certainly not prime, except the trivial case of 11. Eliminating those possibilities gives a lot fewer cases to check. And, since so few palindromes are prime, you can simply check if 2**(n-1) not in (1, n-1) to quickly identify composites, reserving the full primality test for numbers that pass the simpler test. For instance, here is an abbreviated prime test (Henri Cohen, who knows a thing or two about prime numbers, calls numbers that pass this test “industrial grade primes.”)
(define (p? n) ; industrial-grade prime
(let ((x (expm 2 (- n 1) n)))
(or (= x 1) (= x (- n 1)))))
You can run the program at https://ideone.com/4dWXMV. If you want to know more about prime palindromes, look at http://oeis.org/A002385.
Programming Praxis 
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
We discussed Belphegor’s Prime a long time ago.
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]
Squeak Smalltalk has #primesUpTo:do: but not #isPalindrome. The same is true of Pharo.
Drat. What do I do to preserve indentation?
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
@Richard: code tags are what you want here – https://wordpress.com/support/wordpress-editor/blocks/syntax-highlighter-code-block/2/
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:
BTW, Schrage’s method for doing modular exponentiation without overflow is covered here: https://programmingpraxis.com/2014/01/14/
[…] Came across this on the “Programming Praxis” blog […]
A Mathematica version:
Part[Select[Prime /@ Range[10000], PalindromeQ], 100]Evaluates to:
94049
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).
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 retIt’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 retbut 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.
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");And it would seem, if anyone cares, that 9223372002002733229 is the largest palindromic prime less than 2^63.
Answer is 94949
A solution in python:
count = 0 number = 1 while count < 100: number += 1 if str(number) == str(number)[::-1]: is_prime = True for divisor in range(2, number): if number % divisor == 0: is_prime = False break if is_prime: count += 1 if count == 100: print(f'The 100th prime palindrome is: {number}')Examples: