Второе — что в базу можно брать не все простые, а только примерно половину. Потому что если простое число p делит a_k=(m+k)^2-n, то
n=(m+k)^2 mod p,
поэтому n является квадратичным вычетом по модулю p.

Ответ на вопрос "является ли n квадратичным вычетом по модулю p" называется символом Лежандра.
И его вполне можно быстро найти. Вручную я, скорее всего, воспользовался бы мультипликативностью и квадратичным законом взаимности. Но для программирования, пожалуй, проще всего сказать, что символ Лежандра a над p равен a^{(p-1)/2} (ибо малая теорема Ферма). А k-ю степень мы вычисляем за log(k) умножений, ибо череда возведений в квадраты.

Да, в качестве рекламы — про квадратичный закон взаимности (связывающий символы Лежандра p по модулю q и q по модулю p) мне очень нравится рассуждение отсюда — https://users.mccme.ru/smirnoff/papers/qrecip01.pdf

Но так или иначе, символ Лежандра быстро ищется — так что можно быстро понять, какие простые у нас не могут появиться в принципе. Собственно, именно поэтому повторение 2, 19, 37, 47 в примере выше не столь удивительно (хотя нам и повезло — уже три вектора в четырёхмерном пространстве оказались линейно зависимыми — но куча других простых в первых десятках тут не появилась вообще).

Ну и последнее — это выбор размера базы. Возьмём слишком маленькую — и слишком мала будет доля гладких чисел, придётся перебирать много разных k. Скажем, много ли чисел вида
2^{p_1}*3^{p_2}*5^{p_3}
среди первого миллиона?

Возьмём слишком большую базу — и нужно будет много гладких чисел, опять-таки, придётся много перебирать.

Но на формулировке этой задачи оптимизации я, пожалуй, хочу в этом рассказе на сегодня остановиться.

Во-первых, моя работа над ошибками: то, что мы разобрали позавчера, это ещё не весь метод квадратичного решета, а его начальная часть. Её называют алгоритмом Диксона — обычно ссылаясь на его работу 1981 года (http://pages.cs.wisc.edu/~cs812-1/dixon.pdf ), где Диксон рассматривает квадраты случайно выбираемых чисел по модулю N и строго оценивает время работы такого алгоритма.

Если же говорить более аккуратно о том, кто что придумал, то нужно сказать, что сам придумавший квадратичное решето Carl Pomerance в своей "Tale of two fractals" (https://www.ams.org/notices/199612/pomerance.pdf ) ссылается на Крайчика в 1920-е.
По соседству там был подход с получением небольших квадратов по модулю N через цепные дроби (с идеей, что для подходящих дробей P/Q к корню из N величины P^2-N Q^2 будут не очень большими), сначала в 1930-ых в работах Лемера и Пауэрса,
https://projecteuclid.org/download/pdf_1/euclid.bams/1183495051 , а потом в работах Morrison & Brillhart (https://www.ams.org/journals/mcom/1975-29-129/S0025-5718-1975-0371800-5/S0025-5718-1975-0371800-5.pdf ).

И давайте я тут дам слово Дональду Кнуту ("Искусство программирования", том 2, раздел 4.5.4):

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

Представьте себе, что вы раскладываете на множители — и делите число M=1987645202877 на 197. Не делится. Переходите к 199. Опять не делится. Переходите дальше, и тут кто-то, пришедший понаблюдать за этим процессом, выдаёт — "а зачем вы делите на числа, на которые M не делится? Делить M надо на те числа, на которые делится!"
Скорее всего, вы захотите такому советчику много что сказать. Так вот, удивительная вещь состоит в том, что ровно так и работает квадратичное решето — за счёт того, что мы знаем, какие числа Q(x)=x^2-N на какие простые из базы будут делится!

А именно — вот у нас есть простое число p из базы. По модулю p наше N является квадратичным вычетом (иначе бы мы p из базы выкинули) — и мы можем (просто перебором, это всё равно не самая трудозатратная часть) найти корни a и (-a) из N. Но мы знаем, что Q(x) и Q(y) сравнимы, если x и y сравнимы по модулю p. Значит, на p будут делиться Q(x) при x вида a+n*p и (-a)+n*p, и только при таких!
То есть, к примеру, для простых порядка сотни тысяч или миллиона мы знаем, какие Q(x) нужно на них делить, и это должно дать ускорение в эту самую сотню тысяч-миллион раз.

Так что глобальный алгоритм просеивания (действительно напоминающий решето Эратосфена) устроен следующим образом. Вот мы будем смотреть на числа вида
(m+k)^2-N,
где m это целая часть корня из N, а k пробегает какой-то интервал, и хотим из них отобрать гладкие.

Давайте считать произведение всех степеней простых из базы, на которые такие числа делятся. А именно — заведём массив W для логарифма такого произведения ("видишь длинное произведение — прологарифмируй!"); изначально заполним его нулями.

После чего для каждого p в цикле пробежим по k, для которых m+k имеют один из двух нужных остатков по модулю p, и добавим ln p к W[k]. И, конечно, это будет не цикл по k, а движение с шагом в p.
После чего обработаем делящиеся на p^2 числа (движение с шагом p^2, добавим ещё раз ln p), на p^3, и так далее.

В результате такого цикла мы для каждого k найдём, на какое максимальное гладкое число Q(m+k)=(m+k)^2-N делится. И, сравнив его с логарифмом Q(m+k), выясним, гладкое это число или нет.

И поскольку мы, вместо того, чтобы каждое число пытаться делить — добавляли ln p только туда, куда надо, алгоритм действительно очень сильно ускорился.

Давайте я процитирую три абзаца из уже упомянутой "A Tale of Two Sieves" Карла Померанса:

Кстати, зашифрованное в 1977 году создателями RSA сообщение "The Magic Words are Squeamish Ossifrage" (https://en.wikipedia.org/wiki/The_Magic_Words_are_Squeamish_Ossifrage ) было расшифровано как раз с использованием (чуть более продвинутой и приспособленной к распараллеливанию) версии квадратичного решета — позволившей разложить на множители соответствующее 129-значное число RSA-129 (см. https://en.wikipedia.org/wiki/RSA_numbers#RSA-129 )

Ну и на этом, наверное, я этот рассказ [на время?] завершу. Тут ещё осталось много чего нерасказанного (включая, например, распараллеливание вычислений), но вообще-то я собирался сделать всего лишь короткий перерыв между двумя рассказами о результатах Фюрстенберга и Маргулиса. :)

Давайте я начну рассказывать об одном из результатов Маргулиса — о доказательстве гипотезы Оппенгейма (https://en.wikipedia.org/wiki/Oppenheim_conjecture ). И для начала посмотрим на саму эту гипотезу, посвящённую значениям квадратичной формы от нескольких переменных в целых точках.

Вообще, значения квадратичной формы нескольких переменных в целых точках даже для квадратичной формы с целыми коэффициентами это большая, красивая и нетривиальная наука. Сразу можно вспомнить теорему Лагранжа, утверждающую, что любое натуральное число представляется как сумма четырёх квадратов. Иными словами, значения на Z^4 квадратичной формы B(x,y,z,t)=x^2+y^2+z^2+t^2
это все неотрицательные целые числа.

Можно вспомнить критерий того, когда натуральное число n представляется как сумма двух квадратов: простые вида 4k+3 должны входить в разложение n на простые множители только в чётных степенях.