Поиск степеней по модулю простого числа
От: kfmn Россия  
Дата: 17.12.18 14:19
Оценка: 4 (1)
Всем доброго дня!

Дочка принесла из универа задачку с очень простой формулировкой, но уже пару недель не получается сдать (проверяет стандартная система автопроверки). Сейчас сдача уже не актуальна, на зачет баллов набрала, но вопрос, как же надо было решать, по прежнему свербит. Поэтому, если кто-то ткнет носом в правильный вариант, буду весьма признателен.

Задачка такая: дано простое число p<10^9, положительное целое число a<p и диапазон чисел [l;r], 0<=l<=r<p. Требуется вывести в порядке возрастания все числа из диапазона, которые представимы в виде a^k mod p для какого либо k. Гарантируется, что таких чисел (которые составляют ответ) не более 100 штук.
Ограничения по памяти 256Mb, ограничение по времени 3.5 сек.
Все варианты, которые пробовали, рано или поздно вылетают по TimeLimit'у.

Пробовали следующие варианты:
1. Наивный подход: последовательно вычисляем все степени числа a по модулю p, запихивая результаты в std::set или взводя элементы в vector<bool>, как только нашли элемент, который уже встречался, выходим, т.к. дальше пойдет повторение. Этот вариант не проходит тест раньше всех.
2. Модификация первого подхода: посчитать, какая степень а разделяет (r+1) и (p+l), по выходу из диапазона [l,r] сразу домножать на эту степень, а не постепенно, много раз на a. Проходит всего на 1 тест больше.
3. Перебор всех чисел от l до r и проверка того, может ли это число быть степенью a по модулю p. Проверяли с помощью алгоритма поиска (дискретного логарифма в кольце вычетов по простому модулю).

Комбинация 2 и 3 (в зависимости от величины (r-l)) проходит на 4 теста больше, чем наивный подход, но все равно валится по TimeLimit'у.

Варианты 1 и 2 не годятся, если a — первообразный корень по модулю достаточно большого p (просто цикл получается чересчур длинным), вариант 3 — слишком "тяжелый" для больших диапазонов [l;r].
Вероятно, надо как-то использовать информацию о том, что в ответе не более 100 чисел, но такое возможно в обоих перечисленных раскладах (длинный цикл, но узкий диапазон [l;r], либо наоборот, диапазон широкий, но цикл короткий).

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

В общем, буду благодарен за любые подсказки.

P.S. Задача для первокурсников, серьезных знаний по теории чисел не предполагается (даже тех, о которых я упомянул).
P.P.S. Задачка появилась в контесте на стандартные структуры данных (т.е. те, которые доступны в библиотеках типа STL для C++) — возможно, это как раз и есть главная подсказка, но мы ее не поняли.
Re: Поиск степеней по модулю простого числа
От: bzig  
Дата: 18.12.18 00:04
Оценка:
K>Задачка такая: дано простое число p<10^9, положительное целое число a<p и диапазон чисел [l;r], 0<=l<=r<p. Требуется вывести в порядке возрастания все числа из диапазона, которые представимы в виде a^k mod p для какого либо k. Гарантируется, что таких чисел (которые составляют ответ) не более 100 штук.

В приведённой формулировке k сверху вообще неограничен
Re[2]: Поиск степеней по модулю простого числа
От: vsb Казахстан  
Дата: 18.12.18 00:09
Оценка: +2
Здравствуйте, bzig, Вы писали:


K>>Задачка такая: дано простое число p<10^9, положительное целое число a<p и диапазон чисел [l;r], 0<=l<=r<p. Требуется вывести в порядке возрастания все числа из диапазона, которые представимы в виде a^k mod p для какого либо k. Гарантируется, что таких чисел (которые составляют ответ) не более 100 штук.


B>В приведённой формулировке k сверху вообще неограничен


В какой-то момент закольцуется, можно считать, что ограничен p.
Re[2]: Поиск степеней по модулю простого числа
От: bzig  
Дата: 18.12.18 00:13
Оценка:
Здравствуйте, bzig, Вы писали:


K>>Задачка такая: дано простое число p<10^9, положительное целое число a<p и диапазон чисел [l;r], 0<=l<=r<p. Требуется вывести в порядке возрастания все числа из диапазона, которые представимы в виде a^k mod p для какого либо k. Гарантируется, что таких чисел (которые составляют ответ) не более 100 штук.


