Depth-first search method for searching for all “superpalindromes” whose prime factors are all <=N

A question was asked over at math.SE (here) about whether or not there are infinitely many superpalindromes. I'm going to rephrase the definition to a more suitable one for coding purposes:

Definition: A superpalindrome is a product of primes p(1) * p(2) * ... * p(k), where p(i)<=p(i+1) for 1<=i<=k-1, such that p(1) * p(2) * ... * p(r) is a palindrome for all 1 <= r <= k.

A natural question to ask, is whether or not there's an infinite number of superpalindromes with all of its prime factors <=N. Thus, we can use a depth-first search.

Here is my implementation in C using GMP, but I'm wondering if there is some ways to give non-trivial run-time improvements.

#include <stdio.h>
#include <gmp.h>

// search for superpalindromes containing prime factors in {2,3,...,q}
// where q is the MAX_NR_PRIMES-th prime
#define MAX_NR_PRIMES 20000

#define MAX_SEARCH_DEPTH 100
#define MAX_NR_DIGITS 10000

// start backtracking at the (STARTING_PRIME_ID+1)-th prime
#define STARTING_PRIME_ID 0

mpz_t list_of_small_primes[MAX_NR_PRIMES];
mpz_t n_new[MAX_SEARCH_DEPTH];
int max_depth_reached=0;

mpz_t nr_superpalindromes_found[MAX_SEARCH_DEPTH];

// checks if n is a palindrome
// idea "borrowed" from http://gmplib.org/list-archives/gmp-discuss/2012-February/004876.html
int is_palindrome(mpz_t n) {
char m[MAX_NR_DIGITS];
int len=gmp_sprintf(m,"%Zd",n);
for(int i=0;i<len;i++) {
if(m[i]!=m[len-i-1]) return 0;
}
return 1;
}

// depth-first search for superpalindromes
// searches for a prime p:=list_of_small_primes[i_new], with i_new>=i
// such that n*p is a palindrome; continues searching if one is found
int extend_superpalindrome_backtracking_algorithm(mpz_t n,int i,int depth) {
for(int i_new=i;i_new<MAX_NR_PRIMES;i_new++) {
// n_new[depth]:=n*p
mpz_mul(n_new[depth],n,list_of_small_primes[i_new]);

if(is_palindrome(n_new[depth])) {
// increment count of number of superpalindromes with depth+1 prime factors

// print out the first superpalindrome found with >depth+1 prime factors
if(depth>max_depth_reached) {
max_depth_reached=depth;
gmp_printf("superpalindrome found: %Zd at depth %d: ",n_new[depth],depth+1);
mpz_t primes[MAX_SEARCH_DEPTH];
for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_init(primes[i]);
mpz_set(primes[0],n_new[0]);
for(int i=1;i<=depth;i++) mpz_divexact(primes[i],n_new[i],n_new[i-1]);
for(int i=0;i<=depth;i++) gmp_printf("%Zd ",primes[i]);
printf("\n");
for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_clear(primes[i]);
}

// continue the depth-first search if superpalindrome found
extend_superpalindrome_backtracking_algorithm(n_new[depth],i_new,depth+1);
}
}
}

int main() {
mpz_t n,p;
mpz_init_set_ui(n,1);
mpz_init_set(p,list_of_small_primes[STARTING_PRIME_ID]);
for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_init(n_new[i]);
for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_init(nr_superpalindromes_found[i]);
for(int i=0;i<MAX_NR_PRIMES;i++) mpz_init(list_of_small_primes[i]);

// pre-compute small primes
mpz_set_ui(list_of_small_primes[0],2);
for(int i=1;i<MAX_NR_PRIMES;i++) mpz_nextprime(list_of_small_primes[i],list_of_small_primes[i-1]);

extend_superpalindrome_backtracking_algorithm(n,STARTING_PRIME_ID,0);

// output results
gmp_printf("\nfound all superpalindromes with prime factors in {2,3,...,%Zd}\n",list_of_small_primes[MAX_NR_PRIMES-1]);
mpz_t total_nr_superpalindromes;
mpz_init_set_ui(total_nr_superpalindromes,0);
for(int i=0;i<=max_depth_reached;i++) {
gmp_printf("nr superpalindromes with %d prime factors: %Zd\n",i+1,nr_superpalindromes_found[i]);
}
gmp_printf("total nr superpalindromes: %Zd\n",total_nr_superpalindromes);
mpz_clear(total_nr_superpalindromes);

for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_clear(n_new[i]);
for(int i=0;i<MAX_NR_PRIMES;i++) mpz_clear(list_of_small_primes[i]);
mpz_clear(n);
mpz_clear(p);
return 0;
}


Here's the output:

