Error in Karatsuba's algorithm for Long Arithmetic

I'm trying to make a library with long numbers. And there was a problem that I had been trying to fix for a week. The problem is that I stupidly do not understand where the problem is in my code. And more specifically, the problem is in the implementation of the Karatsuba multiplication algorithm

There is an array of Bits that stores the expansion coefficients of a long number in the number system 2^32.

    private uint[] Bits;
    private const Int64 Base = (long)UInt32.MaxValue + 1;

The multiplication method itself takes arrays of Bits of 2 instances and outputs an array (multiplication result)

public uint[] KaratsubaMultiplication(uint[] x, uint[] y)
    {
        int maxLength = 0;

        if (x.Length > y.Length) maxLength = x.Length;
        else maxLength = y.Length;
        uint[] result = new uint[maxLength * 2];
        if (maxLength == 1)
        {
            ulong mul = (ulong)x[0] * (ulong)y[0];
            if (mul >= Base)
            {
                result[0] = (uint)mul;
                result[1] = (uint)(mul >> 32);
            }
            else result[0] = (uint)mul;
        }
        else
        {
            //a
            uint[] BitsX = ExtendArray(x, maxLength);
            //b
            uint[] BitsY = ExtendArray(y, maxLength);

            //а1 первая часть числа a (тк с конца массив у меня)
            uint[] a1 = TakeElements(BitsX, BitsX.Length / 2, BitsX.Length / 2);

            //a2 вторая часть числа a
            uint[] a2 = TakeElements(BitsX, BitsX.Length / 2);

            //b1 первая часть числа b
            uint[] b1 = TakeElements(BitsY, BitsY.Length / 2, BitsY.Length / 2);

            //b2 вторая часть числа b
            uint[] b2 = TakeElements(BitsY, BitsY.Length / 2);

            //Произведение первых частей чисел
            uint[] a1b1 = KaratsubaMultiplication(a1, b1);

            //Произведение вторых частей чисел
            uint[] a2b2 = KaratsubaMultiplication(a2, b2);

            //Сумма произведений первых и вторых частей
            uint[] a1b1Suma2b2 = Addition(a1b1, a2b2);

            // (a1+a2)(b1+b2) - (a1b1 + a2b2), то что перед - это вот эта переменная хах
            uint[] ab = KaratsubaMultiplication(Addition(a1, a2), Addition(b1, b2));

            ab = Substraction(ab, a1b1Suma2b2);

            result = Addition(Finalize(a1b1, BitsX.Length), Finalize(ab, BitsX.Length / 2));
            result = Addition(result, a2b2);
        }
        result = Normalize(result);
        return result;
    }

Subtraction method:

public static uint[] Substraction(uint[] x, uint[] y)
    {
        if (x.Length > y.Length)
            y = NormalizeArrays(y, x.Length);
        else x = NormalizeArrays(x, y.Length);
        uint[] result = new uint[x.Length];
        ulong carry = 0;
        int i = 0;
        for (i = 0; i < x.Length; i++)
        {
            if (x[i] - carry < y[i])
            {
                result[i] = (uint)((ulong)Base + x[i] + carry - y[i]);
                carry = 1;
            }
            else
            {
                result[i] = (uint)(x[i] + carry - y[i]);
                carry = 0;
            }
        }
        result = Normalize(result);
        return result;
    }

Addition method:

public static uint[] Addition(uint[] x, uint[] y, long radix = Base)
    {
        if (x.Length > y.Length) y = NormalizeArrays(y, x.Length);
        else if (y.Length > x.Length) x = NormalizeArrays(x, y.Length);

        uint[] result = new uint[x.Length + y.Length];
        ulong carry = 0;
        int i = 0;

        if (radix == Base)
        {
            for (; i < x.Length; i++)
            {
                ulong sum = x[i] + y[i] + carry;
                result[i] = (uint)sum;
                carry = sum >> 32;
            }
        }
        else
        {
            for (; i < x.Length; i++)
            {
                ulong sum = x[i] + y[i] + carry;
                result[i] = (uint)(sum % (ulong)radix);
                carry = sum / (ulong)radix;
            }
        }
        if (carry != 0) result[i] = (uint)carry;
        result = Normalize(result);
        return result;
    }

I work in the number system 2^32. And the Bits array stores the expansion coefficients of the long number в обратном порядке

All additional methods for the code to work:

private static uint[] ExtendArray(uint[] array, int length)
    {
        uint[] temp = new uint[length + length % 2];
        Array.Copy(array, temp, array.Length);
        return temp;
    }

    public static uint[] TakeElements(uint[] array, int length, int startIndex = 0)
    {
        uint[] temp = new uint[length];
        Array.Copy(array, startIndex, temp, 0, length);
        return temp;
    }
/// <summary>
    /// Подстраивает массив меньшей длины к массиву с большей
    /// </summary>
    /// <param name="array"></param>
    /// <param name="length"></param>
    /// <returns></returns>
    private static uint[] NormalizeArrays(uint[] array, int length)
    {
        uint[] temp = new uint[length];
        Array.Copy(array, temp, array.Length);
        return temp;
    }
    /// <summary>
    /// Удаляет не лидирующие 0 в массиве
    /// </summary>
    /// <param name="value"></param>
    /// <returns></returns>
    private static uint[] Normalize(uint[] value)
    {
        for (int i = value.Length - 1; i >= 0; i--)
        {
            if (value[i] != 0)
            {
                uint[] temp = new uint[i + 1];
                Array.Copy(value, 0, temp, 0, i + 1);
                return temp;
            }
        }
        return new uint[] { 0 };
    }
    static uint[] Finalize(uint[] res, int n)
    {
        uint[] temp = new uint[res.Length + n];
        for (int i = 0; i < res.Length; i++)
            temp[i + n] = res[i];
        return temp;
    }