B>В приведённой формулировке k сверху вообще неограничен


Туплю — модуль-то может зациклиться. Попробуй тупо ограничить k на уровне алгоритма и погонять с разными значениями (10, 20, 30, ...) — может там ошибка в условии
Re: Поиск степеней по модулю простого числа
От: vsb Казахстан  
Дата: 18.12.18 00:40
Оценка:
Как идея. По формуле Эйлера x/a(mod p) = x*a^(p-2) (mod p). Можно перебирать x от l до r и пытаться делить в попытке прийти к a. Если получилось — печатаем. Если нет (зациклились), возможно запоминаем все значения как неудачные чтобы в будущем быстро отбрасываться.
Отредактировано 18.12.2018 0:42 vsb . Предыдущая версия . Еще …
Отредактировано 18.12.2018 0:42 vsb . Предыдущая версия .
Re[2]: Поиск степеней по модулю простого числа
От: kfmn Россия  
Дата: 18.12.18 05:39
Оценка:
Здравствуйте, vsb, Вы писали:

vsb>Как идея. По формуле Эйлера x/a(mod p) = x*a^(p-2) (mod p). Можно перебирать x от l до r и пытаться делить в попытке прийти к a. Если получилось — печатаем. Если нет (зациклились), возможно запоминаем все значения как неудачные чтобы в будущем быстро отбрасываться.


Мне кажется, это тот же проход по циклу, только в обратном направлении, нет? И при длинном цикле можно надолго застрять на конкретном x.
Re: Поиск степеней по модулю простого числа
От: Sergei MO Россия  
Дата: 23.12.18 14:12
Оценка: 4 (1)
Здравствуйте, kfmn, Вы писали:

K>Вероятно, надо как-то использовать информацию о том, что в ответе не более 100 чисел

Если предположить, что элементы циклической группы с образующим элементом a более-менее равномерно распределены в диапазоне [0, p), то по величине r-l можно сделать оценку размера этой группы.

Я по вашему описанию набросал решение скомбинированное из 1 и 3 вариантов — у меня оно укладывается в 2.5 секунды и 130 Мб памяти. Скорее всего, ваш подход верный, но реализация неэффективная.

  Результаты замеров времени для двух вариантов:
R-L: 1248  naive: 4524 ms  discrete log: 1288 ms  ok
R-L: 1437  naive: 3992 ms  discrete log: 1492 ms  ok
R-L: 1593  naive: 3544 ms  discrete log: 1662 ms  ok
R-L: 2126  naive: 2763 ms  discrete log: 2256 ms  ok
R-L: 2459  naive: 2363 ms  discrete log: 2587 ms  ok
R-L: 2834  naive: 1985 ms  discrete log: 2990 ms  ok
R-L: 2853  naive: 2006 ms  discrete log: 3018 ms  ok
R-L: 3822  naive: 1498 ms  discrete log: 4025 ms  ok
R-L: 815730721  naive: 59 ms

  Само решение:
#include <iostream>
#include <vector>
#include <algorithm>
#define NOMINMAX
#include <windows.h>
using namespace std;

vector<int> naive(int A, int P, int L, int R)
{
    vector<int> res;
    vector<bool> bitmap(P, false);

    int x = 1;
    while (true)
    {
        if (bitmap[x])
        {
            break;
        }
        bitmap[x] = true;
        if (x >= L && x <= R)
        {
            res.push_back(x);
        }
        x = (long long)x * A % P;
    }

    sort(res.begin(), res.end());
    return res;
}

vector<int> discretelog(int A, int P, int L, int R)
{
    int h = sqrt((double)P) + 1;
    vector<int> values(h + 1);

    int c = 1;
    values[0] = 1;
    for (int i = 1; i <= h; i++)
    {
        c = (long long)c * A % P;
        values[i] = (long long)values[i - 1] * A % P;
    }

    vector<bool> bitmap(P, false);
    int t = c;
    for (int i = 1; i <= h; i++)
    {
        bitmap[t] = true;
        t = (long long)t * c % P;
    }

    vector<int> res;
    for (int x = L; x <= R; x++)
    {
        for (int i = 0; i <= h; i++)
        {
            if (bitmap[(long long)values[i] * x % P])
            {
                res.push_back(x);
                break;
            }
        }
    }

    return res;
}