superpalindrome found: 4 at depth 2: 2 2
superpalindrome found: 8 at depth 3: 2 2 2
superpalindrome found: 88 at depth 4: 2 2 2 11
superpalindrome found: 2552 at depth 5: 2 2 2 11 29
superpalindrome found: 257752 at depth 6: 2 2 2 11 29 101
superpalindrome found: 67788776 at depth 7: 2 2 2 11 29 101 263
superpalindrome found: 616267762616 at depth 8: 2 2 2 11 29 101 263 9091
superpalindrome found: 6101667117661016 at depth 9: 2 2 2 11 29 101 263 9091 9901
superpalindrome found: 20302629368699686392620302 at depth 10: 2 11 11 11 101 9091 9091 9091 9901 10151
superpalindrome found: 1001004004006006004004001001 at depth 11: 7 11 13 101 101 101 101 9901 9901 9901 9901

found all superpalindromes with prime factors in {2,3,...,224737}
nr superpalindromes with 1 prime factors: 113
nr superpalindromes with 2 prime factors: 428
nr superpalindromes with 3 prime factors: 1022
nr superpalindromes with 4 prime factors: 1539
nr superpalindromes with 5 prime factors: 1603
nr superpalindromes with 6 prime factors: 1137
nr superpalindromes with 7 prime factors: 565
nr superpalindromes with 8 prime factors: 217
nr superpalindromes with 9 prime factors: 50
nr superpalindromes with 10 prime factors: 13
nr superpalindromes with 11 prime factors: 1
total nr superpalindromes: 6688


Currently, I'm mostly concerned about the efficiency of is_palindrome(mpz_t n) since it feels like a "two-pass" method: (a) copy the number to a string, (b) check if the string is a palindrome. But please highlight any other areas I might not be concerned about but should be.

Most of the numbers encountered by is_palindrome(mpz_t n) will not be palindromes, but my attempt at checking only the first and last digits was thwarted by mpz_sizeinbase not giving the exact number of digits in base 10 (and it added too much overhead to add a "check if the number of digits is correct" function).

Things you did well

• You defined i within your for loops, and abided by the C99 standards.

• You used an external library instead of writing your (more likely inefficient) own methods, and thus avoided unnecessarily .

Things you could improve on:

Efficiency:

• Copying a whole number into a external array with gmp_sprintf() in your is_palindrome() method is very time consuming, as you suspected.

int is_palindrome(mpz_t n) {
char m[MAX_NR_DIGITS];
int len=gmp_sprintf(m,"%Zd",n);
for(int i=0;i<len;i++) {
if(m[i]!=m[len-i-1]) return 0;
}
return 1;
}


Is it reasonable to assume that most numbers will be non-palindromes? If not, I would suggest something like this:

1. Determine the approximate number of digits in your number a using mpz_sizeinbase() and call this number $n$.

2. Compute $\dfrac{a}{d^{(n-30)}}$, then compute $a \bmod d^{30}$. Print both in base $d$ and compare. If they are unequal, reject as non-palindrome.

Computing $\dfrac{a}{d^{(n-30)}}$ will be fast using mpz_tdiv_q(). Unfortunately, computing $d^{30}$ is not fast. Perhaps a rough approximation of that would also increase efficiency.

• Using a sieve for an early exit would increase efficiency.

• Combine some of your for loops.

for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_init(n_new[i]);
for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_init(nr_superpalindromes_found[i]);
...
for(int i=1;i<=depth;i++) mpz_divexact(primes[i],n_new[i],n_new[i-1]);
for(int i=0;i<=depth;i++) gmp_printf("%Zd ",primes[i]);
...
for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_init(primes[i]);
for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_clear(primes[i]);


These loops are very similar, and many iterations (even for simple loops) can take more time than you think. Combining them will increase efficiency.

• Try to avoid putting statements that print data to the console inside of a loop.

for(int i=0;i<=depth;i++) gmp_printf("%Zd ",primes[i]);


gmp_printf() (and all printf()-esque statements) can be very taxing on a system. Extracting them to the outside of the loop will increase efficiency.

Style:

• Your code is very compact.

for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_init(primes[i]);
mpz_set(primes[0],n_new[0]);
for(int i=1;i<=depth;i++) mpz_divexact(primes[i],n_new[i],n_new[i-1]);
for(int i=0;i<=depth;i++) gmp_printf("%Zd ",primes[i]);
printf("\n");
for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_clear(primes[i]);


I'm all for compact code, but you still want your code to be readable. Let it breathe a bit. Use spaces after commas and semi-colons, use an empty line to separate and organize code.

• You sometimes write a loop with braces, and sometimes without braces, even though they both have one statement in them.

for(int i=0;i<len;i++) {
if(m[i]!=m[len-i-1]) return 0;
}
...
for(int i=0;i<MAX_SEARCH_DEPTH;i++) mpz_init(primes[i]);


Choose one way and stick with it. Consistency is very important.

Syntax

• If you don't accept any parameters into a function, they should be declared void.

int main(void)


// n_new[depth]:=n*p