Поэтому более правильно определять диофантовость так: подмножество A в Z_+^k называется диофантовым, если существует многочлен P(a_1,...,a_k,x_1,...,x_n), такой, что
(a_1,...,a_k) из A <=> существуют (x_1,...,x_n) из Z_+^n, для которых
P(a_1,...,a_k,x_1,...,x_n)=0.

(а систему из нескольких диофантовых уравнений можно свернуть в одно взятием суммы их квадратов)

Ну и — на идейном уровне теперь можно понять, как именно должен работать тот порождающий простые числа многочлен: остальные переменные это "сертификат простоты" для искомого k+2 (спасибо А. Горденко за рисунок!).

Более того, если чуть-чуть расширить класс многочленов до экспоненциальных многочленов, разрешив ещё возводить одну переменную в степень другой, то такой многочлен строится совсем "вручную".

Для начала вообще представим себе, что мы имеем право использовать ещё операцию "факториал". Как тогда можно задать простоту числа p?
Во-первых, нужно попросить, чтобы p было не меньше 2. Ну, тут мы уже видели, как это сделать: одно из уравнений системы будет
p=k+2.
А как проверить, что p не делится ни на что, кроме 1 и себя? Цикл использовать нельзя — но можно вместо него воспользоваться факториалом. А именно — давайте посмотрим на (p-1)! ; число p простое тогда и только тогда, когда оно взаимно просто с (p-1)!, а это проверяется алгоритмом Евклида:
a*(p-1)!-b*p=1.
На самом деле, теорема Валлиса утверждает, что простота p равносильна сравнимости (p-1)! с (-1) по модулю p, так что одну переменную можно сэкономить, и написать просто
(p-1)!+1-b*p=0.

И пара (k,b) тут — "сертификат простоты" числа p, который есть для всех простых p и только для них.

А как быть, если наш бюрократ не понимает факториалы? Нужно "соотношение факториальности" как-нибудь выразить:
p=k+2
c+1-b*p=0
//а тут какие-нибудь уравнения, задающие c=(p-1)! //

А через что мы можем факториал выразить, где он появляется? Первое, что приходит в голову — биномиальные коэффициенты. И действительно, если нам нужно r!, а число M очень большое, то посмотрим на
C_M^r = M(M-1)...(M-r+1) / r!

Тогда отношение
M^r / (C_M^r)
чуть-чуть больше r! ; поэтому если M^r на C_M^r поделить с остатком, то при всех достаточно больших M мы как раз r! и получим.

А деление с остатком мы выразить можем: чтобы q было неполным частным от деления A на B, нужно, чтобы выполнялись неравенства
A>=B*q
и
A<B*(q+1)
А нестрогое неравенство — это равенство, если добавить к меньшей части ещё одну свободную переменную (а для строгого нужно ещё 1 добавить). Итого соотношение q=r! выражается системой на неотрицательные целые переменные

// M достаточно большое относительно r, //
M^r = q*(C_M^r) + x_1,
(q+1)*(C_M^r) = M^r + 1 + x_2

Чтобы добиться исходной цели — задать простоту с помощью экспоненциальных многочленом — остаётся выразить соотношение c=C_M^r. Но биномиальный коэффициент это коэффициент в биноме Ньютона. Так давайте возьмём ну очень большое (относительно M и r) число N и посмотрим на (N+1)^M.

Если взять неполное частное q' при делении на N^r, а потом остаток от деления q' на N — то как раз C_M^r и будет. А и неполное частное, и остаток мы выражать уже умеем!

И если это всё сделать — мы получим систему экспоненциально-полиномиальных уравнений, таких, что для p найдутся все остальные (неотрицательные целые) переменные тогда и только тогда, когда p простое. Ура!

Собственно, мой рассказ выше это пересказ [части] статьи Ю. В. Матиясевича в "Кванте" в 1975 году — http://kvant.mccme.ru/1975/05/formuly_dlya_prostyh_chisel.htm

Возвращаясь к диофантовости всех перечислимых множеств — собственно говоря, идея "сертификата" работает и тут.

А именно — представим себе, что какая-то машина долго-долго работала и наконец закончила работу. Давайте соберём все состояния её "регистров" во все моменты времени в один "протокол" ("в момент времени t=1 состояние такое-то; в момент t=2 состояние такое-то; ..."). Так вот — машина заканчивает работу за конечное время тогда и только тогда, когда существует конечный "протокол", в котором в начальный момент времени стоит её начальное состояние, корректно выполнен переход от каждого момента к следующему, и конечное состояние соответствует команде "стоп". Этот "протокол" мы и захотим добавить в дополнительные переменные.
Правда, мы не знаем, сколько машина будет работать — но не страшно, просто пусть у нас будет система счисления с большим-большим основанием N, а последовательные состояния машины — последовательными "цифрами" в записи одной из переменных в этой системе счисления. Собственно, мы уже применили почти этот подход к выражению биномиальных коэффициентов — только тогда мы "доставали" из большого числа (N+1)^M одну его цифру, а тут придётся научиться работать со всеми цифрами одновременно, каким-то образом одновременно обеспечивая, что переход от каждого момента времени t к следующему t+1 корректен.
В качестве одного из шагов возникает теорема Куммера о делимости биномиальных коэффициентов на простые:

И если вот этим путём (на который я, конечно, пока только намекнул) пройти — то получается
Теорема (Davis, Putnam, Robinson, 1961; https://doi.org/10.2307/1970289 )
Любое перечислимое множество экспоненциально-диофантово.

Но что делать с экспонентой? Как выразить хоть что-то экспоненциальное, если у нас в руках только многочлены? Ну хоть что-нибудь экспоненциально растущее — скажем, уже было известно, что если есть какое-то отношение, которое растёт быстрее, чем u^k при любом k, но медленнее, чем u^u, то этого уже хватает. Но как это вытащить из полиномов? И полноте, возможно ли это? [Это сейчас мы знаем, что да — а вот обзор в Mathematical Reviews на их работу в этом смысле не горит энтузиазмом...]

И — последний удар!

Теорема (Ю.В. Матиясевич, 1970; см. http://mi.mathnet.ru/dan35274 )
Отношение "v — 2u-е число Фибоначчи" диофантово.

Два кусочка из статьи Матиясевича:

И собственно явное задание —

Собственно, эта статья замечательна тем, что её можно спокойно разобрать сильному школьнику.

Вот основные идеи:
а) Соотношение x^2-xy-y^2=1 на неотрицательные x и y равносильно тому, что x и y это два последовательных числа Фибоначчи, x=F_{2m+1}, y=F_{2m}. А если бы в правой части было (-1), то было бы x=F_{2m}, y=F_{2m-1}.
Утверждение само по себе симпатичное — и мгновенно доказывающееся по индукции (достаточно написать, что x=y+z).
И тут можно посмотреть на равенства (2) и (3) из кусочка выше — гарантирующие, что l, z, g и h это числа Фибоначчи.

(соответствующая пара лемм)

б) Делимость чисел Фибоначчи. Во-первых, F_m делится на F_n тогда и только тогда, когда m делится на n. Что несложно увидеть, приведя числа Фибоначчи по модулю F_n: будет
номер — 0, 1, 2, ..., n-1, n, n+1,...
число — 0, 1, 1, ..., * , 0, * , ...,
и после n будет то же, что и после нуля, только умноженное на F_{n-1} (взаимно простое с F_n).
Но вот — и это ключевой шаг — если F_m делится не просто на F_n, а на F_n^2, то вот это равносильно тому, что номер m делится на nF_n.