Если бы такой способ проверки существовал — например, его применение позволяло бы для любого фиксированного n>2 доказать теорему Ферма (потому что это вопрос о разрешимости полиномиального уравнения x^n+y^n-z^n=0). Было бы возможно за конечное время (и без дополнительных идей) как-то проверять, представимо ли число в виде суммы трёх кубов целых чисел (которые могут быть и отрицательными, так что это не конечная проверка), например, представимы ли 33 (https://www.quantamagazine.org/sum-of-three-cubes-problem-solved-for-stubborn-number-33-20190326 , https://people.maths.bris.ac.uk/~maarb/papers/cubesv1.pdf ) или 42 (https://www.youtube.com/watch?v=ASoz_NuIvP0 ).
Позволять проверить существование решения в натуральных числах у уравнения
x/(y+z) + y/(z+x) + z/(x+y) = 4
(если вы не видели, то это совершенно прекрасная история — см. https://news.1rj.ru/str/cme_channel/198 + https://habr.com/ru/post/335248/ + https://www.quora.com/How-do-you-find-the-integer-solutions-to-frac-x-y+z-+-frac-y-z+x-+-frac-z-x+y-4/answer/Alon-Amit?share=1 ).
И так далее.

Так вот — теорема Дэвиса, Матиясевича, Патнема и Робинсон даёт отрицательное решение проблемы Гильберта — утверждая, что такого алгоритма в принципе существовать не может.

Потому что перечислимые множества это совершенно не то же самое, что разрешимые (принадлежность которым можно проверить за конечное число операций). И первый пример тут — множество кодов программ, которые когда-нибудь останавливаются, перечислимо (ибо можно запустить симулятор + перебор кодов + параллелизацию), но не разрешимо (ибо теорема Тьюринга и парадокс лжеца).

Ну и прежде, чем обсуждать, почему это так, давайте ещё чуть-чуть посмотрим на эту формулу. И заметим, что эта формула для простых чисел... для удобства пользователя разложена на множители!

Но если присмотреться, то видно, что вторая скобка это 1 минус сумма квадратов. Поэтому, если значение P положительное, то вторая скобка должна быть равна 1, а каждый квадрат внутри неё должен обращаться в ноль.
Поэтому эту теорему можно переформулировать так: число k+2 простое тогда и только тогда, когда существуют (все остальные переменные), для которых выполняется (система полиномиальных уравнений):

Поэтому более правильно определять диофантовость так: подмножество 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-е число Фибоначчи" диофантово.

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

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