Версия для печати темы Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Алгоритмы > Быстрое преобразование Фурье

Автор: T0ohtik 28.4.2008, 23:19

Помогите разобраться с этим алгоритмом. В гугл просьба не слать, был не помог. Только формулы... Примеров нет. Задача такова, что есть очень большое число порядка 2 тыс разрядов. Надо перемножать на другое число. Для решения этой задачи я выбрал алгортм Быстрого преобразования Фурье. К примеру надо перемножить вектор А - 6 разрядов(обозначим m) и вектор B - 5 разрядов(обозначим n). Для начала надо посчитать коэфициенты С. Сколько будит их? (max(m,n))  ?

Автор: ksili 29.4.2008, 06:05

Что за бред?

Автор: maxim1000 29.4.2008, 10:34

Цитата(ksili @  29.4.2008,  06:05 ) Что за бред? многие, казалось бы, далёкие области математики иногда соприкасаются в неожиданных местах ;) умножение длинных чисел, если мне не изменяет память, сводится к свёртке, а БПФ может помочь вычислить её быстрее Добавлено @ 10:37 Цитата(T0ohtik @  28.4.2008,  23:19 ) Для начала надо посчитать коэфициенты С. Сколько будит их? (max(m,n))  ? тут нужно учесть, что БПФ рассчитано на работу с векторами, длина которых - степень двойки (есть, правда, модификации, где вектор делится не на два, а на большее количество подпоследовательностей, но там, по-моему, всё менее оптимально и сложно) так что сначала, думаю, нужно округлить оба вектора до степени двойки вверх (для длинных чисел это делается вполне естественным образом - дополнение нулями)

Автор: maxdiver 29.4.2008, 20:04

Размер C будет равен max(N,M), округленному вверх до степени двойки, и ещё умножить на 2. Если, конечно, реализация алгоритма именно такова, что приводит размеры к степеням двойки (а другая реализация, как уже было сказано, менее популярна).

Автор: ksili 30.4.2008, 05:03

И что будет быстрее, чем просто в столбик?

Автор: Mayk 30.4.2008, 05:22

