Как раз все числа от P(1)=17 до P(16)=15*16+17=257 — простые, а первыми составными оказываются P(17)=17^2 и P(18)=17*19.

Ещё таким будет p=11, это даже вполне несложно проверить вручную —

(И p=2, 3 и 5, кстати, но тут картинка получается не очень впечатляющей: для p=5 там всего 4 значения, для p=3 два, для p=2 вообще одно. Но я их тут всё-таки упомяну — потому что они укладываются в общий принцип.)

Так вот — а зачем вообще (кроме общей красоты скатерти Улама) я так внимательно смотрю на эти простые числа? Дело в том, что отсюда открывается дорога к однозначности разложения на простые в кольцах целых квадратичных расширений.

Классический вопрос о том, какие числа n представляются в виде суммы двух квадратов, связан с гауссовыми целыми числами — числами вида a+bi, составляющими кольцо Z[i].
Ибо (a+bi)(a-bi)=a^2+b^2, и условие про то, что все простые вида 4k+3 в разложении n должны быть в чётной степени, как раз связано с тем, что они при переходе в Z[i] остаются простыми, а 2 и простые вида 4k+1 — нет, "разваливаются" (например, 2=(1+i)(1-i), 5=(1+2i)(1-2i)).
И при работе со всем этим ключевую роль играет то, что в гауссовых целых Z[i] есть алгоритм Евклида и, как следствие, однозначность разложения на их простые.

Так будет не всегда — и стандартный пример состоит в том, что если вместо i=\sqrt{-1} к целым числам добавить \sqrt{-5}, то однозначность разложения на множители нарушится:
(1+\sqrt{-5})(1-\sqrt{-5}) = 6 = 2*3.

Кстати, разложение на простые пытались применить и к теореме Ферма: пытаясь доказать, что у уравнения
x^p+y^p=z^p
нет нетривиальных решений — очень естественно добавить к целым числам корень p-й степени из 1,
\zeta=\exp(2πi/p),
разложить левую часть как произведение (x+\zeta^j y) по j=0,1...,p-1, и попробовать что-нибудь вытащить, сопоставляя разложения на множители в правой и в левой части.
И тут нужно начинать говорить слова "идеалы", ключевым именем тут будет Куммер, и я оставлю ссылку вот на этот короткий текст Кита Конрада —
https://kconrad.math.uconn.edu/blurbs/gradnumthy/fltreg.pdf
(Получается так доказать не для всех простых — так что пришлось ждать ещё 150 лет до доказательства Уайлса — но для многих.)

[+Картинка выше: несколько первых абзацев из текста Конрада.]

Кстати, ещё можно спрашивать не какие числа представимы как a^2+b^2 (т.е. какие квадраты длин бывают для отрезков между точками квадратной решётки), а какие числа представимы как a^2+ab+b^2 (т.е. какие квадраты длин бывают для отрезков между точками треугольной решётки).

И тут на помощь придут числа Эйзенштейна Z[w], где w=\exp(2πi/3)= -1/2 + i sqrt{3}/2.
В которых тоже есть алгоритм Евклида — и, значит, однозначность разложения на (их) простые!

Так вот — а как обстоят дела с однозначностью разложения на простые в других квадратичных расширениях? Если добавлять к целым числам только \sqrt{-d}, то плохо. Но если мы рассмотрим квадратичное расширение Q[\sqrt{-d}] поля рациональных чисел — то есть множество чисел вида a+b\sqrt{-d}, где a и b рациональные — то там более естественно посмотреть на алгебраические целые числа. А при d=4k+3 это не только такие комбинации с целыми a и b, но и с одновременно полуцелыми: число z=(-1+\sqrt{-d})/2 является корнем уравнения z^2+z+ (d+1)/4 =0
И где-то мы такой многочлен уже видели, правда, вместо (d+1)/4 стояло p...

А если более аккуратно — то давайте посмотрим, чему равна норма числа (z+n) — то есть произведение (n+z) на его сопряжённый:
P(n)=(z+n) (\bar{z}+n) = (n- 1/2)^2 + (\sqrt{d}/2)^2 = n^2 -n +1/4 + d/4 = n^2 - n +p,
если p=(d+1)/4.
Поэтому если при малых n это произведение не простое, P(n)=ab, то в кольце целых O(\sqrt{-d}) у нас появится такая же неоднозначность с разложением на простые, что и раньше в Z[\sqrt{-5}] — потому что будет
ab= P(n) = (z+n) (\bar{z}+n)
(и из-за того, что норма z+n ещё маленькая, ни на что дальше сомножители в правой части раскладываться не могут).

Оказывается (и в детали я тут вдаваться не буду) — что на самом деле это критерий: все числа n^2-n+p для n=1,2,...,p-1 простые тогда и только тогда, когда в кольце целых O(\sqrt{1-4p}) есть однозначность разложения на простые. А число d, для которого в O(\sqrt{-d}) есть однозначность разложения, называется числом Хегнера ( Heegner number ) ; оказывается, что их таких ровно 9:
1, 2, 3, 7, 11, 19, 43, 67, 163.

А ещё три последних числа в этом списке — совершенно случайно, чистейшей воды совпадение! — появлялись в той байке, которую я рассказывал почти сразу после скатерти Улама: в рассказе про почти целые числа.

