Can you answer my son's fourth-grade homework question: Which numbers are prime, have digits adding to ten and have a three in the tens place? My son Horatio (nine years old, fourth grade) came home with some fun math homework exercises today. One of his problems was the following little question: 


  I am thinking of a number...
  
  
  It is prime.
  The digits add up to $10.$ 
  It has a $3$ in the tens place.
  
  
  What is my number? 


Let us assume that the problem refers to digits in decimal notation. Horatio came up with $37,$ of course, and asked me whether there might be larger solutions with more digits. We observed together that $433$ is another solution, and also $631$ and $1531.$ But also notice that $10333$ solves the problem, based on the list of the first $10000$ primes, and also $100333$, and presumably many others. 

My question is: How many solutions does the problem have? In particular, are there infinitely many solutions? 

How could one prove or refute such a thing? I could imagine that there are very large prime numbers of the decimal form $10000000000000\cdots00000333$, but don't know how to prove or refute this.

Can you provide a satisfactory answer this fourth-grade homework question?
As requested I'm posting this an answer. I wrote a short sage script to check the primality of numbers of the form $10^n+333$ where $n$ is in the range $[4,2000]$. I found that the following values of $n$ give rise to prime numbers:

$$4,5,6,12,53,222,231,416.$$

Edit 3: I stopped my laptop's search between 2000 and 3000, since it hadn't found anything in 20 minutes. I wrote a quick program to check numbers of the form $10^n+3*10^i+33$. Here are a couple 


100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000030000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000033
100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000300033
100000000000000000000000000000000000000000000000000000300000000000000000000000000000000000000033
100000000000000000000000000000000000000000000000030000000000000000000000000000000000000000000033
100000000000000000000000000000000000000000000030000000000000000000000000000000000000000000000033
10000000000000000000000000000000003000000033
10000000000000000000000000000030000000000033
10000000000000000000000030000000000000000033
10000000003000000000000000000000000000000033


There seemed to be plenty of numbers of this form and presumably I could find more if I checked some of the other possible forms as outlined by dr jimbob. 

Note: I revised the post a bit after jimbob pointed out I was actually looking for primes that didn't quite fit the requirements. 

Edit 4: As requested here are the sage scripts I used. 
To check if $10^n+333$ was prime:

for n in range(0,500):
  k=10^n+333
  if(is_prime(k)):
    print n


And to check for numbers of the form $10^n+3*10^i+33$:

for n in range(0,500):
  k=10^n+33
  for i in range(2,n):
    l=k+3*10^i
    if(is_prime(l)):
      print l
This question was found while searching for mathematical recreational problems to use for didactic reasons in programming. More exactly in sage (computer algebra system = CAS used for structural mathematical questions, based on python) for different aspects of using or not using (python native) &quot;iter tools&quot;. The idea was to illustrate instances of the classes like Partitions and Combinations. Since sage also comes with a quick detection of prime number if it has a &quot;decent size&quot;, i tried to write a short code snippet for this problem. I will insert first some dialog with the python interpreter, showing why (and how) the constructors Partitions and Combinations can be useful, then the compact code to search for prime numbers satisfying the conditions in the OP.

First of all, finding partitions of ten with parts $\ge 1$ and at least one piece equal to $3$ is simple:
partitions = [part for part in Partitions(10) if 3 in part]

Let us show them in the interpreter:
sage: for part in partitions:
....:     print(part)
....:
[7, 3]
[6, 3, 1]
[5, 3, 2]
[5, 3, 1, 1]
[4, 3, 3]
[4, 3, 2, 1]
[4, 3, 1, 1, 1]
[3, 3, 3, 1]
[3, 3, 2, 2]
[3, 3, 2, 1, 1]
[3, 3, 1, 1, 1, 1]
[3, 2, 2, 2, 1]
[3, 2, 2, 1, 1, 1]
[3, 2, 1, 1, 1, 1, 1]
[3, 1, 1, 1, 1, 1, 1, 1]


Now i would like to construct - just as an illustration, all prime numbers having $25$ digits that are using the non-zero digits from the last entry, say. We have to only generate the places for the many ones. The numbers to be generated
(and then tested for primality) are thus of the shape
$$
\underbrace{1}_{24}0\dots0
\underbrace{1}_{e}0\dots0
\underbrace{1}_{d}0\dots0
\underbrace{1}_{c}0\dots0
\underbrace{1}_{b}0\dots0
\underbrace{1}_{a}0\dots0
\underbrace{3}_1
\underbrace{1}_0
\ .
$$
The positions $a,b,c,d,e$ in the interval $[2,23]$ can be generated by asking sage for the combinations of five elements in the interval.
So the code to collect the primes is...
candidates = [ 10^24 + 31 + sum([10^a for a in pos])
               for pos in Combinations([2..23], 5) ]
