Finding large primes

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 1000 and 20001 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.

Answer 1

For machine-sized primes (< 2^63 ~= 10^21), here's a simple compiled program to detect a good integer (digits < 4, mostly zeros). You can leave off CompilationTarget -> "C" if you don't have a C compiler; it doesn't help much on my machine anyway.

good = Compile[{{p, _Integer}},
  Module[{n = p, q, r, z = 0, nz = 0},
   While[n > 0,
    q = Quotient[n, 10];
    r = n - 10*q;
    If[r > 3, Return[False]];
    If[r == 0, z++, nz++];
    n = q];
   Return[z > nz]],
  RuntimeAttributes -> {Listable}, Parallelization -> True,
  CompilationTarget -> "C" (* optional; helps only a little *)
  ]

Here's a basic use:

Pick[#, good@#] &@ ParallelMap[Prime, Range[1, 1000000]] // AbsoluteTiming
(*
  {1.01688,
   {100003, 200003, 1000003, 1000033, 1000303, 1001003, 
    1003001, 1003003, 1010003, 1020001, 2000003, 2000303, 2002001, 
    2020001, 2020003, 2100001, 2300003, 3000103, 3000301, 3001001, 
    3001003, 3002003, 3010001, 3200003, 3300001, 10000103, 10000303, 
    10003001, 10003003, 10030003}}
*)

The bottleneck is going to be Prime as far as speed goes. Prime tries to be efficient. It caches results and reuses them as starting points.

Quit[]
Prime[PrimePi[10^12]] // AbsoluteTiming
Prime[PrimePi[10^12]] // AbsoluteTiming
Prime[Range[#, # + 10] &@PrimePi[10^12]] // AbsoluteTiming
(*
  {1.62425, 999999999989}
  {0.005789, 999999999989}
  {0.003345,
   {999999999989, 1000000000039, 1000000000061, 1000000000063,
    1000000000091, 1000000000121, 1000000000163, 1000000000169, 
    1000000000177, 1000000000189, 1000000000193}}
*)

Notice the work that went into computing the first result. When you parallelize Prime, it has to do all that work on each kernel. But since it does it simultaneously, the user doesn't lose. Note that computing each successive prime takes an appreciable amount of time, so this problem is a good candidate for parallelization.

Here is the timing in the neighborhood of the OP's interest.

Quit[]

(* Be sure to reevaluate the definition of good! *)

Pick[#, good@#] &@
  ParallelMap[Prime, Range[PrimePi[10^12], PrimePi[10^12] + 10^7]] // AbsoluteTiming
(*
  {85.3418,
   {1000000000121, 1000000000211, 1000000000303, 1000000000331,
    1000000000333, 1000000001123, 1000000001201, 1000000001213, 
    ... ~1000 primes omitted ..
    1000232020001, 1000233000001, 1000233000031, 1000233000101, 
    1000233000203, 1000233003001, 1000233010003}}
*)

Answer 2

I suspect it is faster to form digit-wise viable candidates, and cull the ones that are PrimeQ, than to iterate over explicit primes. The reason is that the percentage of the latter is going to be much smaller than that of the former, hence we can save considerably on the number of iterations and tests. The code below implements this approach.

Timing[Reap[For[j = 1, j < 10^5, j++,
    digits = IntegerDigits[FromDigits[IntegerDigits[j, 4]]];
    If[Last[digits] == 0 || Last[digits] == 2, Continue[]];
    If[Count[digits, 0] < Length[digits]/2, Continue[]];
    val = FromDigits[digits];
    If[PrimeQ[val], Sow[val]];
    ]][[2, 1]]]


(* Out[118]= {0.494928, {2003, 3001, 100003, 100103, 102001, 103001, 
  130003, 200003, 200023, 200033, 200201, 202001, 230003, 300023, 
  300301, 1000003, 1000033, 1000303, 1001003, 1003001, 1003003, 
  1010003, 1020001, 2000003, 2000303, 2002001, 2020001, 2020003, 
  2100001, 2300003, 3000103, 3000301, 3001001, 3001003, 3002003, 
  3010001, 3200003, 3300001, 10000103, 10000121, 10000223, 10000303, 
  10001203, 10002203, 10003001, 10003003, 10003031, 10010023, 
  10010101, 10012001, 10020013, 10020103, 10021001, 10021003, 
  10023001, 10030003, 10030103, 10033003, 10100011, 10100201, 
  10102003, 10200011, 10200101, 10200301, 10201001, 10203001, 
  10220003, 10300013, 10300201, 11002001, 11020001, 11020003, 
  11200001, 12001001, 12100003, 13000201, 13000301, 13000303, 
  13200001, 20000003, 20000023, 20000033, 20000213, 20000221, 
  20000303, 20000311, 20001001, 20001031, 20001203, 20001301, 
  20002201, 20003003, 20003023, 20003201, 20011001, 20020001, 
  20020103, 20020303, 20023001, 20030011, 20030203, 20100023, 
  20100203, 20101001, 20102003, 20130001, 20200013, 20203003, 
  20300003, 20300011, 20300101, 20300201, 20310001, 21000103, 
  21003001, 21030001, 21100003, 22000001, 22000003, 22000021, 
  22000031, 22000103, 22000201, 22010003, 22200001, 22300001, 
  23000011, 23000101, 23100001, 30000001, 30000023, 30000133, 
  30000301, 30000323, 30000331, 30001003, 30001031, 30001201, 
  30001303, 30002023, 30002033, 30002303, 30003101, 30003301, 
  30010021, 30010033, 30010103, 30010301, 30010303, 30013003, 
  30030023, 30030031, 30030101, 30100003, 30100201, 30100303, 
  30101003, 30103001, 30130003, 30200003, 30200011, 30200021, 
  30200033, 30200201, 30210001, 30300301, 30303001, 30310003, 
  31000003, 31000301, 31000303, 31002001, 31020001, 31030001, 
  32000011, 32000303, 32001001, 32010001, 32030003, 33000001, 
  33000031, 33000103, 100000123, 100000213, 100000223, 100000231, 
  100001203, 100001303, 100002011, 100002013, 100002031, 100002103, 
  100003021, 100003301, 100010021, 100010023, 100010033, 100010203, 
  100020023, 100020103, 100021003, 100022003, 100030001, 100031003, 
  100120003, 100130003, 100200011, 100200013, 100200031, 100200301, 
  100203001, 100210001, 100230001, 100300001, 100303003, 100310003, 
  100330001, 100330003, 101000023, 101000203, 101001001, 102000023, 
  102002003, 102100001, 102100003, 102300001, 103000201, 103000301, 
  103010003, 103300003, 110000201, 110002003, 110020003, 120000031, 
  120000103, 120000203, 120000301, 120001001, 120003001, 120100003}} *)

An iteration to 10^7 takes around 52 seconds and delivers 10151 such values (also 10501 is a prime, though not quite a contender; too bad it wasn't 10301).

Answer 3

For example

pri = IntegerDigits /@ Table[Prime[n], {n, 20000, 90000}];

res = FromDigits /@ Select[pri, Count[#, 0] > Floor[Length[#]/2] && Max[#] < 4 &]

{1000003, 1000033, 1000303, 1001003, 1003001, 1003003, 1010003, 1020001}

AllTrue[PrimeQ /@ res, TrueQ]

True

Answer 4

The other answers here are satisfactory for your claimed range of interest, around $10^{12}$. However, if you want actually large primes, I'd suggest a different method, using NextPrime[].

random0123[length_, zeroFraction_] := Module[{
    numZeroes,
    startDigits, startNumber,
    rightDigits
  },
  numZeroes = Floor[zeroFraction*length];
  startDigits = Join[
    {RandomChoice[{1,2,3}]},
    Permute[Join[
        RandomChoice[{1,2,3},length-numZeroes-1],
        Table[0,{numZeroes}]
      ],
      RandomPermutation[SymmetricGroup[length-1]]
    ]
  ];
  startNumber = NextPrime[FromDigits[startDigits]];

  rightDigits[number_] := SubsetQ[{0,1,2,3},DeleteDuplicates[IntegerDigits[number,10]]];

  While[Not[rightDigits[startNumber]], 
    startNumber = NextPrime[startNumber];
  ];
  startNumber
]

This can generate somewhat larger results rather quickly, but the time can be highly variable. (This variability is not clearly indicated below.) If you're going to generate a lot of large primes, I'd come back to this and occasionally restart with a new set of startDigits. (There don't seem to be any over here. Let's try looking over there.)

RepeatedTiming[random0123[50,0.6]]
RepeatedTiming[random0123[100,0.6]]
RepeatedTiming[random0123[200,0.8]]
RepeatedTiming[random0123[200,0.2]]

(*  First number is average seconds, second number is the first prime it found.
  {0.02 , 10300220300000013020000000300110321000031302033231}
  {0.02 , 2000210102000000023103001031030300001000000203121000203000000120010020100320002300002001103230101021}
  {2.   , 30020300013030000021000000000003010002000030003003000000000000010000000000000033000033000100003000000000000000000030000100300000100000100030300000000000000000300200300000000000000101300120001000032103}
  {0.12 , 33220313233111203030200221002001211011010021331233310031111200312331021002103133323012223221232132101123113121023010111232133002002132222110213111312301013023313121301131322333131312121132021231203311}
*)

This function determines the number of desired 0s and selects a uniformly randomly chosen, with replacement, collection of 1s, 2s, and 3s to fill out to one fewer digit than requested. (They're randomly ordered. They don't need to be. If we were going to generate a few million digits, we shouldn't randomize here; we should let the upcoming permutation do that work.) Then a sequence of digits is constructed that is first, a random choice of nonzero digit, followed by a random permutation of the desired number of zeroes and the collection of nonzero digits just constructed. Next we define a test function that verifies that a given number contains only the desired digits. Finally we start searching from the generated number sequentially through the subsequent primes until we find one containing only the desired digits.

There are a number of obvious elaborations of this function:

• Specify number of zeroes instead of fraction of zeroes.
• Specify range of allowed digits or range of allowed number or fraction of zeroes.
• Specify how many results we want.
• Specify the leading digit to come from a different set of digits.
• Specify weights for the various digits. (Perhaps 1s are as likely as all other digits put together.)
• Supply specific counts ranges for each digit.
• Allow bases other than 10.
• Set a timeout for restarting the search from a new startNumber.
• Set an upper limit to the number of times to march on from a given startNumber before generating another one.
• Setup parallelization to simultaneously run multiple searches from different startNumbers.
• Keep statistics on number of times through the While[] loop and/or number of startNumbers used.
• Perhaps don't stop at the first valid prime in the While[] loop, but collect a sequence of successive valid primes, to get statistics on how frequent primes satisfying the given criteria are.

Answer 5

construct random numbers meeting the critera and use primeQ. This produces a handful of 2000 digit primes in about a minute.

len = 2000;
Union@Select[Table[
   nz = RandomInteger[{Ceiling[(len + 1)/2], len - 2}];
     FromDigits@
    Join[
     RandomChoice[{1, 2, 3}, 1],
     RandomSample[
      ConstantArray[0, nz]~Join~RandomChoice[{1, 2, 3}, len - nz - 2],
       len - 2],
     RandomChoice[{1, 3}, 1]], {2000}], PrimeQ]