Форум программистов, компьютерный форум, киберфорум
С++ для начинающих
Войти
Регистрация
Восстановить пароль
 
0 / 0 / 0
Регистрация: 15.03.2012
Сообщений: 5
1

Оптимизация поиска линейных конгруэнтных генераторов

30.11.2016, 05:03. Показов 487. Ответов 1

Здравствуйте! Хотел попросить помощи со следующим заданием: необходимо вернуть количество возможных линейных конгруэнтных генераторов (https://www.cyberforum.ru/cgi-bin/latex.cgi?{X}_{i+1} = (a{X}_{i} + c) % m), удовлетворяющих условиям Халла и Добелла (C и M должны быть взаимно простыми, A-1 делится без остатка всеми простыми делителями М, если М делится на 4, то и А-1 должно делиться на 4). Дополнительно задано условие: М должно иметь хотя бы D делителей. Всё работает, но слишком медленно. Буду благодарен за любой совет по оптимизации! Спасибо!

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
#include <map>
#include <set>
#include <vector>
#include <iostream>
#include <regex>
 
using namespace std;
 
pair<unsigned, unsigned> a_range;
unsigned a_range_special_first;
pair<unsigned, unsigned> c_range;
pair<unsigned, unsigned> m_range;
unsigned d;
 
vector<char> coprimes_cache;
 
void get_input(istream &input_str) {
 
    input_str >> a_range.first >> a_range.second;
    input_str >> c_range.first >> c_range.second;
    input_str >> m_range.first >> m_range.second;
    input_str >> d;
 
    for (unsigned a = a_range.first; a < a_range.second; a++) {
        if (a % 4 == 1) {
            a_range_special_first = a;
            break;
        }
    }
 
}
 
unsigned pow_modular(unsigned a, unsigned b, unsigned m) {
 
    unsigned result = 1;
    while (b > 0) {
        if (b & 1) result = (result * a) % m;
        b = b >> 1;
        a = (a * a) % m;
    }
    return result;
 
}
 
bool miller_rabin_test(unsigned d, unsigned n) {
 
    unsigned a = rand() % (n - 4) + 2;
    unsigned x = pow_modular(a, d, n);
    if (x == 1 || x == n - 1) return true;
    while (d != n - 1)  {
        x = (x * x) % n;
        d *= 2;
        if (x == 1) return false;
        if (x == n-1) return true;
    }
    return false;
 
}
 
bool is_prime(unsigned n, unsigned k) {
 
    if (n == 4 || n % 2 == 0) return false;
    if (n <= 3 || n == 5 || n == 7 || n == 11 || n == 13 || n == 17 || n == 19 || n == 23) return true;
    unsigned d = n - 1;
    while (d % 2 == 0) d /= 2;
    for (unsigned i = 0; i < k; i++) {
        if (!miller_rabin_test(d, n)) return false;
    }
    return true;
 
}
 
bool has_m_divisors(unsigned n, const unsigned &m) {
 
    if ((m == 2 && n > 1) ||
        (m == 3 && n > 2 && n % 2 == 0)) {
        return true;
    }
    unsigned n_div = 2;
    unsigned sqrt_res = (unsigned) ceil(sqrt((float) n));
    bool sqrt_exception = sqrt_res * sqrt_res == n;
    for (unsigned i = 2; i <= sqrt_res; i++) {
        if (n % i == 0 && (n_div += (sqrt_exception && i == sqrt_res ? 1 : 2)) >= m) return true;
    }
    return false;
 
}
 
unsigned pollard_rho(const unsigned &n, unsigned level = 1) {
 
    if (level >= 20 || n == 1) return n;
    if (n % 2 == 0) return 2;
    unsigned x, y, c = rand() % (n - 1) + 1;
    unsigned d = 1;
    unsigned n_try = 0;
    x = y = rand() % (n - 2) + 2;
    unsigned sN = ceil(sqrt(n));
    if (is_prime(n, 10)) return n;
    while (d == 1) {
        x = (unsigned) (pow((double) x, 2.0) + c) % n;
        y = (unsigned) (pow(pow((double) y, 2.0), 2.0) + c) % n;
        unsigned diff = (unsigned) abs((int) x - (int) y);
        d = __gcd(diff, n);
        if (++n_try >= sN) {
            return pollard_rho(n, level + 1);
        }
        if (d == n) {
            return pollard_rho(n, level + 1);
        }
    }
    if (!is_prime(d, 10)) return pollard_rho(d);
    return d;
 
}
 
vector<unsigned> factorize(const unsigned &n) {
 
    static map<unsigned, vector<unsigned>> factorizations_cache;
 
    vector<unsigned> factorizations;
    if (factorizations_cache.count(n)) {
        factorizations = factorizations_cache[n];
    } else {
        if (!(n % 2) && factorizations_cache.count(n / 2)) {
            factorizations = factorizations_cache[n / 2];
            factorizations.push_back(2);
        } else if (!(n % 3) && factorizations_cache.count(n / 3)) {
            factorizations = factorizations_cache[n / 3];
            factorizations.push_back(3);
        } else if (!(n % 5) && factorizations_cache.count(n / 5)) {
            factorizations = factorizations_cache[n / 5];
            factorizations.push_back(5);
        } else {
            unsigned changed_n = n;
            unsigned divisor = 1;
            while (divisor != changed_n) {
                changed_n /= divisor;
                divisor = pollard_rho(changed_n);
                factorizations.push_back(divisor);
            }
        }
        factorizations_cache[n] = factorizations;
    }
    return factorizations;
 
}
 
bool test_coprimality(const unsigned &c, const unsigned &m, const vector<unsigned int> &mFactor) {
 
    for (auto it = mFactor.rbegin(); it != mFactor.rend(); it++) {
        if (c % *it == 0) return false;
    }
    return true;
 
}
 
int main(int argc, char *argv[]) {
    
    get_input(cin);
    unsigned long long hull_dobell_passed = 0; // not greater than 2^40
    for (unsigned m = m_range.first; m <= m_range.second; m++) {
        coprimes_cache = vector<char>(c_range.second - c_range.first + 1, -1);
        if (d != 2 && is_prime(m, 5)) continue;
        if (!has_m_divisors(m, d)) { // M has at least D divisors
            continue;
        }
        bool m_is_divisible_by_4 = m % 4 == 0; // if 4 divides M then also 4 divides A-1
        unsigned a_increment = m_is_divisible_by_4 ? 4 : 1;
        for (unsigned a = (m_is_divisible_by_4 ? a_range_special_first : a_range.first); a <= a_range.second; a += a_increment) {
            vector<unsigned> factors_of_m = factorize(m);
            unsigned cCounter = 0;
            for (auto it = factors_of_m.rbegin(); it != factors_of_m.rend(); it++) {
                if ((a - 1) % *it != 0) goto continue_a;
            }
            for (unsigned c = c_range.first; c <= c_range.second; c++) {
                if (coprimes_cache[cCounter] == -1) {
                    bool res = test_coprimality(c, m, factors_of_m);
                    if (res) {
                        hull_dobell_passed++;
                    }
                    coprimes_cache[cCounter] = res ? 1 : 0;
                } else if (coprimes_cache[cCounter]) {
                    hull_dobell_passed++;
                }
                cCounter++;
            }
            continue_a:;
        }
    }
    cout << hull_dobell_passed << endl;
    return 0;
 
}
Добавлено через 1 час 41 минуту
Попытался изменить способ проверки взаимной простоты на проверку по теореме Эйлера. Чтобы добиться уникальности значений множителей, сменил вектор, в который сохранял результат факторизации, на set. В результате, скорость упала.
Добавил в начало заполнение сита Эратосфена и в проверку по нему в начале функции is_prime. При размере сита ~100 программа стала незначительно быстрее. При меньших и бóльших значениях скорость падает

Добавлено через 1 час 31 минуту
Прошу прощения, только сейчас заметил, что цикл по значениям C не зависит от значений А и его можно вынести за его пределы. Тем не менее, всё ещё слишком долго исполняется
0

Помощь в написании контрольных, курсовых и дипломных работ здесь.

Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
30.11.2016, 05:03
Ответы с готовыми решениями:

Параллельный генератор псевдослучайных чисел, на базе четырех линейных конгруэнтных генераторов
Всем доброе время суток)) У меня такая задачка: Нужно реализовать параллельный генератор...

