Primey Palindromes V2

This code takes just 0.25 secs to list all the palindromic primes less than 10000000000 that satisfy the given constraints.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243

    #include "iostream"
    #include "string"

    using namespace std;

    #define int64 long long int

    int64 power(int64 a, int64 n, int64 mod) {
        int64 power = a, result = 1;

        while (n) {
            if (n & 1) { 
                result = (result * power) % mod;
            }

            power = (power * power) % mod;
            n >>= 1;
        }

        return result;
    }

    bool witness(int64 a, int64 n) {
        int64 t, u, prev, curr;

        u = n / 2;
        t = 1;

        while (! ( u & 1)) {
            u /= 2;

            ++t;
        }

        prev = power(a, u, n);

        for (int64 i = 1; i <= t; ++i) {
            curr = (prev * prev) % n;

            if ((curr == 1) && (prev != 1) && (prev != (n - 1))) { 
                return true;
            }

            prev = curr;
        }

        if (curr != 1) { 
            return true;
        }

        return false;
    }

    bool is_prime(int64 number) {
        if (((! (number & 1)) && number != 2) || (number < 2) || ((number % 3 == 0) && (number != 3))) {
            return (false);
        }

        if (number < 1373653LL) {
            for (int64 k = 1; 36 * k * k - 12 * k < number; ++k) {
                if ((number % (6 * k + 1) == 0) || (number % (6 * k - 1) == 0)) {
                    return (false);
                }
            }

            return true;
        }

        if (number < 9080191LL) {
            if (witness(31, number)) {
                return false;
            }

            if (witness(73, number)) {
                return false;
            }

            return true;
        }

        if (witness(2, number)) {
            return false;
        }

        if (witness(7, number)) {
            return false;
        }

        if (witness(61, number)) {
            return false;
        }

        return true;
    }

    int64 next_palindrome(string a) {
        bool flag = false;

        int len = a.size();
        int carry = 0;

        for (int i = 0; i < len; i++) {
            if (a[i] == '9') {
            } else {
                flag = true;

                break;
            }
        }

        if (! flag) {
            for(int i = 1; i < len; i++) {
                a[i] = '0';
            }

            a[0] = '1';
            a += '1';

            int64 palin = 0;

            for (int i = 0; i < a.size(); i++) {
                palin = (10 * palin) + (a[i] - '0');
            }

            return palin;
        }

        for (int i = len - 1; i >= 0; i--) {
            if (a[i] == '9') {
                a[i] = '0';
            } else {
                a[i] += 1;

                break;
            }
        }

        int max = (len - 1) / 2;

        for (int i = 0; i <= max; i++) {
            if (a[i] >= a[len - 1 - i] + flag ? 1 : 0) {
                a[len - 1 - i] = a[i];
                flag = false;
            } else {
                int j;

                if (i == max) {
                    a[i] += 1;
                    a[len - 1 - i] = a[i];

                    if (a[i] == 58) {
                        j = i;
                        carry = 1;

                        while (carry) {
                            a[len - 1 - i] = '0';
                            a[i] = '0';

                            i--;

                            a[i] += 1;
                            a[len - i - 1] = a[i];

                            if(a[i] != 58) {
                                carry = 0;
                            }
                        }

                        i = j;
                    }
                }

                a[len - 1 - i] = a[i];
                flag = true;
            }
        }

        int64 palin = 0;

        for (int i = 0; i < a.size(); i++) {
            palin = (10 * palin) + (a[i] - '0');
        }

        return palin;
    }

    int64 digit_sum(int64 num) {
        string num_str = to_string(num);

        int64 sum = 0;

        for (int i = 0; i < num_str.size(); i++) {
            sum += (num_str[i] - '0');
        }

        return sum;
    }

    int main() {
        int t = 1;

        while (t--) {
            int digits[10] = {};

            int64 i = 9;

            for (; i < 10000000000;) {
                i = next_palindrome(to_string(i));

                while (! (i & 1)) {
                    i = next_palindrome(to_string(i));
                }

                if (is_prime(i)) {
                    if (is_prime(i + digit_sum(i))) {
                        digits[to_string(i).size()]++;
                    }
                }
            }

            for (int i = 2; i < 10; i++) {
                if (digits[i]) {
                    printf("Palindromic_Primes(%d digits) = %3d", i, digits[i]);

                    if (! is_prime(digits[i])) {
                        printf("[Not_Prime].\n", i);
                    } else {
                        printf("[Prime].\n", i);
                    }
                }
            }
        }

        return 0;
    }

    /*Output-
    **Palindromic_Primes_Count(2 digits) =   1[Not_Prime].
    **Palindromic_Primes_Count(3 digits) =   5[Prime].
    **Palindromic_Primes_Count(5 digits) =  13[Prime].
    **Palindromic_Primes_Count(7 digits) =  61[Prime].
    **Palindromic_Primes_Count(9 digits) = 364[Not_Prime].
    */

