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) "iter tools". 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 "decent size", 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"There are {len(solutions)} solutions 10...031\n"
"with 25 digits with zeros and ones in the ellipsis.")
print("The 10 biggest solutions among them are:")
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"Candidates for {main_digit}{part}3{units_digit}")
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)) < 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"There are {len(solutions)} solutions "
f"with {LENGTH} digits")
print("The 10 biggest solutions among them are:")
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 "special properties", 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. "Inside" 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"s" 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 "digits of numbers" 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 "easy prime" 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"The number 10...0333 with {len(str(N))} digits is a PRIME")
....:
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 "isolated questions" that are not quickly answered by "cheap tricks" 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 "easy factor" 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 "special $a$" and for the next hard ones...)