Алгоритм Шенкса-Тонелли

Есть две последовательности чисел:
x, x + 1, x + 2, ...
и
x2(mod n), (x + 1)2(mod n), (x + 2)2(mod n), ... где n — любое: n ⩽ x2
Необходимо найти первый элемент из второй последовательности, который делится на простое p
Простым перебором от начала последовательности это делать долго при большой последовательности и больших p
На сколько я понял, тут нужно использовать алгоритм Шенкса-Тонелли. Правда, я не понял, как его тут применить.
Собственно вопрос: как применить к данной задаче данный алгоритм?
Или, если вы знаете способ решить задачу еще быстрее — как это сделать?

Ответ

В общем, не так много где этот вопрос расписан, особенно, в русскоязычной литературе.
Я бы тоже не стал описывать ответ если бы не одно но: данный вопрос актуален при реализации алгоритма квадратичного решета на этапе просеивания. Решение данного сравнения позволит найти элементы, делящиеся на числа из факторной базы без перебора, что значительно ускорит процесс просеивания при больших простых из факторной базы.
Собственно, ответ на мой вопрос:
Сперва теория:
Рассмотрим конечное поле GFp (p > 2), и элемент a, являющийся квадратичным вычетом по модулю p. Требуется найти x такое, что x2 ≡ n (mod p).
Представим число p − 1 в виде p − 1 = 2rs, где s – нечетно. Заметим, что поскольку p − 1 – четно, то r ⩾ 1.
Пусть z – квадратичный невычет по модулю p. Вычислим y ≡ zs (mod p). Поскольку порядок любого элемента является делителем числа 2rs, то порядок y является делителем 2r, откуда y2r ≡ 1 (mod p). Можно также показать, что y2r−1 ≡ −1 (mod p), т.е. порядок элемента y равен в точности 2r
Вычислим далее элементы:
λ0 ≡ as (mod p) ω0 ≡ a(s+1)/2 (mod p)
Заметим, что ω02 ≡ aλ0 (mod p) и x2 ≡ a (mod p) ⟹ x2s ≡ as ≡ λ0 (mod p).
Поскольку порядок элемента xs является делителем 2r, то порядок λ0 является делителем 2r−1
Идея метода Шенкса–Тоннелли состоит в построении последовательности пар чисел (ωi, λi), удовлетворяющих условию:
ωi2 ≡ aλi (mod p) (i = 1, 2, 3, ...)
причем порядок λi+1 является собственным делителем порядка λi, до тех пор, пока порядок очередного λi не окажется равным 0. Тогда для найденного i выполнятся условия λi = 1 и ωi2 ≡ a (mod p) откуда x = ωi является искомым корнем.
Поскольку исходные значения (λ0, ω0) уже определены, то осталось просто описать формулы для вычисления значений (λi+1, ωi+1):
λi+1 ≡ λi · a · y2r−m(mod p) ωi+1 ≡ ωi+1 · y2r−m−1(mod p)
где m - порядок элемента λi. Так же, корнями будут все числа, кратные p.
Реализация же алгоритма крайне проста:
int Shenks_Tonelli(int p,long long n){ n=n%p; int s=p-1,r=0; //получаем разложение p-1 while(s%2==0){ s/=2; r++; } //начальные значения: λ и ω int l = PowMod(n,s,p); int w = PowMod(n,(s+1)/2,p); //находим порядок λ int mod=l,m=0; while(mod!=1){ mod=MulMod(mod,mod,p); m++; } //находим квадратичный невычет int z = quadratic nonresidue(p); //находим коэф-ты, на которые будем умножать int yd_l=PowMod(PowMod(z,s,p),pow(2,r-m),p); int yd_w=PowMod(PowMod(z,s,p),pow(2,r-m-1),p); //находим корень while(l!=1){ l=MulMod(l,yd_l,p); w=MulMod(w,yd_w,p); } return w; }
Правда, этап просеивания можно еще ускорить, если находить корни сравнения x2 ≡ n (mod pβ), и делить числа, отстоящие от них на pβ, постепенно понижая β. Корни такого сравнения находятся по индукции. Подробнее об этом здесь