【笔记】Big Integer - Schönhage–Strassen Multiplication

Schönhage–Strassen Multiplication

Integer FFT

/// <summary>
/// Cooley–Tukey FFT algorithm,
/// </summary>
public static Complex[] DFTReal(int[] data)
{
    var result = data.Select(i => new Complex(i, 0d)).ToArray();

    DFT(result);

    return result;
}

/// <summary>
/// Cooley–Tukey FFT algorithm,
/// </summary>
public static int[] IDFTReal(Complex[] data)
{
    IDFT(data);

    return data.Select(c => (int)Math.Round(c.Real)).ToArray();
}

FFT算法参考《Cooley–Tukey FFT Algorithm - Iterative Edition》

将大数转为长度k的向量，k = CeilingPowerOf2(2max(len(a), len(b)) - 1)。其中前半部分为大数反向阿拉伯数字，余下补0。

然后计算两个向量的卷积，卷积结果的进位计算即为最终结果。优点就是整个计算都不涉及64位以上大数计算。

/// <summary>
/// Schönhage-Strassen multiplication
/// </summary>
public static BigInteger SchönhageStrassen(BigInteger left, BigInteger right)
{
    var leftDigits = BigInteger.Abs(left).ToString().ToCharArray().Select(ch => ch - '0').ToArray();
    var rightDigits = BigInteger.Abs(right).ToString().ToCharArray().Select(ch => ch - '0').ToArray();
    var n = Math.Max(leftDigits.Length, rightDigits.Length);
    // Padding to power of two >= 2N–1 so we can apply FFT.
    var m = (int)CeilingPowerOf2(((uint)n << 1) - 1);
    Array.Reverse(leftDigits);
    Array.Reverse(rightDigits);
    Array.Resize(ref leftDigits, m);
    Array.Resize(ref rightDigits, m);

    var leftDFT = DFTReal(leftDigits);
    var rightDFT = DFTReal(rightDigits);
    for(var i = 0; i < m; ++i)
    {
        leftDFT[i] *= rightDFT[i];
    }

    var idft = IDFTReal(leftDFT);
    idft = idft.Reverse().SkipWhile(num => num == 0).ToArray();
    for(var i = idft.Length - 1; i > 0; --i)
    {
        idft[i - 1] += idft[i] / 10;
        idft[i] %= 10;
    }
    var digits = idft.Select(num => num.ToString());
    var result = BigInteger.Parse(string.Join("", digits));
    if (left.Sign * right.Sign < 0)
        result = -result;

    return result;
}

/// <summary>
/// least power of 2 greater than or equal to <paramref name="n"/>
/// </summary>
public static uint CeilingPowerOf2(uint n)
{
    --n;
    n |= (n >> 1);
    n |= (n >> 2);
    n |= (n >> 4);
    n |= (n >> 8);
    n |= (n >> 16);

    return n + 1;
}