Цитата(ksili @  30.4.2008,  09:03 ) И что будет быстрее, чем просто в столбик?  Да. На кывте говорится о такой оценке: Цитата(http://www.rsdn.ru/forum/message/573145.flat.aspx#573145) Ищи описание алгоритма Шёнхаге-Штрассена для умножения целых чисел. Там используется БПФ. Его сложность O( n*log(n)*log( log(n) ));

Автор: maxdiver 30.4.2008, 08:45

В Кормене описывается алгоритм умножения с использованием БПФ с асимптотикой N log N (авторы алгоритма там не приводятся, возможно, что алгоритм "народный"  ). Кстати, у меня на сайте есть (http://maximal.hocomua.ru/fft_multiply.htm) описание этого алгоритма, (собсно перепечатанное из Кормена  ) и работающий код. Так что желающие могут потестировать и сравнить скорость. Скажу, что например умножение двух 100-тысячных чисел он выполняет примерно за секунду (чего ни столбиком, ни, я думаю, даже Карацубой не добиться).

Автор: maxim1000 30.4.2008, 10:03

Цитата(ksili @  30.4.2008,  05:03 ) И что будет быстрее, чем просто в столбик? для достаточно длинных чисел (точного или даже приблизительного порога не знаю - надо на практике смотреть) БПФ всё-таки значительно сложнее свёртки в смысле структуры алгоритма и разнообразных вычислений и для маленьких длин, скорее всего, будет неэффективно но вот с ростом длины рост сложности у него меньше

Автор: ksili 30.4.2008, 10:07

Так там смысл в том, что каждый разряд числа считается за элемент вектора, и таким образом два длинных числа превращаются в два многомерных вектора, которые потом надо свернуть?

Автор: maxim1000 30.4.2008, 15:01

да кстати, нашёл один pdf-ник по этой теме: http://edu.krasu.ru/file.php/6/articles/LONG.pdf

Автор: T0ohtik 30.4.2008, 21:58

maxim1000, этот pdf был взят с сайта http://algolist.manual.ru/maths/longnum.php. Собственно говоря, оттуда я и узнал об этом алгоритме, но т.к. в математике я слабовато, то решил создать тему, чтоб задавать вопросы. Хочу сам написать реализацию, хотя в том же pdf она уже и приведена. Кстати о скорости выполнения она равна Log N и судя по графику, приведенному в том же pdf, алгоритм оптимально использовать при длине числа -  от 100 до 120 и от 180 разрядов.  Ну а сейчас спатки, завтра начнем разбираться далее

Автор: T0ohtik 2.5.2008, 11:30

Что-то я не пойму. В чем заключается этот алгоритм. Наши числа хранятся в массиве А и В - т.е. получается, они представлены как многочлен.  Далее мы находим коэффициенты свертки, с. Которые равны сумме произведений всех элементов, индексы которых в сумме равно i. И эти коэффициенты будут произведением числа А на В. Правильно? Или я в чем - то ошибся.

Автор: maxim1000 2.5.2008, 13:26

да, только там есть ещё один нюанс: после свёртки некоторые коэффициенты могут оказаться больше максимальной цифры например 5*100+12*10+3 для десятичной так что нужно будет пробежаться и поперекидывать лишние части в старшие разряды: (5+1)*100+2*10+3

Автор: T0ohtik 2.5.2008, 14:29

Т.е. другими словами алгоритм умножения десятичных чиел заключается в следующем: 1. выбираем максимальную длину из 2х векторов А и В. Количество коэффициентов свертки С у нас будит равно длине максимального вектора. 2. Дополняем меньший из векторов 0 до размеров большего из векторов. 3. Производим свертку векторов.  4. После свертки проверяем вектор и в случаи переполнения делаем переносы. Это и все? Не ясен момент с округлением  векторов до степени двойки вверх. Зачем? Обязательно?

Автор: maxim1000 2.5.2008, 19:17

Цитата(T0ohtik @  2.5.2008,  14:29 ) Это и все? Не ясен момент с округлением  векторов до степени двойки вверх. Зачем? Обязательно? БПФ, как его обычно строят, состоит в том, чтобы: 1. разбить последовательность на две подпоследовательности одинаковой длины 2. посчитать спектры для них 3. посчитать спектр полной последовательности из двух спектров подпоследовательностей одинаковость длин подпоследовательностей означает, что длина исходной последовательности должна быть чётной, на каждом уровне длины отличаются в два раза т.к. на самом низком уровне у нас последовательности из одного элемента, то длина исходной должна быть степенью двойки Добавлено @ 19:25 в принципе, если на каком-то шаге у нас получилась последовательность нечётной длины, то можно посчитать спектр обычным способом - суммой произведений т.е. если у нас длина 2*2*2*2*3, может, и можно обойтись без дополнения до степени двойки но вот если 2*2*12345, удар по производительности может быть большой а реализовывать отличение первых случаев от вторых, ИМХО, не будет стоить выгод

Автор: T0ohtik 3.5.2008, 15:26

maxim1000,  что-то я совсем запутался... Давай по порядку. Возьмем вектора A - 5 элементов и вектор В - 6 элементов. Цитата(maxim1000 @  2.5.2008,  19:17 ) 1. разбить последовательность на две подпоследовательности одинаковой длины В нашем случае разбивать то не надо, у нас же уже вектора и так 2. Цитата(maxim1000 @  2.5.2008,  19:17 ) 2. посчитать спектры для них Что это значит - по-простому найти коэффициент свертки? Если да, то зачем нам надо 3 шаг.  maxim1000, если есть ссылки, дай их, пожалуйста.

Автор: maxim1000 3.5.2008, 23:41

это я всё говорил про БПФ - ведь его нужно применять к каждому вектору, потом перемножить результаты поэлементно, потом - обратное БПФ так вот предыдущий текст был объяснением, почему длины векторов, которые подаются на вход БПФ должны быть степенями двойки (ну или, по крайней мере, почему в этом случае всё значительно проще) посмотрел wikipedia (БПФ), там в конце есть ссылка: http://alglib.sources.ru/fft/ на первый взгляд интересная там, кстати, в числе прочего обсуждается, что наиболее распространённый алгоритм - для векторов длины степени двойки, но есть и куча других

Автор: T0ohtik 5.5.2008, 23:24

Спасибо! Поглядим

Автор: DRUID3 31.5.2008, 18:10

Цитата(maxim1000 @  3.5.2008,  23:41 ) это я всё говорил про БПФ - ведь его нужно применять к каждому вектору, потом перемножить результаты поэлементно, потом - обратное БПФ так вот предыдущий текст был объяснением, почему длины векторов, которые подаются на вход БПФ должны быть степенями двойки (ну или, по крайней мере, почему в этом случае всё значительно проще) посмотрел wikipedia (БПФ), там в конце есть ссылка: http://alglib.sources.ru/fft/ на первый взгляд интересная Цитата(maxim1000 @  3.5.2008,  23:41 ) там, кстати, в числе прочего обсуждается, что наиболее распространённый алгоритм - для векторов длины степени двойки, но есть и куча других   Есть обобщенные методы - которые позволяют организовать вычисление БПФ для практически любого основания. Но их "быстрость" серьезно уступает RADIX-алгоритмам (2, 4, 8, 16 ... etc.) RADIX использует периодичность и симметрию матрицы преобразования - которая появляется в следствие использования гармонических функций как базиса. Простыми словами в ДПФ одни и те же отсчеты функции много раз (вот тут то и избыток) перемножаются с одними и теми же значениями базисных функций - а RADIX этот избыток убирает. Но... Вот тут самое интересное. Зачем автору топика гармонические функции??? Корреляция с ними будет вычислена с погрешностью и для float и для int! Метод даст толко приблизительное решение! FFT "проканывает" в радиотехнике и других измерениях где погрешность методики можно свести к величине меньшего порядка нежели погрешность самого измерения (АЦП) - но для численных методов математики это вызывает большие затруднения! Не грамотнее ли будет выбрать другой базис - так называемые теоретико-числовые преобразования. Это специально разработанный базис из дискретных целочисленных(!!!) функций - дает 0-ю погрешность при вычислениях свертки для int!!! Во всем остальном он очень похож на Фурье - подходят теже алгоритмы (RADIX например) но без комплексной арифметики (в ДПФ/БПФ это обусловленно функцией основания exp(w)=cos(w)+j*sin(w)). Если свертка циклическая (т.е. за пределами окра рассмотрения ничего нет - как раз наш случай) то тоже дополняется "0"-ми при несовпадении длин последовательностей... Исходники (готовые) есть в сети и распространяются вместе с очень хорошей книгой fxtbook... Кстати, если погуглить, то можно найти недавнее обсуждение данной методики на форуме электронщиков.

Автор: maxdiver 31.5.2008, 21:25

DRUID3 Во, спасибо, меня как раз интересует модификация ДПФ для целочисленной арифметики. Что-то никак не получалось сделать, чтобы работало. Попробую поискать эту книгу fxtbook...

Автор: DRUID3 31.5.2008, 21:43

Цитата(maxdiver @  31.5.2008,  21:25 ) DRUID3Во, спасибо, меня как раз интересует модификация ДПФ для целочисленной арифметики. Что-то никак не получалось сделать, чтобы работало. Попробую поискать эту книгу fxtbook... Ну не надо путать одно с другим. Теоретико-числовые преобразования - это свертка с "0"-й погрешностью. А челочисленное FFT это целочисленное FFT ))). В БПФ Вы, наверное, просто не делаете масштабирования результатов после каждого прохода бабочек - а для intFFT это принципиально - иначе на любом мало-мальски длинном промежутке вычислений набежит просто ужаснейшая погрешность из-за переполнения типа данных int. Быстрое переполнение типа или задействованного регистра на операции умножение с накоплением (а это основа ЦОС в современном его представлении) это "бич" всех целочисленных вычислений и вычислителей - от микроконтроллеров до ПЛИС - для того в DSP и применяют аккумуляторную архитектуру с 48 битами. Вобщем масштабируйте и будет Вам счастье... Вот книга http://www.jjj.de/fxt/ + уникальнейшая библиотека ЦОС исходников! Не первый год тестируемая и в своем составе с тестовыми алгоритмами, если захотите изменить исходник... Вот FFT на все случаи жизни http://www.jjj.de/fft/fftpage.html Покопайтесь, там есть и целочисленные и очень редкие типа RADIX5...