Оптимизация в пределах линейных участков - исключение лишних переменных
Необходимо произвести оптимизацию в пределах линейных участков (исключить лишние переменные). На...

Алгоритмы поиска в линейных структурах
На железнодорожном вокзале хранится информация о поездах на текущие сутки (№ поезда, время...

Алгоритмы поиска в линейных структурах
На молочной ферме содержится информация о коров (номер, дата рождения), о доярок (номер, ФИО) и о...

1
0 / 0 / 0
Регистрация: 15.03.2012
Сообщений: 5
30.11.2016, 05:03  [ТС] 2
Прошу прощения, только сейчас заметил, что цикл по значениям C не зависит от значений А и его можно вынести за его пределы. Тем не менее, всё ещё слишком долго исполняется
0
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
30.11.2016, 05:03

Помощь в написании контрольных, курсовых и дипломных работ здесь.

Оптимизация поиска в листе
Есть задача Дан лист целых числе и число, написать функцию, возвращающую лист из двух значений в...

Оптимизация функции поиска
Есть функция поиска номеров по шаблону, шаблон задает пользователь. Пример ввода: 097???????. Надо...

Оптимизация поиска файлов
Как можно ускорить поиск файлов по всей файловой системе? Думал сделать что-то типо индексации, но...

Оптимизация поиска в спрвочнике
Здравствуйте, форумчане! У меня есть несколько CSV файлов, где помимо прочего хранятся названия...

Оптимизация формы поиска/фильтрации
Всем привет, Есть форма со всеми закупками (~10строк) и есть стойкое желание производить поиск...

Оптимизация категоризированного поиска товаров
Добрый день. Столкнулся с проблемой скорости загрузки страниц на каталогах подарков...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
2
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.