Finding large primes in Mathematica

Question

I'm quite new to Mathematica and I am trying to find large prime numbers that can be written using only the digits 0, 1, 2 and 3 and more than half of these digits have to be 0. For example 100003 would qualify.

I was thinking of using the Prime[n] function with a bunch if "If" statements, but I was wondering if there is a more efficient way of going about this. Thanks in advance.

Answer 1

For a good answer ask this question at http://mathematica.stackexchange.com

For an answer of sorts, read on ...

The function

myTestQ[num_Integer] := 
   And[DigitCount[num][[10]] > Plus @@ DigitCount[num][[1 ;; 3]], PrimeQ[num]]

returns True when it's passed an integer which (i) has more 0 digits than 1s,2s, and 3s and (ii) is prime. And evaluates its arguments in order and short-circuits so should avoid testing the primality of a number whose digit counts are not right.

In practice your computation time will be dominated by primality testing; generating and throwing away a lot of integers with the wrong digit counts is inexpensive enough that you can ignore the inefficiencies in that part of the program.

Now to test a range of integers composed of only 0,1,2,3 digits. There's a whole family of integers using only those digits, namely all integers when written in base-4. So lets generate those:

IntegerString[Range[lo, hi, 2], 4]

This will generate strings representing all the base-4 integers from the lo-th to the hi-th in steps of 2 (you're not going to be interested in any even numbers). For example

IntegerString[Range[1, 13, 2], 4]

produces

{"1", "3", "11", "13", "21", "23", "31"}

that is, the 1st, 3rd, 5th, ..., 13th base-4 numbers in string form. Of course, you'll need them back as integers to test them, and we'll use ToExpression for that.

Mash it all together and you get something like

Select[ToExpression[IntegerString[Range[lo, hi, 2], 4]], myTestQ]

replacing lo and hi with whatever you like, but making sure lo is an odd number.

Answer 2

You could go fishing. First generate large integers of the required pattern and then test them for primality. I don't have Mathematica, but here is a quick and dirty implementation in Python:

import random

def rand_candidate(n):
    #assumes n > 4
    while True:
        s = [random.choice(['1','2','3'])]
        for i in range(n-2):
            s.append(random.choice(['0','0','0','1','2','3']))
        s.append(random.choice(['1','3']))
        if s.count('0') > n/2: return int(''.join(s))

The previous function starts by selecting a non-zero leading digit, then picks the intermediate digits in a way which is biased to 0, then picks the final digit to make it odd, checking that the overall result has the requisite number of zeros.

For a primality test, I got lazy and just used a Fermat pseudo-prime test. The following test was used (and is maybe still used) for the PGP implementation of RSA. It has a very low probability of certifying a composite number as prime. In Mathematica you would of course have recourse to more sophisticated tests:

def prob_prime(n):
    tests = [2,3,5,7]
    return all(n % p != 0 and pow(p,n-1,n) == 1 for p in tests)

Putting them together:

def go_fishing(n):
    while True:
        p = rand_candidate(n)
        if prob_prime(p): return p

Typical output:

>>> go_fishing(20)
23002001032200000101
>>> go_fishing(100)
1100200103020002000000000020232222230321210000103030031022000000200001001300020010012122020101020113

The computer algebra system Derive (which I do have on my machine) certifies both of the above as prime.