vector<int> solve(int A, int P, int L, int R)
{
    if (R - L < 2400)
    {
        return discretelog(A, P, L, R);
    }
    else
    {
        return naive(A, P, L, R);
    }
}

void test(int A, int P, int L, int R)
{
    LARGE_INTEGER start, end, freq;
    QueryPerformanceFrequency(&freq);
    cout << "R-L: " << R - L;

    QueryPerformanceCounter(&start);
    vector<int> resn = naive(A, P, L, R);
    QueryPerformanceCounter(&end);
    cout << "  naive: " << (end.QuadPart - start.QuadPart) * 1000 / freq.QuadPart << " ms";

    if (R - L > 10000)
    {
        cout << endl;
        return;
    }

    QueryPerformanceCounter(&start);
    vector<int> resd = discretelog(A, P, L, R);
    QueryPerformanceCounter(&end);
    cout << "  discrete log: " << (end.QuadPart - start.QuadPart) * 1000 / freq.QuadPart << " ms";

    bool correct = (resn.size() == resd.size());
    if (correct)
    {
        for (int i = 0; i < resn.size(); i++)
        {
            if (resn[i] != resd[i])
            {
                correct = false;
                break;
            }
        }
    }
    cout << (correct ? "  ok" : "  error") << endl;
}

int main()
{
    test(29, 999999937, 204411352, 204412600); // 124999993 / 1248
    test(7, 999999937, 450141522, 450142959);  // 111111105 / 1437
    test(32, 999999751, 239615046, 239616639); //  99999976 / 1593
    test(41, 999999937, 395909863, 395911989); //  76923073 / 2126
    test(6, 999999751, 217939683, 217942142);  //  66666651 / 2459
    test(6, 999999883, 724460982, 724463816);  //  55555550 / 2834
    test(35, 999999937, 399778083, 399780936); //  55555553 / 2853
    test(2, 999999937, 432617324, 432621146);  //  41666665 / 3822

    test(13, 815730721, 0, 815730721);  //  10 / 815730721
    return 0;
}
Re[2]: Поиск степеней по модулю простого числа
От: kfmn Россия  
Дата: 24.12.18 09:21
Оценка:
Здравствуйте, Sergei MO, Вы писали:

SM>Здравствуйте, kfmn, Вы писали:


K>>Вероятно, надо как-то использовать информацию о том, что в ответе не более 100 чисел

SM>Если предположить, что элементы циклической группы с образующим элементом a более-менее равномерно распределены в диапазоне [0, p), то по величине r-l можно сделать оценку размера этой группы.

SM>Я по вашему описанию набросал решение скомбинированное из 1 и 3 вариантов — у меня оно укладывается в 2.5 секунды и 130 Мб памяти. Скорее всего, ваш подход верный, но реализация неэффективная.


Спасибо за потраченное время, но, увы, этот вариант тоже не устроил систему (вероятно, там более слабое железо ). Признаюсь, я не сильно заморачивался на подбор порога для разности R-L, пробовал 1000, 3000 и 10000. На 1000 удавалось пройти до 39 теста, на 3000 и выше — только до 36-го, так же, как и при использовании дискретного логарифма во всех случаях. И порог 2400 тоже оказался из этой категории (36й тест).
Re: Поиск степеней по модулю простого числа
От: andyp  
Дата: 24.12.18 12:14
Оценка: 4 (1)
Здравствуйте, kfmn, Вы писали:

K>...В общем, буду благодарен за любые подсказки.



Есть еще одна идея. Перед началом вычислений

Я бы разложил p-1 на простые множители:

p-1 = p1^a1 * p2^a2 * pk^ak;

и проверил бы

a^((p-1)/pj) mod p

на равенство 1.

Если ни одно из равенств не выполняется, то a — генерирующий элемент группы (Z/pZ)* и в ответе будут первые 100 чисел интервала (или весь интервал, что короче)

PS Это позволит иметь дело с подгруппой размера по крайней мере в 2 раза меньше, чем p-1.
Отредактировано 24.12.2018 12:17 andyp . Предыдущая версия .
Re[3]: Поиск степеней по модулю простого числа
От: Sergei MO Россия  
Дата: 24.12.18 13:09
Оценка: 4 (1)
Здравствуйте, kfmn, Вы писали:

