Бухгалтерский учёт
0,0
рейтинг
9 января в 13:38

Разработка → Расчет биномиальных коэффициентов с использованием Фурье-преобразований

При решении задач комбинаторики часто возникает необходимость в расчете биномиальных коэффициентов. Бином Ньютона, т.е. разложение image также использует биномиальные коэффициенты. Для их расчета можно использовать формулу, выражающую биномиальный коэффициент через факториалы: image или использовать рекуррентную формулу:image Из бинома Ньютона и рекуррентной формулы ясно, что биномиальные коэффициенты — целые числа.

Одним из методов, позволяющих значительно сократить количество вычислений, является применение Фурье преобразований и дискретных Фурье преобразований.

Наличие большого числа библиотек, реализующих Фурье преобразований (во всевозможных вариантах быстрых версий), делает реализацию алгоритмов не очень сложной задачей для программирования.
Реализованные алгоритмы являются частью библиотеки с открытым исходным кодом FFTTools. Интернет-адрес: github.com/dprotopopov/FFTTools


Преобразование Фурье функции f вещественной переменной является интегральным и задаётся следующей формулой:

Преобразование Фурье

Разные источники могут давать определения, отличающиеся от приведённого выше выбором коэффициента перед интегралом, а также знака «−» в показателе экспоненты. Но все свойства будут те же, хотя вид некоторых формул может измениться.

Кроме того, существуют разнообразные обобщения данного понятия.

Дискретное преобразование Фурье


Дискретное преобразование Фурье (в англоязычной литературе DFT, Discrete Fourier Transform) — это одно из преобразований Фурье, широко применяемых в алгоритмах цифровой обработки сигналов (его модификации применяются в сжатии звука в MP3, сжатии изображений в JPEG и др.), а также в других областях, связанных с анализом частот в дискретном (к примеру, оцифрованном аналоговом) сигнале. Дискретное преобразование Фурье требует в качестве входа дискретную функцию. Такие функции часто создаются путём дискретизации (выборки значений из непрерывных функций). Дискретные преобразования Фурье помогают решать дифференциальные уравнения в частных производных и выполнять такие операции, как свёртки. Дискретные преобразования Фурье также активно используются в статистике, при анализе временных рядов. Существуют многомерные дискретные преобразования Фурье.

Формулы дискретных преобразований


Прямое преобразование:

image

Обратное преобразование:

image

Дискретное преобразование Фурье является линейным преобразованием, которое переводит вектор временных отсчётов в вектор спектральных отсчётов той же длины. Таким образом преобразование может быть реализовано как умножение симметричной квадратной матрицы на вектор:

image

Свёртка двух функций


Согласно определению, свёрткой двух функций 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]
Дмитрий Протопопов @dprotopopov
карма
10,0
рейтинг 0,0
Бухгалтерский учёт
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Спецпроект

Самое читаемое Разработка

