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

[Without Name]

Хочу написать алгоритм Тонелли-Шенкса с ипользованием библиотеки GMP для использования в эллиптической криптографии. Описание алгоритма взял здесь https://ru.wikipedia.org/wiki/%D0%90%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC_%D0%A2%D0%BE%D0%BD%D0%B5%D0%BB%D0%BB%D0%B8_%E2%80%94_%D0%A8%D0%B5%D0%BD%D0%BA%D1%81%D0%B0:

Возникли странные проблемы с этим алгоритмом. Вот например, я хочу вычислить квадратный корень из $13 \mod 73$. Запускаю алгоритм пошагово:

1) $x^2 = 13 \mod 73$
$p = 73$
$p - 1 = 72 = 2^3 \cdot 9$
$s = 3, q = 9$

2) Нашел случайный невычет по указанному модулю:
$z = 17$
$c = z^q = 17^{9} \mod 73 = 63$


3) $r = {13}^{\frac{9+1}{2}} = 13^5 \mod 73 = 15$
$t = n^q = 13^9 \mod 73 = 51$
$m = 3$

Запускаю цикл

Находим $i: 0 < i < 3$

$i = 1, {t}^{2^i} = 51^2 \mod 73 = 46$
$i = 2, {t}^{2^i} = 51^4 \mod 73 = 72$

Перебрали все $i$ и не нашли такой $i$, что ${t}^{2^i} = 1 \mod 73$

Почему так?

Вот код моего алгоритма:

#include <iostream>
#include <ctime>
#include <chrono>
#include <gmp.h>

gmp_randstate_t state;

void ff_rand_init() {
    gmp_randinit_mt(state);
    gmp_randseed_ui(state, time(NULL));
}

void ff_rand(mpz_t n, mpz_t modulus) {
    mpz_urandomb(n, state, 256);
    mpz_mod(n, n, modulus);
}

bool ff_cmp_ui(mpz_t a, long b, mpz_t p) {
    mpz_t t;

    mpz_init_set(t, a);
    mpz_mod(t, t, p);

    int r = mpz_cmp_ui(t, b);

    mpz_clear(t);

    return r == 0;
}

bool ff_sqrt(mpz_t r, mpz_t n, mpz_t p) {
    mpz_t q;
    mpz_t t;
    mpz_t u;
    mpz_t v;
    mpz_t z;
    mpz_t c;
    mpz_t x;
    mpz_t y;
    mpz_t b;

    mpz_init(q);
    mpz_init(t);
    mpz_init(u);
    mpz_init(v);
    mpz_init(z);
    mpz_init(c);
    mpz_init(x);
    mpz_init(y);
    mpz_init(b);

    mpz_sub_ui(q, p ,1);

    long s = 0;
    while (mpz_tstbit(q, 0) == 0) {
        mpz_divexact_ui(q, q, 2);
        s++;
    }

    if (s == 1) {
       mpz_set(t, p);
       mpz_add_ui(t, t, 1);
       mpz_divexact_ui(t, t, 4);
       mpz_powm(r, n, t, p);

       mpz_clear(q);
       mpz_clear(t);

       return true;
    }

    long i = 1;
    while (i >= 0) {
        ff_rand(z, p);
        i = mpz_legendre(z, p);
    }

    mpz_powm(c, z, q, p);

    mpz_add_ui(u, q, 1);
    mpz_divexact_ui(u, u, 2);
    mpz_powm(r, n, u, p);
    mpz_powm(t, n, q, p);

    long m = s;

    while (true) {
        if (ff_cmp_ui(t, 1, p)) {
            mpz_clear(q);
            mpz_clear(t);
            mpz_clear(u);
            mpz_clear(z);
            mpz_clear(c);

            return true;
        }
        else {
            long i = 1;
            mpz_set(x, t);
            while (i < m) {
                mpz_powm_ui(x, x, 2, p);
                if (mpz_cmp_ui(x, 1)) {
                    break;
                }
                i++;
            }

            long e = m - i - 1;
            mpz_set(b, c);
            while (e > 0) {
                mpz_powm_ui(b, b, 2, p);
                e--;
            }
            mpz_mul(c, b, b);
            mpz_mod(c, c, p);

            mpz_mul(t, t, c);
            mpz_mod(t, t, p);

            mpz_mul(r, r, b);
            mpz_mod(r, r, p);

            m = i;
        }
    }

    return false;
}