K>Признаюсь, я не сильно заморачивался на подбор порога для разности R-L, пробовал 1000, 3000 и 10000. На 1000 удавалось пройти до 39 теста, на 3000 и выше — только до 36-го, так же, как и при использовании дискретного логарифма во всех случаях. И порог 2400 тоже оказался из этой категории (36й тест).


Получается, что в 36-м тесте 1000 < R-L < 2400, и время работы дискретного логарифма не укладывается в лимит.

Дискретный логарифм можно ускорить раза в два, если использовать тот факт, что проверяемые числа идут последовательно. Тогда можно поменять местами вложенный цикл с внешним и избавиться от длинного умножения и взятия модуля во внутреннем цикле.
  Вот так:
vector<int> discretelog(int A, int P, int L, int R)
{
    int h = sqrt((double)P) + 1;
    vector<int> values(h + 1);

    int c = 1;
    values[0] = 1;
    for (int i = 1; i <= h; i++)
    {
        c = (long long)c * A % P;
        values[i] = (long long)values[i - 1] * A % P;
    }

    vector<bool> bitmap(P, false);
    int t = c;
    for (int i = 1; i <= h; i++)
    {
        bitmap[t] = true;
        t = (long long)t * c % P;
    }

    set<int> tmp;
    for (int i = 0; i <= h; i++)
    {
        int value = (long long)values[i] * L % P;
        for (int x = 0; x <= R - L; x++)
        {
            if (bitmap[value])
            {
                tmp.insert(x + L);
            }
            value += values[i];
            if (value >= P)
            {
                value -= P;
            }
        }
    }

    vector<int> res;
    for (auto it = tmp.begin(); it != tmp.end(); ++it)
    {
        res.push_back(*it);
    }
    return res;
}
Re[4]: Поиск степеней по модулю простого числа
От: kfmn Россия  
Дата: 25.12.18 13:39
Оценка:
Здравствуйте, Sergei MO, Вы писали:

SM>Здравствуйте, kfmn, Вы писали:


K>>Признаюсь, я не сильно заморачивался на подбор порога для разности R-L, пробовал 1000, 3000 и 10000. На 1000 удавалось пройти до 39 теста, на 3000 и выше — только до 36-го, так же, как и при использовании дискретного логарифма во всех случаях. И порог 2400 тоже оказался из этой категории (36й тест).


SM>Получается, что в 36-м тесте 1000 < R-L < 2400, и время работы дискретного логарифма не укладывается в лимит.


SM>Дискретный логарифм можно ускорить раза в два, если использовать тот факт, что проверяемые числа идут последовательно. Тогда можно поменять местами вложенный цикл с внешним и избавиться от длинного умножения и взятия модуля во внутреннем цикле.


С таким вариантом и порогом на r-l в 2400 TimeLimit случился уже на 61-м тесте! Если увеличить порог до 3500 — то на 107-м!, а если до 4000, то на 37-м. Т.е. оптимум где-то между 3500 и 4000...
Стал перебирать, выяснилось, что еще на 3800 валится на 107-м тесте, а на 3900 — уже на 37-м.

Так что, увы, несмотря на большой шаг вперед, и этого пока не хватает
Re[2]: Поиск степеней по модулю простого числа
От: kfmn Россия  
Дата: 25.12.18 13:41
Оценка:
Здравствуйте, andyp, Вы писали:

A>Здравствуйте, kfmn, Вы писали:


K>>...В общем, буду благодарен за любые подсказки.



A>Есть еще одна идея. Перед началом вычислений


A>Я бы разложил p-1 на простые множители:


A>p-1 = p1^a1 * p2^a2 * pk^ak;


A>и проверил бы


A>a^((p-1)/pj) mod p


A>на равенство 1.


A>Если ни одно из равенств не выполняется, то a — генерирующий элемент группы (Z/pZ)* и в ответе будут первые 100 чисел интервала (или весь интервал, что короче)


A>PS Это позволит иметь дело с подгруппой размера по крайней мере в 2 раза меньше, чем p-1.