I checked the code on the numbers 635115234552377 and 423523523623. Which were parsed from the string.

You can check the code simply by passing 2 arrays to the Karatsuba multiplication method, {1240623673, 147874} and {2616728615, 98} are numbers 635115234552377 and 423523523623

And I get {14491652, 960651648, 835046072, 755855689}

I expect to see 268986242044270826190301871, but I get 59885057380827124403963171489844109316 these numbers are already in the 10 number system.

I have seen a lot of implementations of this method, but none of them was completely clear to me.

EXTENSION:

In the karatsuba multiplication method, at the end, I passed prod1 and prod2 to the Finalize method, and at the end, I added prod1 and not prod2 to the result. (I also changed the code in the question)

This move led me to the right one To my satisfaction, I decided to multiply the number 268986242044270826190301871 by itself. I expected to get 72353598409099050336288860303289482073245694106100641, but I get 72353598409098710053921939364826055592126409756992417. And I still haven't figured out what the mistake is.

EXTENSION:

Based on the article in wikipedia

I changed and corrected the method. (The code in the question is changed)

There were a couple of questions. In Wikipedia, without multiplication, a1b1 is added to the result, that is, the product of the first parts of the numbers.

But the correct answer is I will only get it if I swap a1b1 with a2b2 and back again. But it only works for the example with 635115234552377 and 423523523623.

If I multiply the number 268986242044270826190301871 by itself and calculate the product in 2 ways ( that is, I will first add a1b1 and 2 method a2b2), then in any case I will get the wrong answer.

Author: EOF, 2019-08-23

1 answers

I solved the problem. To begin with, I replaced the subtraction method.

public static uint[] Substraction(uint[] x, uint[] y, long radix = Base)
    {
        if (x.Length > y.Length)
            y = NormalizeArrays(y, x.Length);
        else x = NormalizeArrays(x, y.Length);
        uint[] result = new uint[x.Length];
        long carry = 0;
        int i = 0;
        if (radix == Base)
        {
            for (i = 0; i < x.Length; i++)
            {
                //carry = carry + x[i] - y[i] + 10;
                //result[i] = (uint)carry % 10;
                //if (carry < 10) carry = -1;
                //else carry = 0;
                carry = carry + x[i] - y[i] + radix;
                result[i] = (uint)carry;
                if (carry < radix) carry = -1;
                else carry = 0;
            }
        }
        else
        {
            for(i = 0; i < x.Length; i++)
            {
                carry = carry + x[i] - y[i] + radix;
                result[i] = (uint)(carry % radix);
                if (carry < radix) carry = -1;
                else carry = 0;
            }
        }
        result = Normalize(result);
        return result;
    }

Also, the multiplication method changed and corrected

public uint[] KaratsubaMultiplication(uint[] x, uint[] y, long radix = Base)
    {
        int maxLength = 0;
        if (x.Length > y.Length) maxLength = x.Length;
        else maxLength = y.Length;
        uint[] result = new uint[maxLength * 2];
        if (maxLength == 1)
        {
            ulong mul = (ulong)x[0] * (ulong)y[0];
            if (radix == Base)
            {
                //if(mul > 10)
                if (mul > Base)
                {
                    result[0] = (uint)mul;
                    result[1] = (uint)(mul >> 32);
                    //result[0] = (uint)mul % 10;
                    //result[1] = (uint)mul / 10;
                }
                else result[0] = (uint)mul;
            }
            else
            {
                result[0] = (uint)(mul % (ulong)radix);
                result[1] = (uint)(mul / (ulong)radix);
            }
        }
        else
        {
            //a
            uint[] BitsX = ExtendArray(x, maxLength);
            //b
            uint[] BitsY = ExtendArray(y, maxLength);

            //а1 первая часть числа a (тк с конца массив у меня)
            uint[] a1 = TakeElements(BitsX, BitsX.Length / 2, BitsX.Length / 2);

            //a2 вторая часть числа a
            uint[] a2 = TakeElements(BitsX, BitsX.Length / 2);

            //b1 первая часть числа b
            uint[] b1 = TakeElements(BitsY, BitsY.Length / 2, BitsY.Length / 2);

            //b2 вторая часть числа b
            uint[] b2 = TakeElements(BitsY, BitsY.Length / 2);

            //Произведение первых частей чисел
            uint[] a1b1 = KaratsubaMultiplication(a1, b1, radix);

            //Произведение вторых частей чисел
            uint[] a2b2 = KaratsubaMultiplication(a2, b2, radix);

            //Сумма произведений первых и вторых частей
            uint[] a1b1Suma2b2 = Addition(a1b1, a2b2, radix);

            // (a1+a2)(b1+b2) - (a1b1 + a2b2), то что перед - это вот эта переменная хах
            uint[] ab = KaratsubaMultiplication(Addition(a1, a2, radix), Addition(b1, b2, radix), radix);

            ab = Substraction(ab, a1b1Suma2b2, radix);

            result = Addition(Finalize(a1b1, BitsX.Length), Finalize(ab, BitsX.Length / 2), radix);
            result = Addition(result, a2b2, radix);
        }
        result = Normalize(result);
        return result;
    }

And the most infuriating error in the code. This is for example ulong sum = x[i] + y[i] + carry; tc x[i] and y[i] this is uint, then as far as I understand and the result outputs in uint and then leads to ulong. And if there is an overflow, the sum will not be correct. So I in everywhere in such moments (in addition and subtraction) for example x[i] bring to ulong.

 0
Author: Павел Ериков, 2019-08-24 16:12:18