Хорошо, а чем плох (или, наоборот, хорош), такой подход: вот, мы же знаем, что

, - это лишь фасад, мираж, а
на самом деле рулят

. Мы так же знаем, что

, где

принадлежит определенной последовательности и слишком большим быть не должно (тут конечно вопрос, что такое "слишком", и традиционное "когда остановиться?"). При этом, точность приближения тем лучше, чем больше величина

, или, приближенно, чем меньше

. Тогда, для простого

, можно попробовать перебирать значения

по возрастанию (причем, не все

будут разрешены по соображениям взаимности) и в лоб извлекать из них корни по модулю

, - здесь же есть быстрые алгоритмы (не имел реального опыта реализации, но вот пишут:
Алгоритм Тонелли—Шенкса; есть некоторые оговорки, но, на первый взгляд, алгоритм не выглядит устрашающим/безнадежным). Если же нам еще больше повезло и

, то это вообще сводится к задаче возведения в степень по модулю:

, где точно быстрые алгоритмы есть. Хочу попробовать, когда будет чуть больше досуга, в ближайшую среду, видимо