Tuesday, July 01, 2008

Anatolii Karatsuba meets Leonardo de Pisa

In the posts regarding Fibonacci, it came up that the VisualWorks' large integer multiplication primitive isn't the most suitable to calculate products of truly huge numbers. In this case, the term "truly huge" refers to numbers which can easily have one million bits, so they are way beyond what would be considered normal usage. And nevertheless...

I looked around and found Karatsuba's algorithm for multiplication. While it would be nice to put it in the VM, it is also possible to implement it in the image by letting it use the existing primitive when the number sizes are suitable for it.

After a bit of massaging the implementation, which is not the most efficient one because again this should be in the VM or otherwise in C, I got Karatsuba to tie the VM's performance at about 2k bits. Here are some measurements after that.

  • 3k bits: Karatsuba 12% faster than the VM.
  • 4k bits: Karatsuba 6% faster than the VM.
  • 8k bits: Karatsuba 33% faster than the VM.
  • 16k bits: Karatsuba 31% faster than the VM.
  • 64k bits: Karatsuba 2.51x faster than the VM.
  • 256k bits: Karatsuba 3.96x faster than the VM.
  • 1m bits: Karatsuba 8.85x faster than the VM.
So I subclassed the F(a+b) Fibonacci calculator so that it uses Karatsuba instead of regular multiplication. Result?
  • 1953152 fibonacci: Karatsuba F(a+b) 4.38x faster than the regular F(a+b) calculator.
How about that...

    SmallInteger>>karatsubaTimes: anInteger

      ^anInteger * self


    LargeInteger>>karatsubaTimes: aLargeInteger

      | splitThreshold
        selfHigh selfLow aLargeIntegerHigh aLargeIntegerLow
        highProduct lowProduct mixedProduct lowTerm |
      self digitLength + aLargeInteger digitLength > 768
        ifFalse: [^self * aLargeInteger].
      splitThreshold := self digitLength
        + aLargeInteger digitLength
        * 2.
      selfHigh := self bitShift: splitThreshold negated.
      selfLow := self - (selfHigh bitShift: splitThreshold).
      aLargeIntegerHigh := aLargeInteger
        bitShift: splitThreshold negated.
      aLargeIntegerLow := aLargeInteger
        - (aLargeIntegerHigh bitShift: splitThreshold).
      highProduct := selfHigh karatsubaTimes: aLargeIntegerHigh.
      lowProduct := selfLow karatsubaTimes: aLargeIntegerLow.
      mixedProduct := selfLow + selfHigh karatsubaTimes:
        aLargeIntegerLow + aLargeIntegerHigh.
      lowTerm := mixedProduct - (lowProduct + highProduct).
      ^((highProduct bitShift: splitThreshold)
        + lowTerm bitShift: splitThreshold)
          + lowProduct

4 comments:

nicolas cellier said...

I always wanted to implement, but the idea to modify the (Squeak) VM was boring...
Did not even think i could do it in Smalltalk!

I have a question however, given LargeInteger implementation, wouldn't a byteShift be more efficient than a bitShift?

Andres said...

Nicolas,

IIRC, the bitShift primitives in VW check the lower 3 bits and if they are zero they just push the bytes around. I do not know about other Smalltalks, but it seems to me it wouldn't be difficult to do.

Andres.

Andres said...

Nicolas,

Also, keep in mind that back in the late 90s, Squeak implemented LargeIntegers without primitives.

Andres.

simple said...

Surely a bitshift?