Спасибо, я думал про это, но не увидел тут перспективы существенного ускорения, поэтому не попробовал. Ну и, к сожалению, для первокурсника, которому не читали теорию чисел, этот вариант не катит точно...
Re[3]: Поиск степеней по модулю простого числа
От: andyp  
Дата: 25.12.18 14:33
Оценка:
Здравствуйте, kfmn, Вы писали:

K>Спасибо, я думал про это, но не увидел тут перспективы существенного ускорения, поэтому не попробовал. Ну и, к сожалению, для первокурсника, которому не читали теорию чисел, этот вариант не катит точно...


Если слетает по времени на тестах, где a генерирует всю группу, то должно помочь. А так да, это из теории чисел, конечно.
Re[5]: Поиск степеней по модулю простого числа
От: Sergei MO Россия  
Дата: 25.12.18 15:11
Оценка:
Здравствуйте, kfmn, Вы писали:

K>С таким вариантом и порогом на r-l в 2400 TimeLimit случился уже на 61-м тесте! Если увеличить порог до 3500 — то на 107-м!, а если до 4000, то на 37-м. Т.е. оптимум где-то между 3500 и 4000...

K>Стал перебирать, выяснилось, что еще на 3800 валится на 107-м тесте, а на 3900 — уже на 37-м.

Фигасе у них тестов. Это обычная задача или повышенной сложности для самых умных?

Тем не менее, я понял, где можно сильно ускорить решение.
В дискретном логарифме есть три основных цикла со следующими количествами итераций:
1) sqrt(P)
2) sqrt(P)
3) sqrt(P) * (R-L)

То есть третий цикл сейчас самый затратный. Это потому что я, не подумав, взял H = sqrt(P), как в википедии.
А вот если взять, например, H = 100, то получим количества итераций:
1) 100
2) P/100
3) 100 * (R-L)

Тогда затраты на третий цикл превысят второй цикл только при R-L > 100000. Далее корректируем порог для наивного алгоритма, и всё должно пройти с большим запасом по времени. Ну и предыдущую оптимизацию можно откатить — она теперь не нужна.
Re: Поиск степеней по модулю простого числа
От: foxafox  
Дата: 26.12.18 22:37
Оценка:
Так как p — простое, то кольцо вычетов циклично, и любая подгруппа циклична. Можно попробовать такое решение:


std::set<unsigned> solve_set(unsigned A, unsigned P, unsigned L, unsigned R)
{
    std::set<unsigned> res;

    unsigned x = 1;

    do
    {
        x = (std::uint64_t) x * A % P;

        if (x >= L && x <= R)
        {
          if (!res.insert(x) .second)
            break;
        }
    }
    while (x != 1);

    return res;
}
Re[2]: Поиск степеней по модулю простого числа
От: kov_serg Россия  
Дата: 27.12.18 06:28
Оценка:
Здравствуйте, foxafox, Вы писали:

F>Так как p — простое, то кольцо вычетов циклично, и любая подгруппа циклична. Можно попробовать такое решение:

Они все цикличны, но длинна цикла может быть меньше чем надо. Например A=p-1
Re: Поиск степеней по модулю простого числа
От: AndreyM16  
Дата: 26.05.19 18:35
Оценка:
Здравствуйте, kfmn, Вы писали:

Здесь все числа из диапазона могут быть представимы в виде a^k mod p, т.к. все числа вида a^k представляют циклическую группу порядка p. Это можно сразу получить из следствия теоремы Лагранжа. (Группа порядка p, где p — простое число, циклична.)
Или можно вывести из теоремы Ферма: допустим a^k1 = a^k2 mod p, (k1 > k2), тогда a^k2 * (a^(k1-k2) — 1) mod p или a^(k1-k2) = 1 mod p. А минимальное число удовлетворяющее данному сравнению это p-1 согласно теореме Ферма.
Re[2]: Поиск степеней по модулю простого числа
От: Evgeny.Panasyuk Россия  
Дата: 27.05.19 11:48
Оценка:
Здравствуйте, AndreyM16, Вы писали:

AM>Здесь все числа из диапазона могут быть представимы в виде a^k mod p, т.к. все числа вида a^k представляют циклическую группу порядка p.


