CreatorCray (тоже не бог весть какой источник, но всё же)) В общем, на одном сайте олимпиадных задач есть задача "умножение длинных чисел". Длина там 10000 цифр. Вот, моё решение карацубой работает раз в 10 медленнее пары решений, и они как раз вроде бы писали БПФ. Так что видимо БПФ начинает обыгрывать алгоритм Карацубы уже при меньших размерах.
maxdiver ну, все зависит от имплементации. Может у тебя в карацубе неверно была подобрана минимальная длина множителей, при котором уже не происходит дальнейшая рекурсия. У меня в моей реализации такое уже было - я оставил эту константу с дебаговых времен и в результате дробил числа до размера в 32 бита. Сейчас до 768 бит. Разница в скорости весьма впечатляющая. 10000 десятичных цифр это примерно 33219 бит. т.е. чуть больше чем 2^(2^15). Что в общем то как раз лежит в указанных на вики пределах.
Perec В КОМ все по немецки???? Кого скачать? Описание алгоритма смотри во втором томе Кнута (английский или русский, какой найдешь) или тут http://en.wikipedia.org/wiki/Sch%C3%B6nhage-Strassen_algorithm - там все на англицком... Где взять имплементацию алгоритма - без понятия. Попробуй поискать на http://koders.com/ или http://www.google.com/codesearch?hl=ru
CreatorCray У Кнута как такового нет описание Шенхаге. Там есть "Китайская теорема". Еще посмотрю. Может пропустил? Что врядли Но все равно спасибо что время на ответы потратили!!!! Если что появится сообщи пожалуйста!!!
Что-то у меня не получается написать аналог преобразования Фурье в поле целых чисел Не догоняю что-то, как именно там надо вводить модули и каким образом выбирать p, по которому и берется mod p (например, в Кормене предлагается выбирать p наименьшее среди простых вида K*N+1, где N - как я понял длина). Ни у кого нет исходника? В гугле пытался искать, что-то не получилось найти Заранее спасибо.
http://rapidshare.de/files/38537380/fftmult.rar.html Забираем быстрое умножение на Дельфах с применением FFTW. Сама библиотека проводит только БПФ, а вот умножение - это уже целиком моя заслуга. Там же есть и умножение в столбик для проверки. Важное замечание - замо FFT работает с плавучкой, но сам результат умножения всегда точный, если ошибка округления не превышает 0.4 : В противном случае код отрубается. Код сто пудов рабочий, с помошью него я находил офигенные простые числа и некоторые даже запостил в соответствующие анналы.
Ну это и я могу выложить умножение, использующее double'ы. Пожалуй я так и сделаю Код (Text): typedef complex<double> base; void multiply (const vector<int> & a, const vector<int> & b, vector<int> & res); void fft (vector<base> & a, bool back); int rev (int n, int bits); void multiply (const vector<int> & a, const vector<int> & b, vector<int> & res) { int n = 1; while (n < (int) max (a.size(), b.size()) ) n <<= 1; n <<= 1; vector<base> aa (n), bb (n); copy (a.begin(), a.end(), aa.begin()); copy (b.begin(), b.end(), bb.begin()); fft (aa, false); fft (bb, false); vector<base> cc (n); for (int i=0; i<n; ++i) cc[i] = aa[i] * bb[i]; fft (cc, true); res.resize (n); for (int i=0, carry=0; i<n || carry; ++i) { if (i == n) { res.push_back (carry); break; } res[i] = carry + (int) floor (cc[i].real() + 0.0001); carry = res[i] / 10; res[i] %= 10; } while (res.size() > 1 && res.back() == 0) res.pop_back(); } void fft (vector<base> & a, bool back) { int n = (int) a.size(); int lg_n = 0; for (int i=1; i<n; i<<=1, ++lg_n) ; vector<char> swapped (n); for (int i=0; i<n; ++i) if (!swapped[i]) { int j = rev (i, lg_n); swap (a[i], a[j]); swapped[i] = swapped[j] = true; } for (int s=1; s<=lg_n; ++s) { int m = 1 << s; base wm (cos (M_PI*2/m), sin (M_PI*2/m)); if (back) wm.imag (-wm.imag()); for (int k=0; k<n; k+=m) { base w (1); for (int j=0; j<m/2; ++j) { base t = w * a[k+j+m/2]; base u = a[k+j]; a[k+j] = u + t; a[k+j+m/2] = u - t; if (back) a[k+j] /= 2, a[k+j+m/2] /= 2; w *= wm; } } } } int rev (int n, int bits) { int res = 0; for (int i=0; i<bits; ++i) res |= ((n >> i) & 1) << (bits-i-1); return res; } А вот с целочисленной арифметикой что-то никак не получается. Просто, если я правильно понимаю, с ней должно получиться заметно быстрее...
Везде написано, что целочисленные Фурье всегда медленнее работают, их применяют чтобы уж совсем исключить возможность ошибки из-за округления. По моему опыту - это все фигня, ведь в double можно класть хоть 16 бит, хоть 8-бит. Если совсем припрет и числа уж очень огромные, можно в даблы класть хоть по одному биту и всегда получить желаемую погрешность. Я бы не стал вообще с целочисленным Фурье заморачиваться, все рекордные простые числа включая мерсена с mersenne.org находят именно плавучкой.
Штрассен работает быстрее плавучки только на числах в миллиарды и сотни миллиардов цифр, где пришлось бы набивать даблы по нескольку а то и по одному биту, или даже одного бита было бы слишком много, чтобы ошибки округления все испортили. Это относится только к сверх-сверх большим числам. А вот для обычных больших чисел, типа не открытого еще вожделенного числа Мерсена на 10 миллионов знаков плавучки более чем достаточно.
Вы не ответили на вопрос. Мне интересны более точные оценки количества цифр, при которых один метод превосходит другой.
Ну блин вы как профессор на экзамене... Если сами все знаете, то так и напишите, а мы почитаем лекционный материал =)))
Факт тот, что целочисленное фурье было придумано не просто-так. Но вот выгодно ли его использовать конкретно для целочисленного умножения - я не знаю. Дело в том, что вещественное фурье может быть ускорено в 2 раза, если есть 2 чисто действительных набора входных данных (нужно просто прибавить к первому набору второй набор, домноженный на i, а затем можно будет разделить этот совместный образ на 2). Правда, при целочисленном умножений нужно 3 преобразования фурье, поэтому это дает ускорение лишь в 4/3 раза. В остальном же, целочисленное Фурье по-моему выигрывает за счет более простых операций. Особенно, если использовать ассемблерное деление int64 на int32. А авторы же программ для поиска чисел Мерсенна, скорее всего, никогда о целочисленном Фурье не слыхали и ассемблера не знают. Так что без проведения опыта нельзя сказать, какой алгоритм эффективнее.
а примеры этих алгоритмов есть на ассемблере? или мир к ним еще не готов, и их не успели написать? а то алгоритмы самые быстрые, ассемблер самый быстрый тоже, а вот вместе что то у них не срастаеться
Прошло больше года, а я таки взялся снова за этот алгоритм и таки разобрался, что из себя представляет целочисленная версия И ведь действительно она медленней оказалась Причём в разы. Не, ну если её заоптимайзить по самое не могу, вполне возможно, всё-таки целочисленное окажется быстрее. Там очень много операций взятия по модулю, а вроде как есть метод (Монтгомери, если я не ошибаюсь), с помощью которого можно будет от них избавиться. Может, кто-нибудь займётся этим. Для себя я понял, что по сочетанию длина кода <=> скорость работы плавучное ДПФ рулит Код (Text): const int mod = 7340033; const int root = 5; const int root_1 = reverse (root); const int root_pw = 1<<20; void fft (vector<int> & a, bool invert) { int n = (int) a.size(); for (int i=1, j=0; i<n; ++i) { int bit = n >> 1; for (; j>=bit; bit>>=1) j -= bit; j += bit; if (i < j) swap (a[i], a[j]); } for (int len=2; len<=n; len<<=1) { int wlen = invert ? root_1 : root; for (int i=len; i<root_pw; i<<=1) wlen = int (wlen * 1ll * wlen % mod); for (int i=0; i<n; i+=len) { int w = 1; for (int j=0; j<len/2; ++j) { int u = a[i+j], v = int (a[i+j+len/2] * 1ll * w % mod); a[i+j] = u+v < mod ? u+v : u+v-mod; a[i+j+len/2] = u-v >= 0 ? u-v : u-v+mod; w = int (w * 1ll * wlen % mod); } } } if (invert) { int nrev = reverse (n); for (int i=0; i<n; ++i) a[i] = int (a[i] * 1ll * nrev % mod); } } int powmod (int a, int b) { int res = 1; while (b) if (b & 1) { res = int (res * 1ll * a % mod); --b; } else { a = int (a * 1ll * a % mod); b >>= 1; } return res; } int reverse (int n) { return powmod (n, mod-2); } Здесь root - примитивное значение корня root_pw-ой степени из единицы по модулю mod. Искал подбором, на практике обычно этого значения root_pw=2^20 хватает (это ограничение на макс. длину вектора), mod=2^20 * 7 + 1 ограничивает значения, принимаемые каждым элементом. Как всегда, предполагается, что входной вектор имеет длину - степень двойки. На всякий случай даю реализацию ДПФ на комплексных числах, с которой я и сравнивал целочисленный вариант. Код (Text): const double PI = 3.1415926535897932384626433832795; typedef complex<double> base; void fft (vector<base> & a, bool invert) { int n = (int) a.size(); for (int i=1, j=0; i<n; ++i) { int bit = n >> 1; for (; j>=bit; bit>>=1) j -= bit; j += bit; if (i < j) swap (a[i], a[j]); } for (int len=2; len<=n; len<<=1) { double ang = 2*PI/len * (invert ? -1 : 1); base wlen (cos(ang), sin(ang)); for (int i=0; i<n; i+=len) { base w (1); for (int j=0; j<len/2; ++j) { base u = a[i+j], v = a[i+j+len/2] * w; a[i+j] = u + v; a[i+j+len/2] = u - v; w *= wlen; } } } if (invert) for (int i=0; i<n; ++i) a[i] /= n; } Итого: на числах длины 10000 целочисленное работало максимум 199 мс, комплексное - 73 мс (имеются в виду все затраты на чтение двух чисел, 3 ДПФ для вычисления умножения и вывод произведения).
