Для начала — давайте мысленно поделим последовательность s_n на 2 (благо, что по нечётному модулю M это проблемы не создаёт; ну или можно сначала делить, она останется целой, а потом уже приводить по модулю).
Для новой последовательности t_n реккурентное соотношение это t_{n+1}=2 t_n^2-1.
И это уже начинает что-то напоминать: P(t)=2 t^2-1 это формула для косинуса двойного угла (она же многочлен Чебышева степени 2).

А косинус двойного угла, если перейти в комплексные числа и рассмотреть там z=cos a + i sin a, это просто возведение в квадрат. А череда из k возведений в квадрат это возведение в степень 2^k. И это уже начинает с числами Мерсенна (M=2^p-1) иметь что-то общее — ещё не формально, но на уровне "это явно кусочек от той же головоломки".

Пытаясь сделать что-то чуть более формально, мы замечаем, что всё начинается с t_1=2, а косинус, равный двум, это немного странно. Но давайте всё равно посмотрим, чему бы равнялось соответствующее комплексное z?
Если z=e^{ia} = cos a + i sin a, то cos a= (z+1/z)/2.
Соответственно, уравнение (z+1/z)/2=2 даёт корни z=2+\sqrt{3} и z=2-\sqrt{3}.

И кстати, сразу видно, что удобнее работать не с косинусом, а с удвоенным косинусом, так нет двойки в знаменателе.
Поэтому с этого момента мы вернёмся к исходной последовательности s_n.
Так вот — для уже вовсе не комплексного z=2+\sqrt{3} мы получаем z+1/z =4. И поэтому члены нашей исходной последовательности s_n до приведения по модулю M получаются как
s_n=(2+\sqrt{3})^{2^{n-1}} + (2-\sqrt{3})^{2^{n-1}}

Осталось понять, при чём же тут всё-таки числа Мерсенна, как будет играть степень двойки, и что нам делать с корнем из 3 при работе с вычетами.

На самом деле — ассоциация с комплексными числами тут очень по делу. Комплексные числа получаются из вещественных, если к ним добавить корень уравнения x^2+1=0, у которого в самом поле вещественных чисел корня нет.
Этот корень мы обозначаем через i; при этом в комплексных числах у этого уравнения есть два равноправных (с точки зрения вещественных чисел) корня, i и -i, и переводя один из них в другой, мы получаем автоморфизм комплексных чисел — комплексное сопряжение.

Точно так же вместо вещественных чисел можно взять, например, поле из 7 элементов — вычеты по модулю 7.
В нём у уравнения x^2+1=0 тоже нет корней (упражнение: проверьте это!). И добавив такой корень — и всё, что добавится из-за этого — мы получим поле из 7^2=49 элементов, элементы которого это комбинации вида a+bx (где x — наш корень). Классический вопрос « для каких простых p по модулю p есть корень из (-1) » имеет столь же классический ответ « для двойки и простых вида 4k+1 ».

Так вот — мы для нашей последовательности хотим рассмотреть числа, включающие в свою запись «корень из 3», и работать с ними «по модулю M». Вспомогательная лемма, которую мы сейчас принимаем на веру: для любого M=2^p-1, где простое p>2, из 3 по модулю M корень извлечь нельзя.

Следствие. Если M простое, то мы имеем право рассмотреть не только поле вычетов по модулю M, но и поле из M^2 элементов, получающееся добавлением x=\sqrt{3}. В котором есть «сопряжение» — замена знака перед \sqrt{3}.

Второе — важное само по себе — утверждение, это что мультипликативная группа конечного поля циклична. То есть что все ненулевые его элементы есть степень некоторого одного, поэтому операция «умножить» превращается ( если его знать) в операцию «сложить степени».

А давайте посмотрим, из скольки элементов состоит эта мультипликативная группа для нашего нового поля — то есть по какому модулю будет сложение степеней?

Всех элементов в поле M^2, без нуля
M^2-1=(M-1)(M+1) = 2^p (2^p-2).
И у нас внезапно образовалась группа, порядок которой делится на очень большую степень двойки (на 2^{p+1}, если быть точным).

Давайте закончим доказательство теста Люка-Лемера. Для начала посмотрим, какое у него условие «успешного окончания». Это s_{p-1}=0 — «вещественная (точнее, не содержащая \sqrt{3}) компонента» у числа z^{2^{p-2}} равна нулю.
Если бы речь шла о косинусе, это бы значило, что 2 cos a=0, и соответствующие комплексные числа были бы i и -i; а их квадрат равнялся бы (-1).
Собственно, s_{p-1}=0 равносильно тому, что s_{p}=-2, то есть z^{2^{p-1}}=-1.

Не то, чтобы это было обязательно — но вот пример для p=5 (которое приводит к действительно простому M=31):

Жирным тут выделены значения t_j — "вещественные части", половинки от последовательности s_j.

Итак, от z=2+\sqrt{3} (соответствующего s_1=4 или t_1=2) мы за p-2=3 итерации дошли до s_{p-1}=t_{p-1}=0, а за ещё одну — до z^{2^{p-1}}=-1 (и, соответственно, s_p=-2, t_p=-1). Если бы мы ещё один раз возвели в квадрат — воткнулись бы в единицу, которая для возведения в квадрат неподвижная точка.

Давайте посмотрим, почему мы этого не сделали — и почему всё было настолько « близко », что стоит M оказаться составным, и ничего не получится. И заодно выясним много интересного, что живёт «по соседству».

Во-первых, если у нас есть поле из q элементов (где q нечётное), то есть очень симпатичная формула, позволяющая проверить, является ли его элемент квадратом (в этом поле!). А именно — малая теорема Ферма для простого q и теорема Лагранжа для q=M^k утверждают, что w^{q-1}=1 для всех элементов поля.

Возьмём какое-нибудь ненулевое w, и давайте рассмотрим половинную степень — w^{(q-1)/2}.
Утверждение. w^{(q-1)/2}=1, если w это квадрат, и w^{(q-1)/2}=-1, если не квадрат.