Спасибо за статью, я раньше думал, что memory-mapped файлы не очень хорошо работают на Windows.
База на диске — хороший компромисс между скоростью работы и простотой реализации, но для некоторых задач она не нужна:
— Если у вас в основном последовательный доступ, то есть, вам надо перебрать простые числа на интервале [X… Y], то имеет смысл их каждый раз заново просеивать: достаточно заранее в памяти просчитать и держать список простых чисел до sqrt(N), в вашем случае — до миллиона.
— Если у вас полностью рандомный доступ, то выгоднее воспользоваться критерием простоты чисел: en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases. Аккуратно написанный, он дает скорость 20 микросекунд на число (в отличие от 2 миллисекунд на seek на HDD): ceur-ws.org/Vol-1326/020-Forisek.pdf
Альтернативно, базу можно сократить в два раза, если хранить биты только на каждое нечетное число, или даже в три раза, если хранить только числа, дающие при делении на 6 остатки 1 или 5 (остальные не могут быть простыми, кроме чисел 2 и 3). А можно вообще упороться, и хранить только числа, дающие один из восьми остатков от деления на 30, что как раз занимает один байт. Я это все описывал в habr.com/ru/post/526924, там есть код на C#, который можно использовать для этого.
вот тут я вижу как достигают скорости 2.5 миллиона evaluated cells per second на CPU
все-таки там 2.5 миллиарда клеток в секунду, а не миллиона.
А с использованием SSE-команд можно еще быстрее, до 15 миллиардов в секунду: habr.com/ru/post/505606
А в этой статье показывают, что на GPU можно достичь скорости, если я не ошибся в подсчетах, полтора триллиона клеток в секунду!
примерно после девяти миллиардов итераций и двенадцати часов работы алгоритм нашел score=8 ݻටຈݹ௫ටຯѧ௨پଈݶಏقງݹலටະऄ௨پໃݢ௨ගɑݸҰقະࡨ௨sໃɥ௨ڏІݟ௨ගɱݸతටࢩݹࢴටแɕ௧ڄໃףۑඤÐݸలقลڝ௨ටໃݹ௨ටࠋ
Отличная задача на поиск оптимальной стратегии! Основная сложность — для каждой позиции определение следующей фигурки требует перебора всех фигурок по всем смещениям и вращениям, а количество позиций растет экспоненциально с увеличением глубины колодца. Используя очень эффеткивный алгоритм (перебирающий около двухсот тысяч позиций в секунду), мне удалось решить задачу для глубины 5 (score=4, всего пришлось перебрать около четырёхсот миллионов позиций за 30 минут): ޛටדݹ௨ඞТஈ௧ƣເঅ௩ɦܘݶహइໂђ௨ටໞݹஜɦܘݚذචʅݹਐටࢺݹસقໃݹ௨ටໃݹਆ
И я в процессе решения задачи для глубины 6. Пока что он нашел score=6 (поле перебора трех миллиардов позиций) ܮටลຍ௨ಀඍݸహටາɝɷචŒݹமϺາŦ௨ಀІݟௐقГঅஶѼໃӌ௨ගɱݸƓΟເח௧ڈໃӻ௨ගȣݸಏටງђஶɈժށѻටໃݹ௨ටເ0
но алгоритм еще не закончил работу, пока перебрав восемь миллиардов позиций за 11 часов.
Я попробовал еще больше соптимизировать однопоточную версию, используя вместо буфера на 15 чисел три буфера на 100 чисел, пользуясь тем, что остаток от деления на 15 у чисел n и n + 300 совпадают. Поэтому, если есть буфер чисел, например, [10000..10099], то от него можно перейти к буферу с числами [10300..10399] только поменяв цифру 0 на 3 в позициях, соответствующих числам, не делящимся на 3 и 5. Потом от нее можно перейти к [10600...10699] и т.д. При переполнении младшего разряда предыдущий увеличивается на единичку (получился алгоритм сложения с числом 300 в столбик). Например, переход от [19900..19999] к [20200..20299] получается заменой 9 -> 2 на третей позиции, 9 -> 0 на второй, и 1 -> 2 на первой для всех чисел.
На моем компьютере (у меня Windows/cygwin, поэтому мне сложно сравнить с вашими результатами напрямую) этот алгоритм работает 1.7 секунд. Я Вам послал пул реквест, было бы интересно узнать, будет ли он обгонять reusebuffer в ваших тестах. Правда, сейчас он работает только для LIMIT = 10^k, но, думаю, его можно доделать для любого диапазона чисел без ухудшения производительности.
decimal.c
#include <stdio.h>
#include <string.h>
#define LIMIT 1000000000
char *FIZZ = "Fizz";
char *BUZZ = "Buzz";
struct Segment {
char buffer[4096]; // 100 numbers
int buffer_length;
int positions[100]; // offsets in buffer1 of third digit of all numbers
int positions_length; // number of items in positions
};
int get_number_length(int number) {
int number_length = 1;
while (number >= 10) {
number /= 10;
number_length++;
}
return number_length;
}
char* print_number(int number, int numberLength, char* buf) {
for(int i = 0; i < numberLength; i++) {
int digit = number % 10;
buf[numberLength - i - 1] = (digit + '0');
number /= 10;
}
buf[numberLength] = '\n';
return buf + numberLength + 1;
}
char* print_fizz(char* buf) {
*(int*)buf = *(int*)FIZZ;
buf[4] = '\n';
return buf + 5;
}
char *print_buzz(char *buf) {
*(int *)buf = *(int *)BUZZ;
buf[4] = '\n';
return buf + 5;
}
char *print_fizz_buzz(char *buf) {
*(int *)buf = *(int *)FIZZ;
*(int *)(buf + 4) = *(int *)BUZZ;
buf[8] = '\n';
return buf + 9;
}
char *print(int num, int num_length, char *ptr) {
if (num_length == -1)
num_length = get_number_length(num);
if (num % 15 == 0)
ptr = print_fizz_buzz(ptr);
else if (num % 3 == 0)
ptr = print_fizz(ptr);
else if (num % 5 == 0)
ptr = print_buzz(ptr);
else
ptr = print_number(num, num_length, ptr);
return ptr;
}
char *print_first_segment(char *ptr) {
for(int num = 1; num < 100; num++)
ptr = print(num, -1, ptr);
return ptr;
}
char *print_last_number(char *ptr) {
return print(LIMIT, -1, ptr);
}
void print_initial_segment(int first_number, struct Segment* segment) {
int digits = get_number_length(first_number);
segment->positions_length = 0;
char *ptr = segment->buffer;
for (int i = first_number; i < first_number + 100; i++) {
ptr = print(i, digits, ptr);
if (i % 3 != 0 && i % 5 != 0)
segment->positions[segment->positions_length++] = ptr - segment->buffer - 4;
}
segment->buffer_length = ptr - segment->buffer;
}
void print_initial_segments(int segment, struct Segment* s1, struct Segment* s2, struct Segment* s3) {
print_initial_segment(segment, s1);
print_initial_segment(segment + 100, s2);
print_initial_segment(segment + 200, s3);
}
void next(struct Segment* s) {
char c = s->buffer[s->positions[0]];
char nextC = c + 3;
if(c < '7') {
for (int i = 0; i < s->positions_length; i++)
s->buffer[s->positions[i]] = nextC;
} else if (s->buffer[s->positions[0] - 1] < '9') {
nextC -= 10;
int nextD = s->buffer[s->positions[0] - 1] + 1;
for (int i = 0; i < s->positions_length; i++) {
s->buffer[s->positions[i]] = nextC;
s->buffer[s->positions[i] - 1] = nextD;
}
} else {
nextC -= 10;
for (int i = 0; i < s->positions_length; i++) {
int pos = s->positions[i];
s->buffer[pos--] = nextC;
while(s->buffer[pos] == '9')
s->buffer[pos--] = '0';
s->buffer[pos]++;
}
}
}
int main(void) {
char* ptr;
struct Segment s1, s2, s3;
ptr = print_first_segment(s1.buffer);
fwrite(s1.buffer, 1, ptr - s1.buffer, stdout);
for (int segment = 100; segment < LIMIT; segment *= 10) {
print_initial_segments(segment, &s1, &s2, &s3);
fwrite(s1.buffer, 1, s1.buffer_length, stdout);
fwrite(s2.buffer, 1, s2.buffer_length, stdout);
fwrite(s3.buffer, 1, s3.buffer_length, stdout);
int repeat_count = 3 * (segment / 100);
for(int i = 1; i < repeat_count; i++) {
next(&s1);
next(&s2);
next(&s3);
fwrite(s1.buffer, 1, s1.buffer_length, stdout);
fwrite(s2.buffer, 1, s2.buffer_length, stdout);
fwrite(s3.buffer, 1, s3.buffer_length, stdout);
}
}
ptr = print_last_number(s1.buffer);
fwrite(s1.buffer, 1, ptr - s1.buffer, stdout);
return 0;
}
Спасибо, интересная статья! Хотел только заметить, что для вычисления угла наклона пути лучше использовать Math.Atan2 — он не требует обработки специальных случаев и должен сократить ваши 9 строчек до одной.
Да, Вы правы, я уже после написания статьи заменил логические операции на shuffle, получив три строчки кода вместо шести (причем, более понятные): AdvancedLifeExtensions.cs:73
Правда, ускорение, по моим подсчетам, получилось не очень большое, около 5%.
В его статье приведены оба варианта: в главе 3 «Prime sieves using irreducible binary quadratic forms» — аналогичный приведенному здесь и в википедии, и в главе 4 «Enumerating lattice points» — его оптимизация до O(n / log log n).
В любом случае, именно алгоритм, приведенный в публикации, в Википедии называется решетом Аткина, хотя многие, действительно, с этим не согласны.
К тому же, я в публикации упоминаю утилиту primegen, основанную на оптимизированном решете Аткина, которая в 30 раз быстрее моей реализации.
А, теперь понял суть ваших претензий. Мне кажется, это вопрос терминологии, какой именно алгоритм является «классическим» решетом Аткина, а какой — его оптимизированной версией. Мой вариант — такой же, как и в статье в Википедии, хотя в дискуссии к ней он тоже критикуется.
В решете Сундарама мы вычеркиваем числа i + j + 2ij для любых i, j (не только простых), для оставшихся чисел 2n+1 — будет простое число, не 2i+1. То есть, число 7 мы вычеркнем правильно, потому что 2*7+1=15 не является простым числом. У меня было написано ni, я исправил на n, надеюсь, так понятней будет
То же самое время. Console.WriteLine вызывается, только если новое расстояние больше уже найденного, и так как максимальное расстояние равно 540, будет вызываться очень редко (около 50 раз)
Я правильно понимаю, что задача состоит в том, чтобы для чисел A и M найти такие числа x и y, что x^y = A, 2 <= y <= M?
Задача очень похожа на задачи из Project Euler и, наверно, поэтому не самая лучшая для сравнения CPU и GPU, поскольку, как мне кажется, она решается меньше, чем за секунду на CPU (в один поток) для любых 64-битных A и M. Идея в том, что если A — полный квадрат, то это легко проверить, просто взяв Math.Sqrt, а если A является степенью 3 или больше какого-то числа, то у него будет простой делитель, не превосходящий (примерно) один миллион, и его можно найти простым перебором. После факторизации числа A найти нужные x и y достаточно тривиально (у будет равен наименьшему общему делителю всех степеней простых чисел, входящих в разложение A).
Если абстрагироваться от задачи, то действительно интересно, почему CUDA медленней — обычно такое бывает либо потому, что время работы GPU мало по сравнению с временем подготовки и передачи данных, либо потому, что данные кладутся не в оптимальный для них тип GPU-памяти.
А мне нравится «информационное» определение энтропии для вероятностного распределения:
E = -∑(p*log p), то есть, энтропия — мат. ожидание минус логарифма вероятности события, или, что то же самое, мат. ожидание количества информации, необходимого, чтобы записать событие.
Например, для игрального кубика события — выпадение чисел 1, 2, 3, 4, 5, 6 с равной вероятностью 1/6; количество инфромации, чтобы передать любое из этих событий -log(1/6) = log(6) ~ 2.6 бит. И энтропия игрального кубика будет
E = -∑(1/6*log(1/6)) ~ 2.6 бит.
Если же кубик у нас не честный, и, например, вероятность выпадения единицы — одна вторая, а остальных чисел — по одной десятой, то энтропия будет
E = -0.5 * log(0.5) — ∑(1/10*log(1/10)) ~ 2.2 бита, то есть, меньше, чем у честного кубика.
В прогрммистских терминах, энтропия даных — это ожидаемый размер зип-архива с этими данными (в предположении, что он сожмет их идеально)!
Если Вы имеете в виду Number theoretic transform, то, возможно, в обработке изображений его не используют. Но, насколько я знаю, это самый быстрый способ перемножить два многочлена с целыми коэффициентами или два целых числа произвольной разрядности, и это имеет приложения, например, в криптографии. И научные статьи о низкоуровневых оптимизациях этого метода регулярно появляются, например, Speeding up the Number Theoretic Transform for Faster Ideal Lattice-Based Cryptography.
Там несколько другая задача, так как считаются остатки не обязательно по модулю простого числа, но кроме практических задач есть еще и научные, и для решения некоторых задач приходится перебирать большие диапазоны чисел, поэтому
>что изменится, если Вы будете считать остатки в 10 раз медленнее или в 10 раз быстрее?
от этого может зависить, напишете вы научную статью через месяц или через год!
Доказано, что эта задача — NP-полная, en.wikipedia.org/wiki/Addition-chain_exponentiation:
База на диске — хороший компромисс между скоростью работы и простотой реализации, но для некоторых задач она не нужна:
— Если у вас в основном последовательный доступ, то есть, вам надо перебрать простые числа на интервале [X… Y], то имеет смысл их каждый раз заново просеивать: достаточно заранее в памяти просчитать и держать список простых чисел до sqrt(N), в вашем случае — до миллиона.
— Если у вас полностью рандомный доступ, то выгоднее воспользоваться критерием простоты чисел: en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases. Аккуратно написанный, он дает скорость 20 микросекунд на число (в отличие от 2 миллисекунд на seek на HDD): ceur-ws.org/Vol-1326/020-Forisek.pdf
Альтернативно, базу можно сократить в два раза, если хранить биты только на каждое нечетное число, или даже в три раза, если хранить только числа, дающие при делении на 6 остатки 1 или 5 (остальные не могут быть простыми, кроме чисел 2 и 3). А можно вообще упороться, и хранить только числа, дающие один из восьми остатков от деления на 30, что как раз занимает один байт. Я это все описывал в habr.com/ru/post/526924, там есть код на C#, который можно использовать для этого.
все-таки там 2.5 миллиарда клеток в секунду, а не миллиона.
А с использованием SSE-команд можно еще быстрее, до 15 миллиардов в секунду: habr.com/ru/post/505606
А в этой статье показывают, что на GPU можно достичь скорости, если я не ошибся в подсчетах, полтора триллиона клеток в секунду!
Не совсем: проверка числа на простоту, в отличие от факторизации, решается за полиномиальное время
ݻටຈݹ௫ටຯѧ௨پଈݶಏقງݹலටະऄ௨پໃݢ௨ගɑݸҰقະࡨ௨sໃɥ௨ڏІݟ௨ගɱݸతටࢩݹࢴටแɕ௧ڄໃףۑඤÐݸలقลڝ௨ටໃݹ௨ටࠋޛටדݹ௨ඞТஈ௧ƣເঅ௩ɦܘݶహइໂђ௨ටໞݹஜɦܘݚذචʅݹਐටࢺݹસقໃݹ௨ටໃݹਆИ я в процессе решения задачи для глубины 6. Пока что он нашел score=6 (поле перебора трех миллиардов позиций)
ܮටลຍ௨ಀඍݸహටາɝɷචŒݹமϺາŦ௨ಀІݟௐقГঅஶѼໃӌ௨ගɱݸƓΟເח௧ڈໃӻ௨ගȣݸಏටງђஶɈժށѻටໃݹ௨ටເ0но алгоритм еще не закончил работу, пока перебрав восемь миллиардов позиций за 11 часов.
На моем компьютере (у меня Windows/cygwin, поэтому мне сложно сравнить с вашими результатами напрямую) этот алгоритм работает 1.7 секунд. Я Вам послал пул реквест, было бы интересно узнать, будет ли он обгонять reusebuffer в ваших тестах. Правда, сейчас он работает только для LIMIT = 10^k, но, думаю, его можно доделать для любого диапазона чисел без ухудшения производительности.
Правда, ускорение, по моим подсчетам, получилось не очень большое, около 5%.
В любом случае, именно алгоритм, приведенный в публикации, в Википедии называется решетом Аткина, хотя многие, действительно, с этим не согласны.
К тому же, я в публикации упоминаю утилиту primegen, основанную на оптимизированном решете Аткина, которая в 30 раз быстрее моей реализации.
11 — 7 = 4
29 — 23 = 6
97 — 89 = 8
127 — 113 = 14
541 — 523 = 18
907 — 887 = 20
1151 — 1129 = 22
1361 — 1327 = 34
9587 — 9551 = 36
15727 — 15683 = 44
19661 — 19609 = 52
31469 — 31397 = 72
156007 — 155921 = 86
360749 — 360653 = 96
370373 — 370261 = 112
492227 — 492113 = 114
1349651 — 1349533 = 118
1357333 — 1357201 = 132
2010881 — 2010733 = 148
4652507 — 4652353 = 154
17051887 — 17051707 = 180
20831533 — 20831323 = 210
47326913 — 47326693 = 220
122164969 — 122164747 = 222
189695893 — 189695659 = 234
191913031 — 191912783 = 248
387096383 — 387096133 = 250
436273291 — 436273009 = 282
1294268779 — 1294268491 = 288
1453168433 — 1453168141 = 292
2300942869 — 2300942549 = 320
3842611109 — 3842610773 = 336
4302407713 — 4302407359 = 354
10726905041 — 10726904659 = 382
20678048681 — 20678048297 = 384
22367085353 — 22367084959 = 394
25056082543 — 25056082087 = 456
42652618807 — 42652618343 = 464
127976335139 — 127976334671 = 468
182226896713 — 182226896239 = 474
241160624629 — 241160624143 = 486
297501076289 — 297501075799 = 490
303371455741 — 303371455241 = 500
304599509051 — 304599508537 = 514
416608696337 — 416608695821 = 516
461690510543 — 461690510011 = 532
614487454057 — 614487453523 = 534
738832928467 — 738832927927 = 540
Processed primes up to 1,000,000,000,000 in 00:26:23.3155338: Max Distance=540
Задача очень похожа на задачи из Project Euler и, наверно, поэтому не самая лучшая для сравнения CPU и GPU, поскольку, как мне кажется, она решается меньше, чем за секунду на CPU (в один поток) для любых 64-битных A и M. Идея в том, что если A — полный квадрат, то это легко проверить, просто взяв Math.Sqrt, а если A является степенью 3 или больше какого-то числа, то у него будет простой делитель, не превосходящий (примерно) один миллион, и его можно найти простым перебором. После факторизации числа A найти нужные x и y достаточно тривиально (у будет равен наименьшему общему делителю всех степеней простых чисел, входящих в разложение A).
Если абстрагироваться от задачи, то действительно интересно, почему CUDA медленней — обычно такое бывает либо потому, что время работы GPU мало по сравнению с временем подготовки и передачи данных, либо потому, что данные кладутся не в оптимальный для них тип GPU-памяти.
E = -∑(p*log p), то есть, энтропия — мат. ожидание минус логарифма вероятности события, или, что то же самое, мат. ожидание количества информации, необходимого, чтобы записать событие.
Например, для игрального кубика события — выпадение чисел 1, 2, 3, 4, 5, 6 с равной вероятностью 1/6; количество инфромации, чтобы передать любое из этих событий -log(1/6) = log(6) ~ 2.6 бит. И энтропия игрального кубика будет
E = -∑(1/6*log(1/6)) ~ 2.6 бит.
Если же кубик у нас не честный, и, например, вероятность выпадения единицы — одна вторая, а остальных чисел — по одной десятой, то энтропия будет
E = -0.5 * log(0.5) — ∑(1/10*log(1/10)) ~ 2.2 бита, то есть, меньше, чем у честного кубика.
В прогрммистских терминах, энтропия даных — это ожидаемый размер зип-архива с этими данными (в предположении, что он сожмет их идеально)!
Там несколько другая задача, так как считаются остатки не обязательно по модулю простого числа, но кроме практических задач есть еще и научные, и для решения некоторых задач приходится перебирать большие диапазоны чисел, поэтому
>что изменится, если Вы будете считать остатки в 10 раз медленнее или в 10 раз быстрее?
от этого может зависить, напишете вы научную статью через месяц или через год!