10

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?

5
  • 1
    Could you give some context on what Karatsuba is? Commented Feb 2, 2010 at 19:49
  • 2
    I'm not sure if this will help but maybe you can somehow cast it to a BitArray so that you can count the bits.
    – AaronLS
    Commented Feb 2, 2010 at 19:51
  • @aaronls: That is a lot faster, thanks.
    – PythonPower
    Commented Feb 2, 2010 at 21:05
  • 2
    @Chris: en.wikipedia.org/wiki/Karatsuba_algorithm Commented Feb 6, 2010 at 20:48
  • << has a lower precedence than +/- Commented Feb 20, 2011 at 5:03

2 Answers 2

13

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.

4
  • 1
    Magnificent work Alexander Higgins! +1 for your answer which helped me in my search for perfect numbers...
    – RvdV79
    Commented May 8, 2013 at 13:05
  • Fascinating, but from a brief microbenchmark it seems .Net already uses this optimization; timings are close, with this sometimes being a bit faster, but on average (without doing the math) the default implementation seems to win by a narrow margin.
    – Alex Paven
    Commented Aug 29, 2015 at 16:09
  • 1
    In practice there is a cutoff after which "grade school" loses to "Karatsuba" due to the 4/3 improvement. This is due to the algorithm overhead cburch.com/proj/karat/results.html. There are also other fast multiplication methods, such as Toom-Cook which is 9/5 improvement over "grade school" but in practice Karatsuba will win out up to some cutoff due to overhead. Likewise, there is a cutoff in the range of several thousands of digits where Fast Fourier beats out Toom-Cook. Commented Feb 20, 2019 at 21:54
  • @AlexanderHiggins : i checked out that URL you listed, and I think their methodology unrealistically penalizes grade-school base case. Since their cutoff for Karatsuba is 4 digits, then at least grade-school oughta to be benchmarked exactly at 4 digits to make it at least semi meaningful. Doing it 1-digit by 1-digit then even someone physically multiplying by their fingers can beat it. Commented Apr 6 at 2:30
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.