For every integer $N$ we can generate two unique palindromes. For example for $N=\overline { abc }$ (where $a,b,c$ are the base 10 digits of $N$ ) there exists two unique palindromes $\overline{abccba}$ and $\overline{abcba}$ . By iterating through all possible $N$ we can generate all palindromes.

Since we are dealing with prime palindromes here,we only need to generate palindromes with odd number of digits(The second type $\overline{abcba}$ ). This is because all prime palindromes except $11$ have an odd number of digits.

Using the method described above , we generate all such palindromes,test primality with Miller-Rabin,count them with a hashmap and print the first $n$ such that $g(n)$ is not prime.

Code in Java

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

import java.math.BigInteger;
import java.util.*;
import java.util.HashMap;

public class main {
    public static BigInteger Make_Palindrome(int n){
        StringBuffer Sn = new StringBuffer(""+n);
        String A = new String(Sn) + new String(Sn.reverse());
        int p = Sn.length();
        String P = A.substring(0 , p) + A.substring(p+1,p*2);
        return new BigInteger(P);

    }
    public static BigInteger f(BigInteger x){
        BigInteger Sum = BigInteger.ZERO;
        String str = "" + x;
        for (int i = 0 ; i < str.length() ; i ++){
            Sum  = Sum.add(new BigInteger(""+str.charAt(i)));
        }
        return Sum;
    }

    public static void main(String[] args){
        HashMap<Integer, Integer> cache = new HashMap<Integer, Integer>();
        for (int n = 0 ; n < 100000 ; n++){          
            BigInteger P = Make_Palindrome(n);
            int Size = P.toString().length();
            boolean Chek = P.isProbablePrime(7) && P.add(f(P)).isProbablePrime(7);
            if (Chek){
                if (cache.containsKey(Size)){
                    cache.put(Size,cache.get(Size)+1);
                }
                else{
                    cache.put(Size, 1);
                }
            }
        }
        Iterator<Integer> keySetIterator = cache.keySet().iterator();
        while (keySetIterator.hasNext()){
            Integer Key = keySetIterator.next();
            int Gn = cache.get(Key);
            boolean IsP =  new BigInteger(""+Gn).isProbablePrime(7);
            if (Gn!=0 && !IsP){
                System.out.println("Answer is "+Key);
                break;
            }

        }
    }


}

First build a sieve, then iterate over all palindromes

import math

primes = [2]

def is_prime(number):
  # primes are always 1, 5 mod 6
  if number % 6 not in (1, 5):
    return False
  if number in primes:
    return True
  for prime in primes:
    if prime > math.sqrt(number):
      return True
    if number % prime == 0:
      return False
  return True

def gen_primes(top):
  for candidate in range(2, top):
    if is_prime(candidate):
      primes.append(candidate)

gen_primes(100000)
assert is_prime(383)
assert is_prime(70607)

def sumdigits(number):
  return sum(map(int, str(number)))

# each palindrome is made of a matching left and right

def palindrome_lefts(length):
  start = 10 ** (length-1)
  for x in xrange(start, 10 ** length):
    yield str(x)

