Дочка принесла из универа задачку с очень простой формулировкой, но уже пару недель не получается сдать (проверяет стандартная система автопроверки). Сейчас сдача уже не актуальна, на зачет баллов набрала, но вопрос, как же надо было решать, по прежнему свербит. Поэтому, если кто-то ткнет носом в правильный вариант, буду весьма признателен.
Задачка такая: дано простое число 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++) — возможно, это как раз и есть главная подсказка, но мы ее не поняли.
K>Задачка такая: дано простое число p<10^9, положительное целое число a<p и диапазон чисел [l;r], 0<=l<=r<p. Требуется вывести в порядке возрастания все числа из диапазона, которые представимы в виде a^k mod p для какого либо k. Гарантируется, что таких чисел (которые составляют ответ) не более 100 штук.
В приведённой формулировке k сверху вообще неограничен
K>>Задачка такая: дано простое число p<10^9, положительное целое число a<p и диапазон чисел [l;r], 0<=l<=r<p. Требуется вывести в порядке возрастания все числа из диапазона, которые представимы в виде a^k mod p для какого либо k. Гарантируется, что таких чисел (которые составляют ответ) не более 100 штук.
B>В приведённой формулировке k сверху вообще неограничен
В какой-то момент закольцуется, можно считать, что ограничен p.
K>>Задачка такая: дано простое число p<10^9, положительное целое число a<p и диапазон чисел [l;r], 0<=l<=r<p. Требуется вывести в порядке возрастания все числа из диапазона, которые представимы в виде a^k mod p для какого либо k. Гарантируется, что таких чисел (которые составляют ответ) не более 100 штук.
B>В приведённой формулировке k сверху вообще неограничен
Туплю — модуль-то может зациклиться. Попробуй тупо ограничить k на уровне алгоритма и погонять с разными значениями (10, 20, 30, ...) — может там ошибка в условии
Как идея. По формуле Эйлера x/a(mod p) = x*a^(p-2) (mod p). Можно перебирать x от l до r и пытаться делить в попытке прийти к a. Если получилось — печатаем. Если нет (зациклились), возможно запоминаем все значения как неудачные чтобы в будущем быстро отбрасываться.
Здравствуйте, vsb, Вы писали:
vsb>Как идея. По формуле Эйлера x/a(mod p) = x*a^(p-2) (mod p). Можно перебирать x от l до r и пытаться делить в попытке прийти к a. Если получилось — печатаем. Если нет (зациклились), возможно запоминаем все значения как неудачные чтобы в будущем быстро отбрасываться.
Мне кажется, это тот же проход по циклу, только в обратном направлении, нет? И при длинном цикле можно надолго застрять на конкретном x.
Здравствуйте, 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 / 815730721return 0;
}
Здравствуйте, 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й тест).
Здравствуйте, 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.
Здравствуйте, 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;
}
Здравствуйте, 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-м.
Так что, увы, несмотря на большой шаг вперед, и этого пока не хватает
Здравствуйте, 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.
Спасибо, я думал про это, но не увидел тут перспективы существенного ускорения, поэтому не попробовал. Ну и, к сожалению, для первокурсника, которому не читали теорию чисел, этот вариант не катит точно...
Здравствуйте, kfmn, Вы писали:
K>Спасибо, я думал про это, но не увидел тут перспективы существенного ускорения, поэтому не попробовал. Ну и, к сожалению, для первокурсника, которому не читали теорию чисел, этот вариант не катит точно...
Если слетает по времени на тестах, где a генерирует всю группу, то должно помочь. А так да, это из теории чисел, конечно.
Здравствуйте, 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. Далее корректируем порог для наивного алгоритма, и всё должно пройти с большим запасом по времени. Ну и предыдущую оптимизацию можно откатить — она теперь не нужна.
Так как 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;
}
Здравствуйте, foxafox, Вы писали:
F>Так как p — простое, то кольцо вычетов циклично, и любая подгруппа циклична. Можно попробовать такое решение:
Они все цикличны, но длинна цикла может быть меньше чем надо. Например A=p-1
Здесь все числа из диапазона могут быть представимы в виде 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 согласно теореме Ферма.
Здравствуйте, AndreyM16, Вы писали:
AM>Здесь все числа из диапазона могут быть представимы в виде a^k mod p, т.к. все числа вида a^k представляют циклическую группу порядка p.
В эту мультипликативную группу не входит ноль — так как у него нет обратного элемента. Порядок такой группы: p-1.
Не каждый элемент этой группы является порождающим. Например достаточно рассмотреть степени двойки в группе p=7: 2, 4, 1
Здравствуйте, kfmn, Вы писали: K>Спасибо за потраченное время, но, увы, этот вариант тоже не устроил систему (вероятно, там более слабое железо ).
Вам же не нужно само значение логарифма, достаточно понять, представляется ли. А представляется число тогда и только тогда, когда логарифм числа (по какому-нибудь образующему) кратен логарифму основания (по модулю (p-1)). Итого, сначала вычисляете наименьшую натуральную степень m, в которую надо возвести основание a чтобы получить 1. Затем для каждого числа x от l до r возводите его (быстрым возведением) в степень m. Если получилась единица — x представимо, иначе нет.
Второй, самый тупой вариант, как раз для первокурсников: неспеша возводите a в натуральные степени (одно умножение), пока не получите снова a. Каждый результат проверяете на попадание в отрезок. Выписываете список попавших. В качестве защиты от злыдней, если a — образующий (проверяется мгновенно), то сразу выписываете весь отрезок.
Update 1:
Реализовал первый вариант.
Последний тесткейс от Sergei MO выкинул, так как в нём p=13^8, что ни разу не есть простое число.
В каждом тесткейсе нашлось ровно по 100 чисел.
Слегка не укладывается в миллисекунду, попробуйте, если тестовая система ещё доступна.
Код
#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 / 815730721return 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 гигов оперативы. А вот умножение, даже с делением по модулю, уже гораздо дешевле.
Здравствуйте, kfmn, Вы писали:
K>Дочка принесла из универа задачку с очень простой формулировкой, но уже пару недель не получается сдать (проверяет стандартная система автопроверки). Сейчас сдача уже не актуальна, на зачет баллов набрала, но вопрос, как же надо было решать, по прежнему свербит. Поэтому, если кто-то ткнет носом в правильный вариант, буду весьма признателен.
Ну видимо, ко всем, кому удалось решить эту задачку, потом подходит специально обученный дяденька из ФСБ, и предлагает перейти к ним учиться, на криптоаналитика