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

Question

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
      mpz_add_ui(nr_superpalindromes_found[depth],nr_superpalindromes_found[depth],1);

      // 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++) {
    mpz_add(total_nr_superpalindromes,total_nr_superpalindromes,nr_superpalindromes_found[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).

Answer 1

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 reinventing-the-wheel.

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)
  

Comments:

• You have some obsolete comments that can be removed (old lines of code).

  // n_new[depth]:=n*p
  

• Use more comments to describe what you are doing. Since you are making more calls to an external library, this means that you will need more comments to explain why you are doing things the way you are, and how you accomplish that. You may also need to include comments on what the external functions purpose is.