solutions = [N for N in candidates if N.is_prime()]
solutions.sort()
print(f&quot;There are {len(solutions)} solutions 10...031\n&quot;
       &quot;with 25 digits with zeros and ones in the ellipsis.&quot;)
print(&quot;The 10 biggest solutions among them are:&quot;)
for N in solutions[-10:]:
    print(N)

We obtain after some seconds:
There are 1332 solutions 10...031
with 25 digits with zeros and ones in the ellipsis.
The 10 biggest solutions among them are:
1111000000000100001000031
1111000000010000000000131
1111000000010000010000031
1111000000010000100000031
1111000000010010000000031
1111000100000000000001031
1111010000000010000000031
1111010000001000000000031
1111010000100000000000031
1111011000000000000000031


Now we are in position to search all solutions with a given (decent) amount of digits. Here is the code for numbers with LENGTH=$33$ digits, say:
LENGTH = 33
solutions = []    # so far

for units_digit in [1, 3]:
    N0 = 30 + units_digit
    for main_digit in [1 .. 10 - 3 - units_digit]:
        N1 = main_digit * 10^(LENGTH-1)
        for part in Partitions(10 - 3 - units_digit - main_digit):
            print(f&quot;Candidates for {main_digit}{part}3{units_digit}&quot;)
            digits = list(set(part))
            count_dic = dict([(d, list(part).count(d)) for d in digits])
            places_ranges = [[ tuple(c)
                               for c in Combinations([2..LENGTH-2], count_dic[d])]
                             for d in digits ]
            places_list = cartesian_product(places_ranges)

            for places in places_list:
                 used_places = [place for d_places in places for place in d_places]
                 if len(set(used_places)) &lt; len(used_places):
                     continue
                 N = N0 + N1
                 for k in range(len(digits)):
                     d = digits[k]
                     d_places = places[k]
                     N += d*sum([10^a for a in d_places])

                 if N.is_prime():
                     solutions.append(N)

print(f&quot;There are {len(solutions)} solutions &quot;
      f&quot;with {LENGTH} digits&quot;)
print(&quot;The 10 biggest solutions among them are:&quot;)
for N in solutions[-10:]:
    print(N)

The above delivers the following ten biggest solutions with $33$ digits:
There are 13008 solutions with 33 digits
The 10 biggest solutions among them are:
200000100000000001000000000000033
200000000000010010000000000000033
200000000010000100000000000000033
200000010000001000000000000000033
200000000011000000000000000000033
200000100010000000000000000000033
200001000010000000000000000000033
200000001100000000000000000000033
200010100000000000000000000000033
300000000001000000000000000000033

Some other random solutions among the computed ones are:
sage: import random
sage: for _ in range(10):
....:     print(random.choice(solutions))
....:
100200000100000000100000001000031
200100000101000000000000000010031
100030000011000000000000000000031
200000000100010000100000010000031
100010000000020000010000001000031
110020000000000000000000010100031
100000000000001000001010002000031
100000000000000011000000010110031
100000000010000000021000010000031
100003000000000110000000000000031


Working with 40 instead, we get:
There are 31202 solutions with 40 digits
The 10 biggest solutions among them are:
2000000000000010000000001000000000000033
2000000100000000000010000000000000000033
2000100000000000000010000000000000000033
2000000000100000010000000000000000000033
2000000000000100100000000000000000000033
2000000000001010000000000000000000000033
2100000000000010000000000000000000000033
2000000110000000000000000000000000000033
2110000000000000000000000000000000000033
3000010000000000000000000000000000000033

One can ask for solutions with &quot;special properties&quot;, for instance all solutions with $40$ digits and three times the digit 3 are...
sage: [N for N in solutions if N.digits().count(3) == 3]
[3000000000000000000000000000030000000031,
 1000000000030000000000000000000000000033,
 1000000030000000000000000000000000000033,
 3000010000000000000000000000000000000033]


Some words on the search method, so that the code can be digested without problems. Let us consider for instance the solution
100000000000000000000000002100231