int main()
{
    mpz_t r, m, x, y;

    mpz_init(r);
    mpz_init_set_str(m, "73", 10);
    mpz_init(x);
    mpz_init_set_str(y, "13", 10);

    ff_rand_init();

    if (ff_sqrt(x, y, m)) {
        gmp_printf("%Zd\n", x);
    }

    return 0;
}

cmake для сборки проекта:

cmake_minimum_required(VERSION 3.5)

project(finite-fields LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

set(GMP_LIB_PATH "/usr/local/lib")

add_executable(ff main.cpp)

target_link_directories(ff PRIVATE ${GMP_LIB_PATH})
target_link_libraries(ff gmp)

Алгоритм у меня не работает - уходит в бесконечный цикл, когда $p$ не сравнимо с 3 $\mod 4$

[dgwuqtj]

Зачем вообще писать с GMP, не имея рабочую версию с int? Числа маленькие, читать код было бы заметно проще...

[Without Name]

Вроде отловил место, где функция зависала. Вот исправленный код:

bool ff_sqrt(mpz_t r, mpz_t n, mpz_t p) {
    mpz_t q;
    mpz_t t;
    mpz_t u;
    mpz_t v;
    mpz_t z;
    mpz_t c;
    mpz_t x;
    mpz_t y;
    mpz_t b;

    mpz_init(q);
    mpz_init(t);
    mpz_init(u);
    mpz_init(v);
    mpz_init(z);
    mpz_init(c);
    mpz_init(x);
    mpz_init(y);
    mpz_init(b);

    mpz_sub_ui(q, p ,1);

    long s = 0;
    while (mpz_tstbit(q, 0) == 0) {
        mpz_divexact_ui(q, q, 2);
        s++;
    }

    if (s == 1) {
       mpz_set(t, p);
       mpz_add_ui(t, t, 1);
       mpz_divexact_ui(t, t, 4);
       mpz_powm(r, n, t, p);

       mpz_clear(q);
       mpz_clear(t);

       return true;
    }

    /*long i = 1;
    while (i >= 0) {
        ff_rand(z, p);
        i = mpz_legendre(z, p);
    }*/

    mpz_init_set_str(z, "17", 10);

    gmp_printf("Random non-residue: %Zd\n", z);

    mpz_powm(c, z, q, p);

    gmp_printf("c = %Zd\n", c);

    mpz_add_ui(u, q, 1);
    mpz_divexact_ui(u, u, 2);
    mpz_powm(r, n, u, p);
    mpz_powm(t, n, q, p);

    gmp_printf("r = %Zd\n", r);
    gmp_printf("t = %Zd\n", t);

    long m = s;

    while (true) {
        if (ff_cmp_ui(t, 1, p)) {
            mpz_clear(q);
            mpz_clear(t);
            mpz_clear(u);
            mpz_clear(z);
            mpz_clear(c);

            return true;
        }
        else {
            long i = 1;
            mpz_set(x, t);
            while (i < m) {
                mpz_powm_ui(x, x, 2, p);
                if (mpz_cmp_ui(x, 1) == 0) {
                    break;
                }
                i++;
            }

            gmp_printf("i = %d\n", i);

            long e = m - i - 1;
            mpz_set(b, c);
            while (e > 0) {
                mpz_powm_ui(b, b, 2, p);
                e--;
            }
            mpz_mul(c, b, b);
            mpz_mod(c, c, p);

            mpz_mul(t, t, c);
            mpz_mod(t, t, p);

            mpz_mul(r, r, b);
            mpz_mod(r, r, p);

            m = i;
        }
    }

    return false;
}

 if (mpz_cmp_ui(x, 1) == 0) {
                    break;
                }

-- 15.03.2025, 10:12 --

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

[RIP]

Without Name в сообщении #1678642 писал(а):

Вот например, я хочу вычислить квадратный корень из $13 \mod 73$.

$13$ — квадратичный невычет по модулю $73$.

[Without Name]

Да, уже заметил, спасибо. Моя функция правильно вычислила корень из $19 \mod 73$