def even_palindromes(length):
  assert length % 2 == 0
  for left in palindrome_lefts(length/2):
    yield left + ''.join(reversed(left))

def palindromes(length):
  if length % 2 == 0:
    for item in even_palindromes(length):
      yield int(item)
    return

  for palindrome in even_palindromes(length-1):
    for i in xrange(10):
      mid = length/2
      yield int(palindrome[:mid] + str(i) + palindrome[mid:])

def matches_in_n(length):
  for item in palindromes(length):
    if is_prime(item) and is_prime(item + sumdigits(item)):
      yield item

for index, match in enumerate(matches_in_n(9)):
  print index+1, match
print 'finished'

I used two codes (both written in Python), first is the code for Miller-Rabin Primality Test (used to check if number is prime):

#Filename: MillerRabin.py

import random

def decompose(n):
   exponentOfTwo = 0

   while n % 2 == 0:
      n = n/2
      exponentOfTwo += 1

   return exponentOfTwo, n

def isWitness(possibleWitness, p, exponent, remainder):
   possibleWitness = pow(possibleWitness, remainder, p)

   if possibleWitness == 1 or possibleWitness == p - 1:
      return False

   for _ in range(exponent):
      possibleWitness = pow(possibleWitness, 2, p)

      if possibleWitness == p - 1:
         return False

   return True

def probablyPrime(p, accuracy=100):
   if p == 2 or p == 3: return True
   if p < 2: return False

   numTries = 0
   exponent, remainder = decompose(p - 1)

   for _ in range(accuracy):
      possibleWitness = random.randint(2, p - 2)
      if isWitness(possibleWitness, p, exponent, remainder):
         return False

   return True

and the second is for the solution itself:

import MillerRabin as PrimalityTest #used for checking if number is prime

def sumOfDigits(num):
    num = (str)(num)
    sum = 0

    for digit in num:
        numberValue = (int)(digit)
        sum += numberValue

    return sum

def numOfNDigitPrimeyPalindromes(numDigits):
    nDigitPalindromes = getNDigitPalindromes(numDigits)

    numPrimeyPalindromes = 0

    for num in nDigitPalindromes:
        if isPrimey(num):
            numPrimeyPalindromes += 1

    return numPrimeyPalindromes

def getNDigitPalindromes(numDigits):
    #Palindrome = leftNo + (midNo) + rightNo
    #rightNo = mirrorString(leftNo)

    nDigitPalindromes = []
    digits = [0,1,2,3,4,5,6,7,8,9]

    numDigitsLeftNo = numDigits/2

    min = 10**(numDigitsLeftNo-1)
    max = min*10

    for leftNo in range (min,max):
        if numDigits % 2 == 0:
            palindrome = (str)(leftNo) + mirrorString(leftNo)
            nDigitPalindromes.append((int)(palindrome))
        else:
            for digit in digits:
                palindrome = (str)(leftNo) + (str)(digit) + mirrorString(leftNo)
                nDigitPalindromes.append((int)(palindrome))

    return nDigitPalindromes

def mirrorString(string):
    string = (str)(string)
    temp = ''

    for character in string:
        temp = character + temp

    return temp

def isPrimey(num):
    a = num
    f = sumOfDigits(a)

    return PrimalityTest.probablyPrime(a) and PrimalityTest.probablyPrime(a+f)

n = 3

while True:
    g = numOfNDigitPrimeyPalindromes(n)
    isPrime = PrimalityTest.probablyPrime(g)

    print ''
    print 'n:\t',n
    print 'g(n):\t',g
    print 'Prime:\t',isPrime
    print ''

    if g != 0 and not isPrime:
        break

    n += 1

Total[Total[ Boole[PrimeQ[#] && PrimeQ[# + Total[IntegerDigits[#]]]] & /@ FromDigits /@ Table[Flatten[{IntegerDigits[#], i, Reverse[IntegerDigits[#]]}], {i, 0, 9}] & /@ Range[1001, 9999, 2]]]

Gives an output of 177 9 digit palindromes in 0.36seconds.