В эту мультипликативную группу не входит ноль — так как у него нет обратного элемента. Порядок такой группы: p-1.
Не каждый элемент этой группы является порождающим. Например достаточно рассмотреть степени двойки в группе p=7: 2, 4, 1
Re[3]: Поиск степеней по модулю простого числа
От: cures Россия cures.narod.ru
Дата: 27.05.19 16:22
Оценка:
Здравствуйте, kfmn, Вы писали:

K>Спасибо за потраченное время, но, увы, этот вариант тоже не устроил систему (вероятно, там более слабое железо ).


Вам же не нужно само значение логарифма, достаточно понять, представляется ли. А представляется число тогда и только тогда, когда логарифм числа (по какому-нибудь образующему) кратен логарифму основания (по модулю (p-1)). Итого, сначала вычисляете наименьшую натуральную степень m, в которую надо возвести основание a чтобы получить 1. Затем для каждого числа x от l до r возводите его (быстрым возведением) в степень m. Если получилась единица — x представимо, иначе нет.

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

Update 1:

Реализовал первый вариант.
Последний тесткейс от Sergei MO выкинул, так как в нём p=13^8, что ни разу не есть простое число.
В каждом тесткейсе нашлось ровно по 100 чисел.

  Для остальных тайминги:
a=29, p=999999937, l=204411352, r=204412600
0.000429843s
a=7, p=999999937, l=450141522, r=450142959
0.000491738s
a=32, p=999999751, l=239615046, r=239616639
0.000556290s
a=41, p=999999937, l=395909863, r=395911989
0.000637724s
a=6, p=999999751, l=217939683, r=217942142
0.000783899s
a=6, p=999999883, l=724460982, r=724463816
0.000936862s
a=35, p=999999937, l=399778083, r=399780936
0.000885491s
a=2, p=999999937, l=432617324, r=432621146
0.001143331s


Слегка не укладывается в миллисекунду, попробуйте, если тестовая система ещё доступна.

  Код
#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <time.h>
#include <vector>
#include <algorithm>

typedef std::vector<uint32_t> uint32_v;


static double timediff(const timespec *tp0, const timespec *tp1)
{
    return double (tp1->tv_sec - tp0->tv_sec) +
        1e-9 * (tp1->tv_nsec - tp0->tv_nsec);
}

static inline uint32_t mul(uint32_t a, uint32_t b, uint32_t p)
{
    return uint64_t(a) * b % p;
}

static inline uint32_t pow(uint32_t a, uint32_t n, uint32_t p)
{
    uint32_t res = 1, b = a;
    while (n) {
        if (n & 1)
            res = mul(res, b, p);
        b = mul(b, b, p);
        n >>= 1;
    }
    return res;
}

struct Factor
{
    uint32_t l;
    uint32_t ps[9];
    uint32_t ks[9];
};


void factor(uint32_t a, Factor &f)
{
    uint32_t l = 0;

    assert (a);
    if (a % 2 == 0) {
        uint32_t k = 0;
        while (a % 2 == 0) {
            a /= 2;
            k++;
        }
        f.ps[l] = 2;
        f.ks[l] = k;
        l++;
    }

    for (uint32_t p = 3; ; p += 2) {
        uint32_t q = a / p, r = a % p;
        if (q < p)
            break;
        if (r == 0) {
            uint32_t k = 0;
            while (a % p == 0) {
                a /= p;
                k++;
            }
            f.ps[l] = p;
            f.ks[l] = k;
            l++;
        }
    }

    if (a > 1) {
        f.ps[l] = a;
        f.ks[l] = 1;
        l++;
    }

    f.l = l;
}

uint32_t find_deg(uint32_t a, uint32_t p, Factor fq)
{
    uint32_t q = p - 1;
    assert (pow(a, q, p) == 1);
    for (uint32_t l = 0; l < fq.l; l++) {
        uint32_t pl = fq.ps[l], kl = fq.ks[l];
        while (kl) {
            uint32_t qq = q / pl;
            if (pow(a, q / pl, p) != 1)
                break;
            q /= pl;
            kl--;
        }
    }
    //printf("a=%u, p=%u, deg=%u\n", a, p, q);
    assert (q);
    return q;
}