В той затравке предлагалось посмотреть на три числа вида e^{π\sqrt{d}), где d=163, 67 или 43. И во всех трёх случаях оказывалось, что после запятой заметное количество (а в первом случае даже просто много!) девяток.

Давайте я покажу кадр из видео Ричарда Борчердса, https://youtu.be/a9k_QmZbwX8?t=197 .

Думаю, при взгляде на всё написанное становится понятно, что "таких совпадений не бывает".
И, как и в случае с числом (1+\sqrt{2})^100, механизм "почти целости" тут — что есть по-настоящему целое число, которое складывается из "почти целого" — и чего-то очень маленького.

(В скобках — каюсь, я довольно долго считал, что e^{π\sqrt{163}) это такое случайное совпадение, ну вот написали что-то, и попали близко к целому; и даже то, что девяток после запятой ну слишком много, меня в те годы не настораживало.)

Так вот — есть такая функция, j-инвариант. Это "хороший" инвариант эллиптической кривой — то есть фактора C/Г комплексной плоскости по решётке Г; иными словами, "хорошая" функция на множестве решёток на плоскости, рассматриваемых с точностью до умножения на комплексные числа (потому что умножение на λ переводит решётку C/Г в решётку C/(λГ)).

Её можно строить разными способами, но я сначала просто поверю в то, что она есть, и скажу про пару её свойств.
В частности — можно в решётке выбрать положительно ориентированный базис, и поделить на первый его элемент; тогда решётка примет вид <1, τ>, где τ — точка из верхней полуплоскости {Im τ>0}. Соответственно, тогда нужно строить функцию j(τ) в верхней полуплоскости — но инвариантную относительно замены базиса. А это не только замена τ-> τ+1 (что очевидно), но и вообще все действие SL(2,Z) перевыборами базиса:
τ-> (aτ+b)/(cτ+d),
где a,b,c,d целые и ad-bc=1.

Но уж раз j(τ) инвариантна относительно замены τ->τ+1, то можно рассматривать j как функцию на факторе верхней полуплоскости по этому сдвигу. А это — проколотый диск с координатой q=exp(2πi*τ); проколотый, потому что чем больше мнимая часть τ, тем меньше модулю q, но нулевым он никогда не будет.

Так вот: если взять j-инвариант для решётки O(\sqrt{-d}), и в O(\sqrt{-d}) есть однозначное разложение на множители — то получающееся число (по нетривиальным причинам) в точности целое!

Но — функцию j(q) в координате q можно в окрестности q=0 (отвечающему стремящейся к бесконечности мнимой части τ — то есть очень "неравномерным", "вертикально-вытянутым" решёткам) раскладывать в ряд. И как оказывается, ряд для j(q) начинается с 1/q (чем более решётка неравномерная, тем больше инвариант) — а дальше уже идут неотрицательные степени:
j(q)=1/q + 744 + 196844*q +... .
Для решётки целых O(\sqrt{-d}) второй вектор это τ=-1/2 + i \sqrt{d}/2, и соответственно
q= exp(2πi τ) = exp(-πi - π\sqrt{d}) = - 1/ exp(π\sqrt{d}) —
как раз почти те числа, что нам нужны!

И, соответственно, для каждого d=43, 67, 163 целое (отрицательное) число j(q) это огромная величина
1/q=-exp(π\sqrt{d}),
отличающаяся от нашей знаком, к которой добавляется 744, и ряд по положительным степеням q — которые, наоборот, уже очень маленькие. Вот отсюда и получаются
...743,9999...,
которые мы видим до и после запятой — и поэтому девяток у exp(π\sqrt{163}) больше, чем у exp(π\sqrt{67}).

Ещё один кадр из того же видео Борчердса: https://youtu.be/a9k_QmZbwX8?t=623 ; кстати — не так часто у фидлсовского медалиста (и не только — https://en.wikipedia.org/wiki/Richard_Borcherds ) есть свой канал на YouTube :)

Вообще-то говоря, рассказ выше немного "повисает в воздухе", потому что я продекларировал функцию — j-инвариант, которую совсем не построил (и потому про которую пока в принципе непонятно, как что-либо проверять). Поэтому давайте я про неё скажу пару слов.

Вот если у нас есть решётка Г на комплексной плоскости и целое число k, то можно попробовать посчитать сумму обратных k-ых степеней ненулевых элементов решётки. И обозначить эту сумму через G_k(Г).
При нечётном k мы получим ноль — просто всё сократится. При k=2 сумма не будет абсолютно сходящейся, так что давайте её просто не рассматривать. А вот при чётных k>2 получается величина, хорошо зависящая от решётки. Правда, не инвариантная при перемасштабировании:
G_k(λГ) = λ^{-k} G_k(Г).

Ну так и давайте возьмём G_4^3(Г) и G_6^2(Г). Оба при умножении решётки на λ будут умножаться на λ^{-12} — а значит, их отношение это инвариант!

Это ещё не j-инвариант, но j-инвариант от этого отличается всего лишь на дробно-линейное преобразование: например, мы хотим, чтобы инвариант стремился к бесконечности, когда решётка становится всё более вырожденной, и поэтому в знаменателе будет не чистое G_6^2 или G_4^3, а некоторая их линейная комбинация (и даже понятно, какая — когда τ убегает на бесконечность, G_4 стремится к удвоенному \zeta(4), а G_6 к удвоенному \zeta(6), ну и остаётся подобрать коэффициенты, чтобы знаменатель обратился в ноль).