Finding large primes

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}}
*)

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).

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