void solve(uint32_t a, uint32_t p, uint32_t l, uint32_t r, uint32_v &v)
{
    assert (a < p);
    assert (l <= r);
    uint32_t q = p - 1;
    Factor fp, fq;
    factor(p, fp);
    assert (fp.l == 1 && fp.ks[0] == 1);
    factor(q, fq);
    
    v.clear();
    if (a <= 1) {
        if (l <= a && a <= r)
            v.push_back(a);
        return;
    }
    
    uint32_t m = find_deg(a, p, fq);
    for (uint32_t x = l; x <= r; x++) {
        if (pow(x, m, p) == 1)
            v.push_back(x);
    }
    printf("size=%u\n", (uint32_t) v.size());
    std::sort(v.begin(), v.end());
}

void test(uint32_t a, uint32_t p, uint32_t l, uint32_t r)
{
    uint32_v v;
    v.reserve(100);
    
    printf("a=%u, p=%u, l=%u, r=%u\n", a, p, l, r);
    timespec tp0, tp1;
    clock_gettime(CLOCK_MONOTONIC, &tp0);
    solve(a, p, l, r, v);
    clock_gettime(CLOCK_MONOTONIC, &tp1);
    printf("%.9fs\n", timediff(&tp0, &tp1));    
}


int main()
{
    test(29, 999999937, 204411352, 204412600); // 124999993 / 1248
    test(7, 999999937, 450141522, 450142959);  // 111111105 / 1437
    test(32, 999999751, 239615046, 239616639); //  99999976 / 1593
    test(41, 999999937, 395909863, 395911989); //  76923073 / 2126
    test(6, 999999751, 217939683, 217942142);  //  66666651 / 2459
    test(6, 999999883, 724460982, 724463816);  //  55555550 / 2834
    test(35, 999999937, 399778083, 399780936); //  55555553 / 2853
    test(2, 999999937, 432617324, 432621146);  //  41666665 / 3822
    //test(13, 815730721, 0, 815730721);  //  10 / 815730721
    return 0;
}


Update 2:

Понял, что жёстко туплю. 3.5 секунды — это явно на второй вариант, для которого ничего не надо знать. Поправил код Sergei MO, убрав жуткие битмапы и, на всякий случай, сведя проверку на попадание в интервал к одному сравнению.

  Код
void naive(uint32_t a, uint32_t p, uint32_t l, uint32_t r, uint32_v &v)
{
    assert (0 < a && a < p);
    v.clear();

    uint32_t x = 1;
    for (;;) {
        if ((x - l) <= (r - l) /*l <= x && x <= r*/)
            v.push_back(x);
        x = mul(x, a, p);
        if (x == 1)
            break;
    }
    std::sort(v.begin(), v.end());
}


На ideone получил на тестах Sergei MO (опять же, не считая последнего — неправильного) максимальный тайминг 1.113 секунд (первый тест). Здесь mul — из первого куска кода, naive можно вызывать в test вместо solve.
Заводить вектор на миллиард булов — слишком жёстко, тут никто не успеет. Повторение случится именно с единицы, если a не равно нулю по модулю p, что верно по условию, но можно и отдельно разобрать случай a=0, он тривиален.
По той же причине не срослось и с сетом, он на самом деле есть уравновешенное бинарное дерево — очень дорогое удовольствие, сотни тактов на добавление или проверку, каковых может быть порядка миллиарда, да ещё и байт по 20 на каждый добавленный элемент — совсем печалька, 20 гигов оперативы. А вот умножение, даже с делением по модулю, уже гораздо дешевле.
Отредактировано 27.05.2019 18:39 cures . Предыдущая версия . Еще …
Отредактировано 27.05.2019 17:59 cures . Предыдущая версия .
Re: Поиск степеней по модулю простого числа
От: Pzz Россия https://github.com/alexpevzner
Дата: 28.05.19 08:59
Оценка:
Здравствуйте, kfmn, Вы писали:

K>Дочка принесла из универа задачку с очень простой формулировкой, но уже пару недель не получается сдать (проверяет стандартная система автопроверки). Сейчас сдача уже не актуальна, на зачет баллов набрала, но вопрос, как же надо было решать, по прежнему свербит. Поэтому, если кто-то ткнет носом в правильный вариант, буду весьма признателен.


Ну видимо, ко всем, кому удалось решить эту задачку, потом подходит специально обученный дяденька из ФСБ, и предлагает перейти к ним учиться, на криптоаналитика
 
Подождите ...
Wait...
Пока на собственное сообщение не было ответов, его можно удалить.