Big integer Karatsuba algorithm

[oscbyspro]

I'm adding Karatsuba multiplication to UIntXL (#33).

---

Computing the 10,000,000th Fibonacci element now takes 2.6 seconds on my M1 MacBook Pro. With long multiplication only, the same calculation takes 34.4 seconds. The threshold for using the Karatsuba algorithm is 40 words. Edit: I have reduced it to 2.3 seconds and 20 words, with some improvements.

/// https://www.wolframalpha.com/input?i2d=true&i=fibonacci+10000000
///
/// - Note: The 10,000,000th element contains 2,089,877 decimal digits.
///
func test10000000() {
    NBK.blackHole(NBKFibonacciXL(NBK.blackHoleIdentity(10_000_000)))
}

[oscbyspro]

GNU says the threshold can be reduced to 10 words, in case I feel like optimizing it.

[oscbyspro]

Hm. I think I could do a bit better if I had composable (perhaps generic) methods like:

public static func increment(
_ base: inout UnsafeMutableBufferPointer<UInt>,
by rhs: UnsafeBufferPointer<UInt>, 
times lhs: UnsafeBufferPointer<UInt>) { 
    if  lhs.count < 40 || rhs.count < 40 { 
        ...
    }   else {
        ... 
    }
}

public static func incrementByLongAlgorithm(
_ base: inout UnsafeMutableBufferPointer<UInt>,
by rhs: UnsafeBufferPointer<UInt>, 
times lhs: UnsafeBufferPointer<UInt>) { 
    ...
}

public static func incrementByKaratsubaAlgorithm(
_ base: inout UnsafeMutableBufferPointer<UInt>,
by rhs: UnsafeBufferPointer<UInt>, 
times lhs: UnsafeBufferPointer<UInt>) { 
    ...
}

By that I mean, I wrote something quite neat and it'll take a day to clean up and push.

---

Edit: initialize(_:to:times:) seems more appropriate for multiplication algorithms.

[oscbyspro]

I'll be doing fewer allocations, so it appears that I can lower the threshold to 20 words.

[oscbyspro]

I can remove an additional allocation and some additions if I make a decrementing version too. Hm. Dunno whether it's worth it, but it'd be neat insofar as every level of recursion allocating 2(lhs.count + rhs.count) contiguous memory.

[oscbyspro]

Welp. It's possible to add the middle product in-place, but when I tried it I got random results. So I suppose there might be some finicky life-time thing to think about when doing it. Or maybe the recursion breaks exclusive access. It's not an issue I want to tackle at the moment, however.

[oscbyspro]

I want to eventually end up with a non-recursive algorithm, but I'll towards it incrementally.

[oscbyspro]

Hm. The middle product can reuse the 1st product's allocation. It's at least something 🤷‍♂️

square: 3 * a.count, non-square: 2 * (a.count + b.count)

[oscbyspro]

I chose to normalize each subsequence because I saw this in 10_000_000_000 ^ 1_000:

[silly string] - [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]

[oscbyspro]

I suppose there is still some work to do. Zzz.

• 10_000_000_000 ^ 1_000 works
• 10 ^ 10_000 triggers a debug assertion

[oscbyspro]

It's late. I'll fix it tomorrow. Also, I think I can write z0 snd z1 directly to the product and then allocate z2.

[oscbyspro]

Observation. It works if I floor the partition, not ceil it. Is that a coincidence? Feels like it. Maybe not? Zzz.

[oscbyspro]

Hm. I rather confident that flooring solves it, because then the high part is larger and it must be trimmable to fit with a half shift because otherwise it would not fit with a full shift, which it must because that's just z1. I already have an assertion for it. floor(_:) good, ceil(_:) bad, on floor can sleep.

[oscbyspro]

Test Case '-[NBKFlexibleWidthKitBenchmarks.NBKFibonacciXLBenchmarks testNoLoop10000000]' passed (2.290 seconds).