также использует биномиальные коэффициенты. Для их расчета можно использовать формулу, выражающую биномиальный коэффициент через факториалы: 
или использовать рекуррентную формулу:
Из бинома Ньютона и рекуррентной формулы ясно, что биномиальные коэффициенты — целые числа. Одним из методов, позволяющих значительно сократить количество вычислений, является применение Фурье преобразований и дискретных Фурье преобразований.
Наличие большого числа библиотек, реализующих Фурье преобразований (во всевозможных вариантах быстрых версий), делает реализацию алгоритмов не очень сложной задачей для программирования.
Реализованные алгоритмы являются частью библиотеки с открытым исходным кодом FFTTools. Интернет-адрес: github.com/dprotopopov/FFTTools
Преобразование Фурье функции f вещественной переменной является интегральным и задаётся следующей формулой:
Разные источники могут давать определения, отличающиеся от приведённого выше выбором коэффициента перед интегралом, а также знака «?» в показателе экспоненты. Но все свойства будут те же, хотя вид некоторых формул может измениться.
Кроме того, существуют разнообразные обобщения данного понятия.
Дискретное преобразование Фурье
Дискретное преобразование Фурье (в англоязычной литературе DFT, Discrete Fourier Transform) — это одно из преобразований Фурье, широко применяемых в алгоритмах цифровой обработки сигналов (его модификации применяются в сжатии звука в MP3, сжатии изображений в JPEG и др.), а также в других областях, связанных с анализом частот в дискретном (к примеру, оцифрованном аналоговом) сигнале. Дискретное преобразование Фурье требует в качестве входа дискретную функцию. Такие функции часто создаются путём дискретизации (выборки значений из непрерывных функций). Дискретные преобразования Фурье помогают решать дифференциальные уравнения в частных производных и выполнять такие операции, как свёртки. Дискретные преобразования Фурье также активно используются в статистике, при анализе временных рядов. Существуют многомерные дискретные преобразования Фурье.
Формулы дискретных преобразований
Прямое преобразование:
Обратное преобразование:
Дискретное преобразование Фурье является линейным преобразованием, которое переводит вектор временных отсчётов в вектор спектральных отсчётов той же длины. Таким образом преобразование может быть реализовано как умножение симметричной квадратной матрицы на вектор:
Свёртка двух функций
Согласно определению, свёрткой двух функций F и G называется функция FхG:
FхG(t) = SUM F(i)*G(j)|i+j=t
Фурье-преобразования для вычисления свёртки
Одним из замечательных свойств преобразований Фурье является возможность быстрого вычисления корреляции двух функций определённых, либо на действительном аргументе (при использовании классической формулы), либо на конечном кольце (при использовании дискретных преобразований).
И хотя подобные свойства присущи многим линейным преобразованиям, для практического применения, для вычисления операции свёртки, согласно данному нами определению, используется формула
FхG = BFT ( FFT(F)*FFT(G) )
Где
- FFT – операция прямого преобразования Фурье
- BFT – операция обратного преобразования Фурье
Проверить правильность равенства довольно легко – явно подставив в формулы Фурье-преобразований и сократив получившиеся формулы
Биномиальные коэффициенты
Рассмотрим полином F(x)=1+x и его свёртку с самим собой n раз
Fx..xF = SUM С( i, n-1 )*x^i = BFT ( FFT(F)*...*FFT(F) ) = BFT ( FFT(F)^(n-1) )
То есть биномиальные коэффициенты С( i, n-1 ) могут быть получены из значений коэффициентов полинома (1+x)^(n-1)
Программируем:
using System;
using System.Drawing;
using System.Linq;
using System.Numerics;
namespace FFTTools
{
public class BinomialBuilder : BuilderBase, IBuilder
{
/// <summary>
/// Performs application-defined tasks associated with freeing, releasing, or resetting unmanaged resources.
/// </summary>
public void Dispose()
{
}
public static void GetLongs(long[] array, long x = 1)
{
var n = array.Length - 1;
if (array.Length > 0) array[0] = 1;
for (var i = 0; i < n; i++)
array[i + 1] = x*array[i]*(n - i)/(i + 1);
}
public static void GetDoubles(double[] array, double x = 1.0)
{
var complex = new Complex[array.Length];
if (array.Length > 0) complex[0] = Complex.One;
if (array.Length > 1) complex[1] = x;
if (array.Length > 0)
{
Fourier(complex, FourierDirection.Forward);
complex = complex.Select(
value => Complex.Pow(value, array.Length - 1)/array.Length).ToArray();
Fourier(complex, FourierDirection.Backward);
}
var index = 0;
foreach (var value in complex) array[index++] = value.Real;
}
public Bitmap ToBitmap(Bitmap source)
{
throw new NotImplementedException();
}
}
}
Проверяем:
using System;
using System.Linq;
using Microsoft.VisualStudio.TestTools.UnitTesting;
namespace FFTTools.UnitTest
{
[TestClass]
public class BinomialUnitTest
{
[TestMethod]
public void BinomialTestMethod()
{
const int count = 10;
var doubles = new double[count];
var longs = new long[count];
BinomialBuilder.GetLongs(longs);
BinomialBuilder.GetDoubles(doubles);
Console.WriteLine(
string.Join(Environment.NewLine,
longs.Zip(doubles, (x, y) => string.Format("{0} - {1} = {2}", x, y, x - y))) +
Environment.NewLine);
Assert.IsTrue(doubles.Zip(longs, (x, y) => x - y).All(x => Math.Abs(x) < 0.001));
}
}
}
Зачем?
При вычислении с помощью треугольника Паскаля трудоёмкость имеет оценку O(n^2).
При вычислении с помощью быстрых Фурье-преобразований трудоёмкость имеет оценку O(n*log n).
Примечание:
В статье заимствованы фрагменты из статьи Расчет биномиальных коэффициентов на Си (С++)
P.S.
Трудоёмкость расчета биномиальных коэффициентов может быть уменьшена до O(n):
Cn[0]=1
Cn[i+1] = Cn[i]*(n-i)/(i+1)
Доказательство:
Cn[i]*(n-i)/(i+1) = n!/((n-i)!i!) * (n-i)/(i+1) = n!/((n-i-1)!(i+1)!) = Cn[i+1]
Комментарии (30)
koldyr
09.01.2016 14:45меня терзают жуткие сомнения. можно посмотреть расчет погрешностей?
dprotopopov
09.01.2016 14:49расчет погрешностей не делал.
при расчёте надо смотреть на упаковку double в соответствии с IEEE — это мантиса+экспонета = 64 бита
brainick
09.01.2016 15:58-3>>Для их расчета можно использовать формулу, выражающую биномиальный коэффициент через факториалы: <<
really?
chersanya
09.01.2016 18:51Было бы значительно лучше, если бы вы добавили реальную скорость этого алгоритма и какого-нибудь из упомянутого выше поста. Интересно, начиная с каких N получается выигрыш?
dprotopopov
09.01.2016 19:44-4реальную скорость алгоритма
Вот вам и тема для исследования.
Хотя всё это можно оценить и теоретически, зная спецификацию процессора.
В статье есть небольшая обманка — возведение в степень тоже может быть выполнено за O(log n) операций умножения на дискретном процессоре.
Bodigrim
09.01.2016 21:20> При вычислении с помощью треугольника Паскаля трудоёмкость имеет оценку O(n^2).
> При вычислении с помощью быстрых Фурье-преобразований трудоёмкость имеет оценку O(n*log n).
Сам массив значений биномиальных коэффициентов от С_n^0 до C_n^n занимает что-то типа O(n^2) бит. Как вы умудрились их посчитать за асимптотически меньшее время?chersanya
09.01.2016 21:24+1Конечно здесь имеется в виду приближённое вычисление, без длинной арифметики — вон, в коде double повсюду. Преобразование фурье с длинными числами вообще редкая в применении вещь.
Bodigrim
09.01.2016 21:29+2Тогда с тем же успехом можно сложить первые эннадцать рядом треугольника Паскаля в массив и объявить о революционном алгоритме, который вычисляет биномиальные коэффициенты за O(1).
Для приближенного вычисления больших коэффициентов уж лучше воспользоваться рядом Стирлинга.chersanya
09.01.2016 22:12Можно конечно, в double по прикидкам влезает около 500 рядов треугольника — то есть массив размером в мегабайт. Приведённый автором поста алгоритм показывает другой подход, который возможно в чём-то лучше — хотя тоже согласен про ряд Стирлинга.
dprotopopov
09.01.2016 22:32при x=0.5 влезет больше (см. параметры метода GetDoubles)
только полученные коэффициенты надо будет умножать на 2^i
halyavin
10.01.2016 01:25+3Чем вас не устроил стандартный алгоритм за O(n)?
dprotopopov
10.01.2016 07:51Согласен, что для приближенного вычисления больших коэффициентов лучше воспользоваться рядом Стирлинга.
Просто мне интересны сейчас темы с использованием Фурье-преобразований.halyavin
10.01.2016 09:15Ряд Стирлинга асимптотический и на контесте вы с ним Accepted никогда не получите. Особенно, если биномиальные коэффициенты нужно по простому модулю считать. Я говорю про стандартный алгоритм, который все пишут за 2-5 минут.
dprotopopov
10.01.2016 09:41все пишут за 2-5 минут — это треугольник Паскаля.
При вычислении с помощью треугольника Паскаля трудоёмкость имеет оценку O(n^2).
или покажите формулу/алгоритм — можно в виде публикации.
… я не волшебник — я только учусь…halyavin
10.01.2016 10:43+2Теперь, когда вы знаете, что он существует и тривиален, вам не составит труда его придумать. Разве можно лишать вас этого удовольствия?
dprotopopov
10.01.2016 10:50умножить-разделить? — это имеется ввиду?
Просто мне интересны сейчас темы с использованием Фурье-преобразований.
dprotopopov
10.01.2016 21:14внёс изменения в статью и в библиотеку FFTTools.
спасибо.
public static void GetLongs(long[] array, long x = 1)
{
var n = array.Length — 1;
if (array.Length > 0) array[0] = 1;
for (var i = 0; i < n; i++)
array[i + 1] = x*array[i]*(n — i)/(i + 1);
}
masai
10.01.2016 18:09Учитывая интересы автора, наверное, тем, что в нём нет преобразования Фурье. :)
Sayonji
10.01.2016 07:13Рассмотрим полином F(x)=1+x и его свёртку с самим собой n раз
Fx..xF = SUM С( i, n-1 )*x^i
Я не понял это место.
По каким значениям t надо суммировать при свёртке, чтобы получить то что написано?
dprotopopov
10.01.2016 07:42Fx..xF(t) = SUM F(i1(0))*...*F(i(n-1)) | i1(0)+...+i(n-1)=t
dprotopopov
10.01.2016 08:18описался — вот правильно:
Fx..xF(t) = SUM F(i(1))*...*F(i(n-1)) | i(1)+...+i(n-1)=t
x=1
F(0) == 1
F(1) == xSayonji
10.01.2016 17:43По-прежнему непонятно.
F(x)=1+x
F(0) == 1
F(1) == x
Вот в этом нет смысла, к примеру.dprotopopov
10.01.2016 18:04Да, сам я виноват — использовал обозначения у полинома и у вектора одинаковые.
полином от переменной x может быть записан как вектор коэффициентов полинома/многочлена.
чтобы не плодить обозначений через F обозначаю и полином и вектор коэффициентов. Прочтение зависит от контекста.
F(x)=1+x — это полином/многочлен
x=1 — это присвоение переменной x значения
F[0] = 1 — это элемент с индексом 0 в векторе
F[1] = x — это элемент с индексом 1 в векторе
Fx..xF[t] = SUM F[i[1]]*...*F[i[n-1]] | i[1]+...+i[n-1]=t
dcc0
10.01.2016 09:14«в расчете биномиальныХ коэффициентов.»
У Вас опечатка.
Вставлю немного своих дилетантских соображений:
Интересно, все ли обращали внимание на, казалось бы, элементарную вещь с выводом треугольника Паскаля (если через массивы), что считать можно только половину строки: ), а массив можно просто переворачивать, правда, нужна проверка на четность, чтобы переворачивать и такие строки +6+15+20+15+6+.dprotopopov
10.01.2016 09:33Спасибо, поправил текст.
Да, конечно все понимают, что биномиальные коэффициенты симметричны, причём с ростом n график приближается к графику нормального распределения (после соответствующего нормирования).
причина по которой все смотрят не на половину — формула (a+b)^n — при подстановке реальных чисел — там вычисленные значения разложения различаются — то есть вычисляют значения C(i,n)*a^i*b^(n-i)
nata16k8
10.01.2016 20:44+1Замечания по обозначениям.
1. В преобразовании Фурье для аналоговых сигналов используется коэффициент 1/sqrt(2pi), а в ДПФ — 1/2pi. Логично использовать что-то одно.
2. А с крышечкой — это симметричная квадратная матрица, а не вектор (не точно скопирован текст из википедии)
3. Обратное БПФ — это обычно ifft. Как расшифровать BFT? В яндекс-поиске первое после ссылок на Ваши тексты — BFT: Bit Filtration Technique.dprotopopov
11.01.2016 10:08Где
FFT – операция прямого преобразования Фурье
BFT – операция обратного преобразования Фурье
zelserg
Интересно. В статье прямо заимствованы фрагменты из моей статьи Расчет биномиальных коэффициентов на Си (С++) , но об этом даже не упомянуто.
Я полагаю, что расчет C(67,33) по треугольнику Паскаля с запоминанием (дальше они не влезают в unsigned64) потребует меньше времени, чем этот Фурье несмотря на «худший» O (n-то не велико, а алгоритм расчета намного сложнее). А уж расчет всех коэффициентов для n<=67 по треугольнику вообще самый быстрый.
PapaBubaDiop
Для n=67 было бы интересно сравнить время Паскаль vs Фурье. Можно в тактовых единицах.
dprotopopov
Добавил ссылку на статью habrahabr.ru/post/274689