Calculation of binomial coefficients using Fourier transforms

When solving combinatorics problems, it is often necessary to calculate binomial coefficients. Newton binomial i.e. expansion also uses binomial coefficients. To calculate them, you can use the formula expressing the binomial coefficient in terms of factorials: or use the recurrence formula: It is clear from Newton’s binomial and recurrence formula that the binomial coefficients are integers.

One of the methods to significantly reduce the number of calculations is the use of Fourier transforms and discrete Fourier transforms.

The presence of a large number of libraries that implement Fourier transforms (in all possible versions of the fast versions) makes the implementation of algorithms not a very difficult task for programming.
The implemented algorithms are part of the open source library FFTTools. Internet address: github.com/dprotopopov/FFTTools


The Fourier transform of a function f of a real variable is integral and is given by the following formula:



Different sources can give definitions that differ from the one given above by choosing the coefficient before the integral, as well as the “-” sign in the exponent. But all the properties will be the same, although the appearance of some formulas may change.

In addition, there are various generalizations of this concept.

Discrete Fourier Transform

Discrete Fourier transform (in the English literature DFT, Discrete Fourier Transform) is one of the Fourier transforms widely used in digital signal processing algorithms (its modifications are used in audio compression to MP3, image compression in JPEG, etc.), as well as in others areas related to the analysis of frequencies in a discrete (for example, digitized analog) signal. The discrete Fourier transform requires a discrete function as an input. Such functions are often created by discretization (sampling values ​​from continuous functions). Discrete Fourier transforms help solve partial differential equations and perform operations such as convolutions. Discrete Fourier transforms are also actively used in statistics, in the analysis of time series. There are multidimensional discrete Fourier transforms.

Discrete Transform Formulas

Direct Transform:



Inverse Transform: The



Discrete Fourier Transform is a linear transform that translates a vector of time samples into a vector of spectral samples of the same length. Thus, the transformation can be implemented as multiplying a symmetric square matrix by a vector:

Convolution of two functions

By definition, the convolution of two functions F and G is the function FxG:

FхG (t) = SUM F (i) * G (j) | i + j = t

Fourier transforms for convolution calculation

One of the remarkable properties of Fourier transforms is the ability to quickly calculate the correlation of two functions defined, either on a real argument (using the classical formula) or on a finite ring (when using discrete transforms).

And although similar properties are inherent in many linear transformations, for practical use, to calculate the convolution operation, according to our definition, the formula is used

FxG = BFT (FFT (F) * FFT (G))

Where

• FFT - direct Fourier transform operation
• BFT - Inverse Fourier Transform Operation

It is quite easy to verify the correctness of the equality by explicitly substituting the Fourier transform formulas and reducing the resulting formulas

Binomial Coefficients

Consider the polynomial F (x) = 1 + x and its convolution with itself n times
Fx..xF = SUM С (i, n-1) * x ^ i = BFT (FFT (F) * ... * FFT ( F)) = BFT (FFT (F) ^ (n-1))
That is, the binomial coefficients C (i, n-1) can be obtained from the values ​​of the coefficients of the polynomial (1 + x) ^ (n-1)

Programming:

using System;
using System.Drawing;
using System.Linq;
using System.Numerics;
namespace FFTTools
{
    public class BinomialBuilder : BuilderBase, IBuilder
    {
        /// 

        ///     Performs application-defined tasks associated with freeing, releasing, or resetting unmanaged resources.
        /// 
        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();
        }
    }
}

We check:

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));
        }
    }
}

What for?

When calculating using the Pascal triangle, the complexity is estimated as O (n ^ 2).
When calculating using fast Fourier transforms, the complexity is estimated as O (n * log n).

Note:

The article borrowed fragments from the article Calculation of binomial coefficients in C (C ++)

PS The
complexity of calculating binomial coefficients can be reduced to O (n):
Cn [0] = 1
Cn [i + 1] = Cn [i] * (ni ) / (i + 1)
Proof:
Cn [i] * (ni) / (i + 1) = n! / ((ni)! i!) * (ni) / (i + 1) = n! / ( (ni-1)! (i + 1)!) = Cn [i + 1]