11
$\begingroup$

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.

Improve this question
$\endgroup$
7
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Michael E2
    Jan 9, 2016 at 14:12
  • $\begingroup$ Is there any theory to help you find them? A direct search is apt to be slow. $\endgroup$
    – Michael E2
    Jan 9, 2016 at 14:27
  • $\begingroup$ what do you consider large? $\endgroup$
    – george2079
    Jan 9, 2016 at 14:29
  • $\begingroup$ somewhere in the range of 10^12+ $\endgroup$
    – aL_eX
    Jan 9, 2016 at 14:50
  • $\begingroup$ You need to use @george2079 to make sure the commenter is notified of your response. Authors of posts are always notified. $\endgroup$
    – Michael E2
    Jan 9, 2016 at 15:34

5 Answers 5

Reset to default
7
$\begingroup$

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}}
*)
Improve this answer
$\endgroup$
2
  • $\begingroup$ Amazing! Thank you very much. $\endgroup$
    – aL_eX
    Jan 9, 2016 at 16:20
  • $\begingroup$ @AlexanderMilanovic You're welcome. $\endgroup$
    – Michael E2
    Jan 9, 2016 at 16:48
4
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • $\begingroup$ Yeah, I get my 85-sec list, plus more, in 1.2 sec, with Do[..., {j, 4^12 + 1, 4^12 + 4^9, 2}] instead of For and skipping even numbers. $\endgroup$
    – Michael E2
    Jan 10, 2016 at 3:24
  • $\begingroup$ @MichaelE2 Originally I coded mod 3, where one cannot skip evens. When I realized mod 4 was required I failed to notice I could now skip the evens. Oh well. $\endgroup$
    – Daniel Lichtblau
    Jan 10, 2016 at 16:19
  • $\begingroup$ Well, I dismissed your approach at first as infeasible because I didn't think of the IntegerDigits[j, 4]] trick. I hope people eventually recognize how superior this answer this. $\endgroup$
    – Michael E2
    Jan 10, 2016 at 16:55
3
$\begingroup$

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

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for the answer. Is there any way to parallel compute the problem as only one Mathematica Kernel is used during the computation? $\endgroup$
    – aL_eX
    Jan 9, 2016 at 14:36
  • 2
    $\begingroup$ I found that using ParallelTable instead of Table in the first line seems to speed up the process $\endgroup$
    – aL_eX
    Jan 9, 2016 at 14:51
2
$\begingroup$

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.
Improve this answer
$\endgroup$
7
  • $\begingroup$ Actually the method I showed will handle ranges over large numbers just fine. The percentage of hits will drop of course, in accordance with the Prime Number Theorem. (If it is not in accordance with PNT, that would itself be extremely interesting, to say the least). $\endgroup$
    – Daniel Lichtblau
    Jan 10, 2016 at 16:35
  • $\begingroup$ @DanielLichtblau : A difficulty with PrimeQ[] is that it caches results, but not its sieve between calls. (It expects random access over the integers.) In contrast, NextPrime[] caches its "small primes" sieve so that subsequent calls can skip sieving. (It expects sequential access over the integers.) This may not show up in the above examples (10-200 digits), but I would expect to see a difference somewhere above 500 digits on recent hardware. $\endgroup$
    – Eric Towers
    Jan 10, 2016 at 17:15
  • $\begingroup$ Prime uses a sieve. PrimeQ, which is the engine behind NextPrime, does not sieve anything. I'm not sure what you mean by a "small primes" sieve, but regardless I don't see an advantage to using NexPrime for this problem since it ends up doing just as many PrimeQ tests, if not more, than the method I had posted (and those are going to be a bottleneck). $\endgroup$
    – Daniel Lichtblau
    Jan 10, 2016 at 18:01
  • $\begingroup$ @DanielLichtblau : Rather than argue... The documentation states: "NextPrime[n] finds the smallest prime $p$ such that $p>n$. For n less than 20 digits, the algorithm does a direct search using PrimeQ on the odd numbers greater than $n$. For $n$ with more than 20 digits, the algorithm builds a small sieve and first checks to see whether the candidate prime is divisible by a small prime before using PrimeQ. This seems to be slightly faster than a direct search." So NextPrime[] for big primes. $\endgroup$
    – Eric Towers
    Jan 10, 2016 at 23:46
  • $\begingroup$ (1) I should check that. I suspect it just means a GCD check is done with a product of several small primes before going to PrimeQ. It's been a few years since I've looked closely, but I don't recall anything really sievelike in there. $\endgroup$
    – Daniel Lichtblau
    Jan 11, 2016 at 4:42
1
$\begingroup$

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]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.