Хабр Курсы для всех
РЕКЛАМА
Практикум, Хекслет, SkyPro, авторские курсы — собрали всех и попросили скидки. Осталось выбрать!
Хабр доступен 24/7 благодаря поддержке друзей
Комментарии 46
А не хотите графики совместить? было бы на много удобнее сравнивать
Вряд ли станет удобнее, посмотрите на шкалу обоих графиков
В том же и суть. Пока на первый взгляд алгоритмы не сильно отличаются, а если графики совместить, то сразу будет видна разница.
Добавил совмещенный график
Вот теперь — найс. Сразу понятна разница между алгоритмами(
В логарифмической шкале было бы нагляднее. А сколько будет работать второй алгоритм, если включить dump()?
Функцию dump я написал, чтобы убедиться, что второй алгоритм возвращает те же результаты, что и первый. Вывод там выполняется в шестнадцатеричном виде, что, в принципе, не должно занять много времени. Но сама запись данных на жесткий диск будет узким местом для второго алгоритма.
Да, если в 16-ричном то не интересно.
Получается, что 50-мегабайтные числа они перемножают за 5 секунд. Как?
Получается, что 50-мегабайтные числа они перемножают за 5 секунд. Как?
Это уже вопрос к разработчикам GMP.
Может быть. И по какому же модулю надо считать БПФ для такого размера чисел? Если брать байтовое представление — не меньше, чем 2^42. И сколько там будет операций вычисления остатка от деления 84-битных чисел на 42-битные? Долго получится.
Если по вашему модулю, то (100 * 1024 * 1024 * 8) / 42 = 19 972 876.2
Не понял. Цифры тоже надо возводить в квадрат. Если взять двоичную запись числа, а в качестве «цифры» рассматривать байты, то переполнения можно избежать, если взять P не меньше, чем 256^2*50*2^20= 3.4*10^12 — около 42 бит. При этом P-1 должно делиться на 2^27, чтобы модулярное БПФ могло работать. Впрочем, при запасе в 15 бит это не проблема, можно взять, например, P=15*2^38+1. Но как эффективно реализовать умножение по этому модулю?
Можно посмотреть исходники GMP
Очень странно. Они пишут
(файл mpn/generic/mul_fft.c).
Но ведь 2^(2*N) mod (2^N+1)=1, так что стерени двойки принимают не более 2*N значений. Откуда там примитивный элемент. Да еще и по непростому модулю.
Кто-нибудь в этом разобрался?
where 2^omega is a primitive root mod 2^N+1
(файл mpn/generic/mul_fft.c).
Но ведь 2^(2*N) mod (2^N+1)=1, так что стерени двойки принимают не более 2*N значений. Откуда там примитивный элемент. Да еще и по непростому модулю.
Кто-нибудь в этом разобрался?
Получается, что они вообще не используют модулярную арифметику над числами. Все вычисления идут с многочленами, коэффициенты которых являются чуть более короткими, но все равно длинными числами. И умножение работает рекурсивно. Все равно надо будет понять, как работает их «fast convolution» в этих условиях.
Умножение по константному модулю вряд ли выполняется медленно. Это моно сделать за 2 команды целочисленного умножения, которые на Core i7 и AMD K10 выполняются за 3 такта каждая
Как? Вот у нас есть
long long a,b;
const long long P=0x3c000000001ll;
Хочется найти (a*b)%P. Как это сделать за «2 команды целочисленного умножения»?
long long a,b;
const long long P=0x3c000000001ll;
Хочется найти (a*b)%P. Как это сделать за «2 команды целочисленного умножения»?
Сначала вычисляем a * b, а затем находим остаток от деления на константу. Деление на константу можно заменить умножением. В руководстве по оптимизации для процессоров AMD описан алгоритм такой замены.
support.amd.com/us/Processor_TechDocs/40546.pdf
раздел «8.1 Replacing Division with Multiplication»
support.amd.com/us/Processor_TechDocs/40546.pdf
раздел «8.1 Replacing Division with Multiplication»
по сути это замена N / C = N * (1 / C). Поскольку C — константа, то 1 / C можно вычислить заранее
Хотя в том же руководстве написано, что для семейства 12h такая замена неактуальна, т.к. деление и так выполняется быстро
А там есть остаток от деления 128 бит на 64?
В x86/x64 на вход операции деления подается частное 2 регистрах и делитель в одном регистре. Для x64 как раз получаем 128-битное частное и 64-битный делитель
При использовании деления (а оно сейчас быстрое, как выяснилось) получается вот такой код умножения по модулю 64-битных чисел:
; a = RAX
; b = RBX
; c = RCX
MUL RBX
DIV RCX
; в RDX у нас модуль (a * b) % c
a*b будет уже 128-битным. Деление заменяется на умножение, чтобы найти частное, и то с ограниченной точностью. Чтобы найти остаток, придется умножить P*d — это будет уже третье умножение (а еще перед нимпотребуется сдвиг 128-битного числа, чтобы получить нужные биты). И потом, возможно, надо будет подкорректировать результат. Наверное, за разумное время это можно сделать.
Вот только в GMP делают не так, они берут модуль 2^N+1 — который не простой, а в качестве omega используют степень двойки. Как оно может работать, пока не понимаю.
Вот только в GMP делают не так, они берут модуль 2^N+1 — который не простой, а в качестве omega используют степень двойки. Как оно может работать, пока не понимаю.
Логарифмический масштаб:
habrastorage.org/storage2/67a/b64/5e6/67ab645e6b8a42f956b4d1b3d36500fb.png
habrastorage.org/storage2/67a/b64/5e6/67ab645e6b8a42f956b4d1b3d36500fb.png
Они построены на разных данных, т.к. там где алгоритм №1 в состоянии посчитать число Фибоначчи за разумное количество времени, алгоритм №2 отрабатывает менее, чем за секунду. Если вы посмотрите на таблицу, то увидите, что самое большое число, обработанное алгоритмом №1, равно 3.000.000, а самое маленькое, обработанное алгоритмом №2, — 50.000.000. При том, что число в 16 раз больше, затрачиваемое время — в 50 раз меньше.
Неделя чисел Фибоначчи на хабре?
Респект автору. Тут было много измышлений по поводу алгоритмов вычислений чисел Фибоначчи, а это топик с правильной реализацией (к сожалению в Python не встроены быстрые алгоритмы для умножения длинных чисел). Конечно есть некоторые недочеты, например, как уже указывали, что нужен логарифмический масштаб для графиков.
Конечно, общественное сознание очень сильно волнует вычисленное количество знаков после запятой числа Π или быстрое вычисление чисел Фибоначчи. Я, правда, лично не знаю для чего общественное сознание так сильно ВОЛНУЕТ эти задачи. Может его заинтересует и эта тривиальная задача. Сколько корней имеет следующая полиномиальная система и если сможите найдите все ее корни системы
x^3 — y^2 + z — 1 = 0
y^3 — z^2 + x — 1 = 0
z^3 — x^2 + y — 1 = 0
Конечно, общественное сознание очень сильно волнует вычисленное количество знаков после запятой числа Π или быстрое вычисление чисел Фибоначчи. Я, правда, лично не знаю для чего общественное сознание так сильно ВОЛНУЕТ эти задачи. Может его заинтересует и эта тривиальная задача. Сколько корней имеет следующая полиномиальная система и если сможите найдите все ее корни системы
x^3 — y^2 + z — 1 = 0
y^3 — z^2 + x — 1 = 0
z^3 — x^2 + y — 1 = 0
А использовать Maple разрешается?
В ответе по-моему нет 27 корней.
Зато есть кнопочка «More solutions».
Я ее и нажимал, увеличилось но все равно мало.
Странно, у меня показывает 26 комплексных и одно целое решение
habrastorage.org/storage2/786/9e8/7a7/7869e87a7227a907c09f73d8434707b0.png
habrastorage.org/storage2/786/9e8/7a7/7869e87a7227a907c09f73d8434707b0.png
Мне по more solutions не удалось получить больше 21 корня. Мало жал, наверное.
Мой ответ.
Эта программа получает следующий эквивалентный вид системы
z^27 + (-9)*z^24 + 29*z^21 + 6*z^19 + (-53)*z^18 + 22*z^17 + (-63)*z^16 + 96*z^15 + (-149)*z^14 + 242*z^13 + (-261)*z^12 + 484*z^11 + (-545)*z^10 + 740*z^9 + (-908)*z^8 + 972*z^7 + (-1220)*z^6 + 1047*z^5 + (-1045)*z^4 + 943*z^3 + (-535)*z^2 + 422*z + (-216)*1
userTime = 0.0
sysTime = 0.01
realTime = 0.01
hilbertPolynomial = 27
Значит система имеет 27 корней с учетом кратности. Их можно теперь легко найти, найдя корни полинома от z, например с помощью Maxima
to_poly_solve([z^27 + (-9)*z^24 + 29*z^21 + 6*z^19 + (-53)*z^18 + 22*z^17 + (-63)*z^16 + 96*z^15 + (-149)*z^14 + 242*z^13 + (-261)*z^12 + 484*z^11 + (-545)*z^10 + 740*z^9 + (-908)*z^8 + 972*z^7 + (-1220)*z^6 + 1047*z^5 + (-1045)*z^4 + 943*z^3 + (-535)*z^2 + 422*z + (-216)*1], [z]);
Программа на Python
# -*- coding: utf-8 -*-
import ginv
st = ginv.SystemType("Polynomial")
im = ginv.MonomInterface("Lex", st, ['x', 'y', 'z'])
eqs = [
'x^3 - y^2 + z - 1',
'y^3 - z^2 + x - 1',
'z^3 - x^2 + y - 1',
]
ic = ginv.CoeffInterface("GmpZ", st)
ip = ginv.PolyInterface("PolyList", st, im, ic)
iw = ginv.WrapInterface("CritPartially", ip)
id = ginv.DivisionInterface("Janet", iw)
basis = ginv.basisBuild("TQ", id, eqs)
for p in basis.iterGB():
print p
print " userTime = ", basis.userTime()
print " sysTime = ", basis.sysTime()
print " realTime = ", basis.realTime()
print "hilbertPolynomial = ", basis.hilbertPolynomial()
Эта программа получает следующий эквивалентный вид системы
1395367452523974847088496*x +
(-36974012043720606602747)*z^26 + (-18042387149981405949931)*z^25 + 33030642427988789483293*z^24 + 290848876458135913225112*z^23 + 129076392042710083597936*z^22 + (-219952910543490007651408)*z^21 + (-735043497452966299206159)*z^20 + (-310272794617140912662271)*z^19 + 160633475001590756076807*z^18 + 982342916167428059233558*z^17 + (-20920432695844226462772)*z^16 + 1237879512449968788217761*z^15 + (-649364739492104109294111)*z^14 + 1770783798251395989600640*z^13 + (-3269570785655865115267646)*z^12 + 1017137392573548569160425*z^11 + (-5380290668716432422480115)*z^10 + 2511318414925313271358736*z^9 + (-5820765926743222198514292)*z^8 + 3916234105006809425362144*z^7 + (-1271740986960481649374276)*z^6 + 4526571866128344763445768*z^5 + 4225399407124480064647579*z^4 + (-1702488139911346623054830)*z^3 + 2599321588208462787828165*z^2 + (-6499050198343010356877518)*z + 1138531641035453825071656*1
1395367452523974847088496*y +
(-24263245905633422666089)*z^26 + (-2564101797577402706105)*z^25 + (-7612716934398042952513)*z^24 + 183708473843508669676888*z^23 + 27382590493064706561728*z^22 + 93382883659052502917872*z^21 + (-445019735963407736728965)*z^20 + (-124580325693827052967653)*z^19 + (-496913514381305731286115)*z^18 + 686028568407960758274770*z^17 + (-293264534395468068345180)*z^16 + 1754916317739869655913083*z^15 + (-1494161044420719326758149)*z^14 + 3166864848665922491507168*z^13 + (-4324148866718107366719898)*z^12 + 5118193353820594354340467*z^11 + (-9669988041948198759925505)*z^10 + 8234964392530939187961184*z^9 + (-14800290986146424015670204)*z^8 + 14407709710372641542274656*z^7 + (-17110238558169309305535500)*z^6 + 19289855195567445653291224*z^5 + (-14071675248685297091156215)*z^4 + 17208495316564910541540422*z^3 + (-9780261032315595871168233)*z^2 + 5684543607089979672443254*z + (-4606430757804595389204888)*1
z^27 + (-9)*z^24 + 29*z^21 + 6*z^19 + (-53)*z^18 + 22*z^17 + (-63)*z^16 + 96*z^15 + (-149)*z^14 + 242*z^13 + (-261)*z^12 + 484*z^11 + (-545)*z^10 + 740*z^9 + (-908)*z^8 + 972*z^7 + (-1220)*z^6 + 1047*z^5 + (-1045)*z^4 + 943*z^3 + (-535)*z^2 + 422*z + (-216)*1
userTime = 0.0
sysTime = 0.01
realTime = 0.01
hilbertPolynomial = 27
Значит система имеет 27 корней с учетом кратности. Их можно теперь легко найти, найдя корни полинома от z, например с помощью Maxima
to_poly_solve([z^27 + (-9)*z^24 + 29*z^21 + 6*z^19 + (-53)*z^18 + 22*z^17 + (-63)*z^16 + 96*z^15 + (-149)*z^14 + 242*z^13 + (-261)*z^12 + 484*z^11 + (-545)*z^10 + 740*z^9 + (-908)*z^8 + 972*z^7 + (-1220)*z^6 + 1047*z^5 + (-1045)*z^4 + 943*z^3 + (-535)*z^2 + 422*z + (-216)*1], [z]);
Ответ
%union([z=1],[z=-1.433140352978328*%i-0.87723099406978],[z=-1.180643959560538*%i-
0.30394363865802],[z=-1.095904638315721*%i-0.43038640721464],[z=-1.092531540289019*%i-
0.52912398837362],[z=-1.031468814711028*%i-0.71140740529994],[z=0.39699714792632-
0.99846864859532*%i],[z=0.71315631594634-0.85827126307494*%i],[z=-0.76039132493366*%i-
1.028797119340559],[z=1.29963755259568-0.47886775883786*%i],[z=-0.44847182216701*%i-
1.384073601510903],[z=0.82323873561898-0.1943816172483*%i],[z=1.531933402380145-
0.027231818844858*%i],[z=0.027231818844858*%i+1.531933402380145],[z=0.1943816172483*%i+
0.82323873561898],[z=0.44847182216701*%i-1.384073601510903],[z=0.47886775883786*%i+
1.29963755259568],[z=0.76039132493366*%i-1.028797119340559],[z=0.85827126307494*%i+
0.71315631594634],[z=0.99846864859532*%i+0.39699714792632],[z=-%i],[z=%i],[z=
1.031468814711028*%i-0.71140740529994],[z=1.092531540289019*%i-0.52912398837362],[z=
1.095904638315721*%i-0.43038640721464],[z=1.180643959560538*%i-0.30394363865802],[z=
1.433140352978328*%i-0.87723099406978])
0.30394363865802],[z=-1.095904638315721*%i-0.43038640721464],[z=-1.092531540289019*%i-
0.52912398837362],[z=-1.031468814711028*%i-0.71140740529994],[z=0.39699714792632-
0.99846864859532*%i],[z=0.71315631594634-0.85827126307494*%i],[z=-0.76039132493366*%i-
1.028797119340559],[z=1.29963755259568-0.47886775883786*%i],[z=-0.44847182216701*%i-
1.384073601510903],[z=0.82323873561898-0.1943816172483*%i],[z=1.531933402380145-
0.027231818844858*%i],[z=0.027231818844858*%i+1.531933402380145],[z=0.1943816172483*%i+
0.82323873561898],[z=0.44847182216701*%i-1.384073601510903],[z=0.47886775883786*%i+
1.29963755259568],[z=0.76039132493366*%i-1.028797119340559],[z=0.85827126307494*%i+
0.71315631594634],[z=0.99846864859532*%i+0.39699714792632],[z=-%i],[z=%i],[z=
1.031468814711028*%i-0.71140740529994],[z=1.092531540289019*%i-0.52912398837362],[z=
1.095904638315721*%i-0.43038640721464],[z=1.180643959560538*%i-0.30394363865802],[z=
1.433140352978328*%i-0.87723099406978])
Если решать с нуля, то проще использовать метод гомотопий. То, что корней 27 видно, поскольку система уже является базисом гребнера. Если коэффициенты при 1 и 2 ступенях обнулить, а свободные члены чуть изменить, сделав комплексными, корни будет найти легко. Потом достаточно плавно двигать коэффициенты к правильным, а на каждом шаге уточнять корни методом ньютона. Главное чтобы траектории не сливались. Получается 27 корней.
Согласен, только базис Грёбнера в таком-то допустимом относительно умножения линейного порядка. В данном случае в любом порядке сначала по полной степени, поскольку, тогда лидирующие мономы взаимно просты, то согласно второму критерию Бухбергера все S-полиномы будут редуцироваться к нулю и значит это базис Грёбнера. Число мономов не имеющих делителей среди лидирующих мономов {x^3, y^3, z^3} будет в точности 27=#{1_M, x, y, z, x^2,…, x^2y^2z^2} и значит полином Гильберта равен 27 и именно столько корней с учетом кратности имеет система. А корни после построения базиса Грёбнера в лексикографическом порядке можно найти и методами отделения корней. Также можно действовать и по другому, перейти к рассмотрению фактора кольца по идеалу и к коммутативному кольцу матриц, представляющих элементы данного фактора. Тогда вопрос о числе различных корней и числе различных действительных корней решается с помощью аппарата Trace-матриц без их явного нахождения. А собственные значения матриц соответствующих переменных будут давать их корни.
А можно не строить базис в лексикографическом базисе, а прямо в текущем найти редукции x^3, x^4,..., x^27. Линейная зависимость результатов редукции и даст минимальный полином для x. Аналогично ищутся выражения для y и z (как линейная зависимость y,1,x,x^2,...,x^26).
Метод гомотопий хорош тем, что в нем не надо искать базис Гребнера вообще — достаточно поиграться с многогранниками Ньютона и суммами Минковского. Сразу находятся все корни исходной системы (включая кратные, нулевые и бесконечные). Но вернуть результат в алгебраическую форму может не получиться (хотя, найдем (x-x1)(x-x2)...(x-xk) — многочлен готов. Не обязательно минимальный).
Метод гомотопий хорош тем, что в нем не надо искать базис Гребнера вообще — достаточно поиграться с многогранниками Ньютона и суммами Минковского. Сразу находятся все корни исходной системы (включая кратные, нулевые и бесконечные). Но вернуть результат в алгебраическую форму может не получиться (хотя, найдем (x-x1)(x-x2)...(x-xk) — многочлен готов. Не обязательно минимальный).
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Сравнение алгоритмов вычисления чисел Фибоначчи