its units_digit is 1. It has position zero in the number. (Counting from right.)
its main_digit is 1. It has position $32$ in the number. (Counting from right.)
the search is done now for number of the shape $1\dots 31$ with $33$ digits. &quot;Inside&quot; the dots we have many occurrences of zero, but some few digits are not zero and have to add to $5=10-1-3-1$. For this reason we consider all partitions of $5$. At some point we consider also the partition part = [2, 2, 1]. The digits inside this part are building the list [2, 1]. (We take the set associated to the list part, then make it a list for programmatic reasons.)
it is handy to record how many digits are in the part for each encountered digit, this gives the dictionary {2: 2, 1: 1}, so two comes twice, the one once.
we are searching now for possible positions for the two digits 2, and for the one digit 1. First we consider the places for the 2. We build the Combinations([2..31], 2). An element in these combinations object is [2, 6]. (It leads to the above solution.)
Then we consider the place&quot;s&quot; for the 1. We build the Combinations([2..31], 1). An element in these combinations object is [5]. (It leads to the above solution.) We consider the cartesian product. The tuple $(\ $[2, 6]$\ ,\ $[5]$\ )$ is somewhere in the cartesian product. The used places are [2, 6, 5], and we have the full length, i.e. there is no collision of the places dedicated for the 1 digit and the 2 digits.
then we consider the corresponding number
$N=1\cdot 10^{32}+2\cdot (10^6+10^2)+1\cdot(10^5)+30 +1$ and check if it is a prime number. Yes, it is. So we append it to the list of (already found) solutions (and search further and possibly append further...)


Some words on the true mathematical problem behind the puzzle. The &quot;digits of numbers&quot; do not have a structural meaning. Methods to find infinitely many solutions are a matter of human inventive force. A more concrete question would be if there are infinitely many prime numbers in a specific infinite subset like the one mentioned in the OP:
$$
N(a) := 
10^a + 333\ .
$$
Well, this is also a complicated problem. Some cheap exclusions would be...