Автор: maxdiver 31.5.2008, 22:14

DRUID3 Не, у меня вообще чушь получалась. даже для размеров 2x2 ) А вообще, раз тут (в intFFT) тоже возникают проблемы, может, лучше уж писать обычное FFT и не париться... ) В любом случае, спасибо за ссылки.

Автор: DRUID3 31.5.2008, 23:35

обычное в смысле где? если там где требуется свертка целочисленных отсчетов функций - то FFT вообще применять не стоит. А автору топика нужна именно она. А если где-то лежит выбор между intFFT и floatFFT - то следует выбрать плавающую точку если это позволяет процессор. А если не позволяет - ARM7, BlackFin, AVR32 etc. то intFFT вполне реализуемый алгоритм. И так просто от него отказываться не стОит, написав его самостоятельно получите уникальные навыки - я сам начал заниматься ЦОС именно с написания своего FFT

Автор: ksili 2.6.2008, 07:34

Ссылка хорошая... наверно )) Скачал fxt-2008.04.10.tgz, там внутри файл с расширением 10. Что с ним делать-то? Скачал PDF в архиве fxtbook.pdf.gz, Acrobat Reader говорит что файл повреждён...

Автор: DRUID3 2.6.2008, 09:12

Цитата(ksili @  2.6.2008,  07:34 ) Ссылка хорошая... наверно ))Скачал fxt-2008.04.10.tgz, там внутри файл с расширением 10. Что с ним делать-то?Скачал PDF в архиве fxtbook.pdf.gz, Acrobat Reader говорит что файл повреждён...  linuxовый архив поврежден, а там есть вообще без архива pdf-ка ее и качать...  Файл с расширением .10 это такая шутка от юниксоидов. Но 7-zip его элементарно открывает, исходники внутри этого файла...

Автор: ksili 2.6.2008, 10:19

Точно. Усё получилось

Powered by Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)