Комментарии (31)

  • +2
    Интересно. В статье прямо заимствованы фрагменты из моей статьи Расчет биномиальных коэффициентов на Си (С++) , но об этом даже не упомянуто.
    Я полагаю, что расчет C(67,33) по треугольнику Паскаля с запоминанием (дальше они не влезают в unsigned64) потребует меньше времени, чем этот Фурье несмотря на «худший» O (n-то не велико, а алгоритм расчета намного сложнее). А уж расчет всех коэффициентов для n<=67 по треугольнику вообще самый быстрый.
    • 0
      Для n=67 было бы интересно сравнить время Паскаль vs Фурье. Можно в тактовых единицах.
    • 0
      Добавил ссылку на статью habrahabr.ru/post/274689
  • 0
    меня терзают жуткие сомнения. можно посмотреть расчет погрешностей?
    • 0
      расчет погрешностей не делал.
      при расчёте надо смотреть на упаковку double в соответствии с IEEE — это мантиса+экспонета = 64 бита
  • НЛО прилетело и опубликовало эту надпись здесь
  • 0
    Было бы значительно лучше, если бы вы добавили реальную скорость этого алгоритма и какого-нибудь из упомянутого выше поста. Интересно, начиная с каких N получается выигрыш?
    • –4
      реальную скорость алгоритма

      Вот вам и тема для исследования.
      Хотя всё это можно оценить и теоретически, зная спецификацию процессора.
      В статье есть небольшая обманка — возведение в степень тоже может быть выполнено за O(log n) операций умножения на дискретном процессоре.
  • 0
    > При вычислении с помощью треугольника Паскаля трудоёмкость имеет оценку O(n^2).
    > При вычислении с помощью быстрых Фурье-преобразований трудоёмкость имеет оценку O(n*log n).

    Сам массив значений биномиальных коэффициентов от С_n^0 до C_n^n занимает что-то типа O(n^2) бит. Как вы умудрились их посчитать за асимптотически меньшее время?
    • +1
      Конечно здесь имеется в виду приближённое вычисление, без длинной арифметики — вон, в коде double повсюду. Преобразование фурье с длинными числами вообще редкая в применении вещь.
      • +2
        Тогда с тем же успехом можно сложить первые эннадцать рядом треугольника Паскаля в массив и объявить о революционном алгоритме, который вычисляет биномиальные коэффициенты за O(1).
        Для приближенного вычисления больших коэффициентов уж лучше воспользоваться рядом Стирлинга.
        • 0
          Можно конечно, в double по прикидкам влезает около 500 рядов треугольника — то есть массив размером в мегабайт. Приведённый автором поста алгоритм показывает другой подход, который возможно в чём-то лучше — хотя тоже согласен про ряд Стирлинга.
          • 0
            при x=0.5 влезет больше (см. параметры метода GetDoubles)
            только полученные коэффициенты надо будет умножать на 2^i
  • +3
    Чем вас не устроил стандартный алгоритм за O(n)?
    • 0
      Согласен, что для приближенного вычисления больших коэффициентов лучше воспользоваться рядом Стирлинга.
      Просто мне интересны сейчас темы с использованием Фурье-преобразований.
      • 0
        Ряд Стирлинга асимптотический и на контесте вы с ним Accepted никогда не получите. Особенно, если биномиальные коэффициенты нужно по простому модулю считать. Я говорю про стандартный алгоритм, который все пишут за 2-5 минут.
        • 0
          все пишут за 2-5 минут — это треугольник Паскаля.
          При вычислении с помощью треугольника Паскаля трудоёмкость имеет оценку O(n^2).
          или покажите формулу/алгоритм — можно в виде публикации.
          … я не волшебник — я только учусь…
          • +2
            Теперь, когда вы знаете, что он существует и тривиален, вам не составит труда его придумать. Разве можно лишать вас этого удовольствия?
            • 0
              умножить-разделить? — это имеется ввиду?
              Просто мне интересны сейчас темы с использованием Фурье-преобразований.
            • 0
              внёс изменения в статью и в библиотеку 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);
              }
    • 0
      Учитывая интересы автора, наверное, тем, что в нём нет преобразования Фурье. :)
  • 0
    Рассмотрим полином F(x)=1+x и его свёртку с самим собой n раз
    Fx..xF = SUM С( i, n-1 )*x^i

    Я не понял это место.
    По каким значениям t надо суммировать при свёртке, чтобы получить то что написано?
  • 0
    Fx..xF(t) = SUM F(i1(0))*...*F(i(n-1)) | i1(0)+...+i(n-1)=t
    • 0
      описался — вот правильно:
      Fx..xF(t) = SUM F(i(1))*...*F(i(n-1)) | i(1)+...+i(n-1)=t
      x=1
      F(0) == 1
      F(1) == x
      • 0
        По-прежнему непонятно.
        F(x)=1+x
        F(0) == 1
        F(1) == x

        Вот в этом нет смысла, к примеру.
        • 0
          Да, сам я виноват — использовал обозначения у полинома и у вектора одинаковые.
          полином от переменной 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
  • 0
    «в расчете биномиальныХ коэффициентов.»
    У Вас опечатка.
    Вставлю немного своих дилетантских соображений:
    Интересно, все ли обращали внимание на, казалось бы, элементарную вещь с выводом треугольника Паскаля (если через массивы), что считать можно только половину строки: ), а массив можно просто переворачивать, правда, нужна проверка на четность, чтобы переворачивать и такие строки +6+15+20+15+6+.
    • 0
      Спасибо, поправил текст.
      Да, конечно все понимают, что биномиальные коэффициенты симметричны, причём с ростом n график приближается к графику нормального распределения (после соответствующего нормирования).
      причина по которой все смотрят не на половину — формула (a+b)^n — при подстановке реальных чисел — там вычисленные значения разложения различаются — то есть вычисляют значения C(i,n)*a^i*b^(n-i)
  • +1
    Замечания по обозначениям.
    1. В преобразовании Фурье для аналоговых сигналов используется коэффициент 1/sqrt(2pi), а в ДПФ — 1/2pi. Логично использовать что-то одно.
    2. А с крышечкой — это симметричная квадратная матрица, а не вектор (не точно скопирован текст из википедии)
    3. Обратное БПФ — это обычно ifft. Как расшифровать BFT? В яндекс-поиске первое после ссылок на Ваши тексты — BFT: Bit Filtration Technique.
    • 0
      Где
      FFT – операция прямого преобразования Фурье
      BFT – операция обратного преобразования Фурье
      • 0

        FFT = Fast Fourier Transform
        BFT = ?


        И если есть общепринятая нотация, не стоит изобретать своё.

Только зарегистрированные пользователи могут оставлять комментарии. Войдите, пожалуйста.