Optimizing Karatsuba Implementation

Question

So, I'm trying to improve some of the operations that .net 4's BigInteger class provide since the operations appear to be quadratic. I've made a rough Karatsuba implementation but it's still slower than I'd expect.

The main problem seems to be that BigInteger provides no simple way to count the number of bits and, so, I have to use BigInteger.Log(..., 2). According to Visual Studio, about 80-90% of the time is spent calculating logarithms.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Numerics;

namespace Test
{
    class Program
    {
        static BigInteger Karatsuba(BigInteger x, BigInteger y)
        {
            int n = (int)Math.Max(BigInteger.Log(x, 2), BigInteger.Log(y, 2));
            if (n <= 10000) return x * y;

            n = ((n+1) / 2);

            BigInteger b = x >> n;
            BigInteger a = x - (b << n);
            BigInteger d = y >> n;
            BigInteger c = y - (d << n);

            BigInteger ac = Karatsuba(a, c);
            BigInteger bd = Karatsuba(b, d);
            BigInteger abcd = Karatsuba(a+b, c+d);

            return ac + ((abcd - ac - bd) << n) + (bd << (2 * n));
        }

        static void Main(string[] args)
        {
            BigInteger x = BigInteger.One << 500000 - 1;
            BigInteger y = BigInteger.One << 600000 + 1;
            BigInteger z = 0, q;

            Console.WriteLine("Working...");
            DateTime t;

            // Test standard multiplication
            t = DateTime.Now;
            z = x * y;
            Console.WriteLine(DateTime.Now - t);

            // Test Karatsuba multiplication
            t = DateTime.Now;
            q = Karatsuba(x, y);
            Console.WriteLine(DateTime.Now - t);

            // Check they're equal
            Console.WriteLine(z == q);

            Console.Read();
        }
    }
}

So, what can I do to speed it up?

Answer 1

commenting on @Alexander Higgins's response :

ac + ((abcd - ac - bd) << n) + (bd << (2 * n))

This is still doing duplicative work when could use nesting layers :

ac + ((abcd - ac - bd + (bd << n)) << n)

Resulting in 1 less pair of ( ) and 1 less multiplication.

But personally I don't really think Karatsuba saves all that much, since at the bigint level, subtractions are slightly more costly than additions (and you're doing 2 of them at each recursion level),

and don't forget that to even create abcd, it requires YET another extra pair of bigint additions.

Factor in all that extra processing it's almost too much hassle compared to just going with the naive approach :

(using AH`s same letters for consistency) :

    BigInteger b = x >> n;
    BigInteger a = x - (b << n);
    BigInteger d = y >> n;
    BigInteger c = y - (d << n);

---

return ac + ((ad + bc + (bd << n)) << n)

which doesn't look all that bad when compared to

return ac + ((abcd - ac - bd) << n) + (bd << (2 * n))

with the extra benefit of not requiring storing any partial products for the subtraction parts.

I usually pre-determine the sign before any recursion begins, perform all ops in unsigned realm, and simply do 1 single final negate if necessary.

Answer 2

Why count all of the bits?

In vb I do this:

<Runtime.CompilerServices.Extension()> _
Function BitLength(ByVal n As BigInteger) As Integer
    Dim Data() As Byte = n.ToByteArray
    Dim result As Integer = (Data.Length - 1) * 8
    Dim Msb As Byte = Data(Data.Length - 1)
    While Msb
        result += 1
        Msb >>= 1
    End While
    Return result
End Function

In C# it would be:

public static int BitLength(this BigInteger n)
{
    byte[] Data = n.ToByteArray();
    int result = (Data.Length - 1) * 8;
    byte Msb = Data[Data.Length - 1];
    while (Msb != 0) {
        result += 1;
        Msb >>= 1;
    }
    return result;
}

Finally...

    static BigInteger Karatsuba(BigInteger x, BigInteger y)
    {
        int n = (int)Math.Max(x.BitLength(), y.BitLength());
        if (n <= 10000) return x * y;

        n = ((n+1) / 2);

        BigInteger b = x >> n;
        BigInteger a = x - (b << n);
        BigInteger d = y >> n;
        BigInteger c = y - (d << n);

        BigInteger ac = Karatsuba(a, c);
        BigInteger bd = Karatsuba(b, d);
        BigInteger abcd = Karatsuba(a+b, c+d);

        return ac + ((abcd - ac - bd) << n) + (bd << (2 * n));
    }

Calling the extension method may slow things down so perhaps this would be faster:

int n = (int)Math.Max(BitLength(x), BitLength(y));

FYI: with the bit length method you can also calculate a good approximation of the log much faster than the BigInteger Method.

bits = BitLength(a) - 1;
log_a = (double)i * log(2.0);

As far as accessing the internal UInt32 Array of the BigInteger structure, here is a hack for that.

import the reflection namespace

Private Shared ArrM As MethodInfo
Private Shard Bits As FieldInfo
Shared Sub New()
    ArrM = GetType(System.Numerics.BigInteger).GetMethod("ToUInt32Array", BindingFlags.NonPublic Or BindingFlags.Instance)
    Bits = GetType(System.Numerics.BigInteger).GetMember("_bits", BindingFlags.NonPublic Or BindingFlags.Instance)(0)

End Sub
<Extension()> _
Public Function ToUInt32Array(ByVal Value As System.Numerics.BigInteger) As UInteger()
    Dim Result() As UInteger = ArrM.Invoke(Value, Nothing)
    If Result(Result.Length - 1) = 0 Then
        ReDim Preserve Result(Result.Length - 2)
    End If
    Return Result
End Function

Then you can get the underlying UInteger() of the big integer as

 Dim Data() As UInteger = ToUInt32Array(Value)
 Length = Data.Length 

or Alternately

Dim Data() As UInteger = Value.ToUInt32Array()

Note that _bits fieldinfo can be used to directly access the underlying UInteger() _bits field of the BigInteger structure. This is faster than invoking the ToUInt32Array() method. However, when BigInteger B <= UInteger.MaxValue _bits is nothing. I suspect that as an optimization when a BigInteger fits the size of a 32 bit (machine size) word MS returns performs normal machine word arithmetic using the native data type.

I have also not been able to use the _bits.SetValue(B, Data()) as you normally would be able to using reflection. To work around this I use the BigInteger(bytes() b) constructor which has overhead. In c# you can use unsafe pointer operations to cast a UInteger() to Byte(). Since there are no pointer ops in VB, I use Buffer.BlockCopy. When access the data this way it is important to note that if the MSB of the bytes() array is set, MS interprets it as a Negative number. I would prefer they made a constructor with a separate sign field. The word array is to add an addition 0 byte to make uncheck the MSB

Also, when squaring you can improve even further

 Function KaratsubaSquare(ByVal x As BigInteger)
    Dim n As Integer = BitLength(x) 'Math.Max(BitLength(x), BitLength(y))

    If (n <= KaraCutoff) Then Return x * x
    n = ((n + 1) >> 1)

    Dim b As BigInteger = x >> n
    Dim a As BigInteger = x - (b << n)
    Dim ac As BigInteger = KaratsubaSquare(a)
    Dim bd As BigInteger = KaratsubaSquare(b)
    Dim c As BigInteger = Karatsuba(a, b)
    Return ac + (c << (n + 1)) + (bd << (2 * n))

End Function

This eliminates 2 shifts, 2 additions and 3 subtractions from each recursion of your multiplication algorithm.