considering $N(a)$ modulo $7$, we have $N(a)=0$ mod seven for $a=1$ mod six. (So $N(a)$ has the factor $7$ for $a=7, 14, 21,$...
considering $N(a)$ modulo $p=31$, we have $N(a)=0$ mod $p$ for $a=3$ mod $(p-1)$. In fact even for $a=3$ mod $15$.

and so on. Here are the factorizations of the first few numbers in the sequence:
sage: for a in [3..60]:
....:     N = 10^a + 333 
....:     print('PRIME' if N.is_prime() else '     ', N, '=', N.factor())
....:

      1333 = 31 * 43
PRIME 10333 = 10333
PRIME 100333 = 100333
PRIME 1000333 = 1000333
      10000333 = 7 * 857 * 1667
      100000333 = 79 * 1265827
      1000000333 = 17 * 139 * 423191
      10000000333 = 19 * 526315807
      100000000333 = 347 * 288184439
PRIME 1000000000333 = 1000000000333
      10000000000333 = 7 * 1428571428619
      100000000000333 = 23 * 103 * 42211903757
      1000000000000333 = 74761 * 13375958053
      10000000000000333 = 2383 * 31727 * 132265613
      100000000000000333 = 29 * 3448275862068977
      1000000000000000333 = 31 * 59 * 1573679 * 347432263
      10000000000000000333 = 7 * 113 * 157 * 9049 * 8898632591
      100000000000000000333 = 167265053 * 597853515761
      1000000000000000000333 = 79 * 402859 * 31420988107753
      10000000000000000000333 = 439 * 22779043280182232347
      100000000000000000000333 = 2350339 * 42547053850529647
      1000000000000000000000333 = 43 * 11083 * 2098332035864691157
      10000000000000000000000333 = 7 * 17 * 84033613445378151260507
      100000000000000000000000333 = 244129 * 5233003 * 78276183759559
      1000000000000000000000000333 = 61 * 267139 * 13269967 * 4624481284981
      10000000000000000000000000333 = 19 * 5387 * 6550519 * 14915015630858219
      100000000000000000000000000333 = 4270499 * 23416467255934259673167
      1000000000000000000000000000333 = 4517 * 1159027 * 191010110705909287187
      10000000000000000000000000000333 = 7 * 431 * 3314550878355982764335432549
      100000000000000000000000000000333 = 10729 * 86117699364577 * 108230168748901
      1000000000000000000000000000000333 = 31 * 32258064516129032258064516129043
      10000000000000000000000000000000333 = 79 * 179 * 15443 * 45791851773167882549468491
      100000000000000000000000000000000333 = 3058170639611 * 32699287183242339516503
      1000000000000000000000000000000000333 = 23^2 * 4129417913867017 * 457778604072487381
      10000000000000000000000000000000000333 = 7 * 47 * 1579 * 19249611639085181456464115836463
      100000000000000000000000000000000000333 = 191 * 1223 * 26683 * 59000101 * 133534237 * 2036386152511
      1000000000000000000000000000000000000333 = 5183468525235985007 * 192921013242667183619
      10000000000000000000000000000000000000333 = 163 * 148250309 * 413825061582392670299853092099
      100000000000000000000000000000000000000333 = 17 * 41809 * 140695853552499954273847595437514861
      1000000000000000000000000000000000000000333 = 1049011 * 57358632811 * 16619622950403176530634173
      10000000000000000000000000000000000000000333 = 7^2 * 877 * 1699 * 2731 * 50152114413015514081748431875409
      100000000000000000000000000000000000000000333 = 1303 * 461596713173393557 * 166261952579593717268623
      1000000000000000000000000000000000000000000333 = 29 * 43 * 75344436707208992173 * 10643448330526537714943
      10000000000000000000000000000000000000000000333 = 19 * 526315789473684210526315789473684210526315807
      100000000000000000000000000000000000000000000333 = 79 * 21493 * 2140363 * 7291466811954023 * 3773753526074490611
      1000000000000000000000000000000000000000000000333 = 31 * 103 * 313185092389602254932665205136235515189476981
      10000000000000000000000000000000000000000000000333 = 7 * 669998127191476652301479 * 2132202122056322773969661
      100000000000000000000000000000000000000000000000333 = 173909 * 84555100781 * 6800457497533245880372446409581277
      1000000000000000000000000000000000000000000000000333 = 44319249889 * 863620780381437833 * 26126697381028729821509
      10000000000000000000000000000000000000000000000000333 = 467 * 919 * 12457 * 2145379 * 871866755809868463760709656795160507
PRIME 100000000000000000000000000000000000000000000000000333 = 100000000000000000000000000000000000000000000000000333
      1000000000000000000000000000000000000000000000000000333 = 2087 * 2423 * 15013 * 2594461187 * 6381048309061 * 795641525842044878663
      10000000000000000000000000000000000000000000000000000333 = 7 * 139 * 192106702699 * 923265594035567 * 57945269915555160600013037
      100000000000000000000000000000000000000000000000000000333 = 4603 * 803977 * 34682531183 * 779120470267938849650272081658642321
      1000000000000000000000000000000000000000000000000000000333 = 17 * 823 * 971 * 1658571172243 * 7209852591239633 * 6155615937870395875787
      10000000000000000000000000000000000000000000000000000000333 = 23 * 405720054989371 * 1071632036299123913967888194095887739225601
      100000000000000000000000000000000000000000000000000000000333 = 816271 * 122508333629395139604371587377231336161642395723969123
      1000000000000000000000000000000000000000000000000000000000333 = 79 * 35423 * 357344884625843825276933351962913318493991424437458749

The &quot;easy prime&quot; mentioned last as PRIME in the above has $54$ digits. The next ones are for $223$, then $232$ digits...
sage: for a in [3..2000]: 
....:     N = 10^a + 333 
....:     if not N.is_prime(): 
....:         continue 
....:     print(f&quot;The number 10...0333 with {len(str(N))} digits is a PRIME&quot;) 
....: 
The number 10...0333 with 5 digits is a PRIME
The number 10...0333 with 6 digits is a PRIME
The number 10...0333 with 7 digits is a PRIME
The number 10...0333 with 13 digits is a PRIME
The number 10...0333 with 54 digits is a PRIME
The number 10...0333 with 223 digits is a PRIME
The number 10...0333 with 232 digits is a PRIME
The number 10...0333 with 417 digits is a PRIME

We have a hard question behind the puzzle, i have no method to attack it.
But experimenting with such questions may be a good idea. Didactically, or as a programming exercise. Even &quot;isolated questions&quot; that are not quickly answered by &quot;cheap tricks&quot; or computer power can be hard. For instance, letting the above code go a little further, and printing the value of $a$ shows that the computer does not find an &quot;easy factor&quot; in case of $a=2111$. The isolated question is then:
$$
\text{Is }
10^{2111}+333
\text{ a prime?}
$$
(I will maybe try to use elliptic curves to have some more insight for this &quot;special $a$&quot; and for the next hard ones...)