From dd62d0e94fa1d46d7568597ad7f2d35374d2b514 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Oscar=20Bystr=C3=B6m=20Ericsson?= Date: Tue, 7 Nov 2023 13:19:21 +0100 Subject: [PATCH] Big integer multiplication algorithms (#108). --- .../{NBK+Allocation.swift => NBK+Data.swift} | 20 +- ...gnedInteger+Multiplication+Karatsuba.swift | 226 ++++++++++++++++++ ...tUnsignedInteger+Multiplication+Long.swift | 112 +++++++++ ...StrictUnsignedInteger+Multiplication.swift | 47 ++++ .../NBKStrictUnsignedInteger+Partition.swift | 31 +++ .../{NBK+Allocation.swift => NBK+Data.swift} | 4 +- .../NBKStrictUnsignedInteger+Partition.swift | 84 +++++++ 7 files changed, 521 insertions(+), 3 deletions(-) rename Sources/NBKCoreKit/Private/{NBK+Allocation.swift => NBK+Data.swift} (75%) create mode 100644 Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication+Karatsuba.swift create mode 100644 Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication+Long.swift create mode 100644 Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Partition.swift rename Tests/NBKCoreKitTests/Private/{NBK+Allocation.swift => NBK+Data.swift} (95%) create mode 100644 Tests/NBKCoreKitTests/Private/NBKStrictUnsignedInteger+Partition.swift diff --git a/Sources/NBKCoreKit/Private/NBK+Allocation.swift b/Sources/NBKCoreKit/Private/NBK+Data.swift similarity index 75% rename from Sources/NBKCoreKit/Private/NBK+Allocation.swift rename to Sources/NBKCoreKit/Private/NBK+Data.swift index ed42af35..8b195feb 100644 --- a/Sources/NBKCoreKit/Private/NBK+Allocation.swift +++ b/Sources/NBKCoreKit/Private/NBK+Data.swift @@ -8,11 +8,29 @@ //=----------------------------------------------------------------------------= //*============================================================================* -// MARK: * NBK x Allocation +// MARK: * NBK x Data //*============================================================================* extension NBK { + //=------------------------------------------------------------------------= + // MARK: Initializers + //=------------------------------------------------------------------------= + + /// Initializes `base` to `source` then returns the initialized `count`. + @inlinable public static func initializeGetCount( + _ base: UnsafeMutablePointer, to source: UnsafeBufferPointer) -> Int { + base.initialize(from: source.baseAddress!, count: source.count) + return source.count as Int + } + + /// Initializes `base` to `source` then returns the initialized `count`. + @inlinable public static func initializeGetCount( + _ base: UnsafeMutablePointer, repeating element: T, count: Int) -> Int { + base.initialize(repeating: element, count: count) + return count as Int + } + //=------------------------------------------------------------------------= // MARK: Utilities //=------------------------------------------------------------------------= diff --git a/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication+Karatsuba.swift b/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication+Karatsuba.swift new file mode 100644 index 00000000..341aa436 --- /dev/null +++ b/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication+Karatsuba.swift @@ -0,0 +1,226 @@ +//=----------------------------------------------------------------------------= +// This source file is part of the Numberick open source project. +// +// Copyright (c) 2023 Oscar Byström Ericsson +// Licensed under Apache License, Version 2.0 +// +// See http://www.apache.org/licenses/LICENSE-2.0 for license information. +//=----------------------------------------------------------------------------= + +//*============================================================================* +// MARK: * NBK x Strict Unsigned Integer x Mul. x Karatsuba x Sub Sequence +//*============================================================================* +//=----------------------------------------------------------------------------= +// MARK: + where Base is Unsafe Buffer Pointer +//=----------------------------------------------------------------------------= + +extension NBK.StrictUnsignedInteger.SubSequence where Base: MutableCollection { + + //=------------------------------------------------------------------------= + // MARK: Initializers + //=------------------------------------------------------------------------= + + /// Initializes `base` to the [Karatsuba][algorithm] product of `lhs` and `rhs`. + /// + /// + /// high : z1 low : z0 + /// ┌───────────┴───────────┬───────────┴───────────┐ + /// │ x1 * y1 │ x0 * y0 │ + /// └───────────┬───────────┴───────────┬───────────┘ + /// add │ x0 * y0 │ + /// ├───────────────────────┤ + /// add │ x1 * y1 │ + /// ├───────────────────────┤ + /// sub │ (x1 - x0) * (y1 - y0) │ : z2 + /// └───────────────────────┘ + /// + /// + /// - Parameter base: A buffer of size `lhs.count` + `rhs.count`. + /// + /// - Note: The `base` must be uninitialized or `pointee` must be trivial. + /// + /// [algorithm]: https://en.wikipedia.org/wiki/karatsuba_algorithm + /// + @inline(never) @inlinable public static func initializeByKaratsubaAlgorithm( + _ base: inout Base, to lhs: UnsafeBufferPointer, times rhs: UnsafeBufferPointer) + where Base == UnsafeMutableBufferPointer { + //=--------------------------------------= + Swift.assert(base.count == lhs.count + rhs.count, NBK.callsiteOutOfBoundsInfo()) + //=--------------------------------------= + let (k0c): Int = Swift.max(lhs.count, rhs.count) + let (k1c): Int = k0c &>> 1 + let (k2c): Int = k1c &<< 1 + + let (x1s, x0s) = NBK.SUISS.partitionTrimmingRedundantZeros(lhs, at: k1c) + let (y1s, y0s) = NBK.SUISS.partitionTrimmingRedundantZeros(rhs, at: k1c) + let (z1c, z0c) = (x1s.count + y1s.count, x0s.count + y0s.count) + + let count = 2 * (z0c + z1c) as Int + Swift.withUnsafeTemporaryAllocation(of: Base.Element.self, capacity: count) { buffer in + var pointer: UnsafeMutablePointer + //=----------------------------------= + // pointee: deferred deinitialization + //=----------------------------------= + defer { + buffer.baseAddress!.deinitialize(count: count) + } + //=----------------------------------= + // temporary memory regions + //=----------------------------------= + pointer = buffer.baseAddress! as UnsafeMutablePointer + var x0 = UnsafeMutableBufferPointer(start: pointer, count: x0s.count); pointer += x0.count + var x1 = UnsafeMutableBufferPointer(start: pointer, count: x1s.count); pointer += x1.count + var y0 = UnsafeMutableBufferPointer(start: pointer, count: y0s.count); pointer += y0.count + var y1 = UnsafeMutableBufferPointer(start: pointer, count: y1s.count); pointer += y1.count + var z0 = UnsafeMutableBufferPointer(start: pointer, count: z0c); pointer += z0.count + var z1 = UnsafeMutableBufferPointer(start: pointer, count: z1c); pointer += z1.count + Swift.assert(buffer.baseAddress!.distance(to: pointer) == count) + //=----------------------------------= + // pointee: initialization 1 + //=----------------------------------= + x0.baseAddress!.initialize(from: x0s.baseAddress!, count: x0s.count) + x1.baseAddress!.initialize(from: x1s.baseAddress!, count: x1s.count) + y0.baseAddress!.initialize(from: y0s.baseAddress!, count: y0s.count) + y1.baseAddress!.initialize(from: y1s.baseAddress!, count: y1s.count) + + NBK.SUISS.initialize(&z0, to: UnsafeBufferPointer(x0), times: UnsafeBufferPointer(y0)) + NBK.SUISS.initialize(&z1, to: UnsafeBufferPointer(x1), times: UnsafeBufferPointer(y1)) + + let xs = NBK.SUISS.compare(x1, to: x0).isLessThanZero + if xs { + Swift.swap(&x0, &x1) + } + + let ys = NBK.SUISS.compare(y1, to: y0).isLessThanZero + if ys { + Swift.swap(&y0, &y1) + } + + NBK.SUISS.decrement( &x1, by: UnsafeBufferPointer(x0)) + NBK.SUISS.decrement( &y1, by: UnsafeBufferPointer(y0)) + //=----------------------------------= + // product must fit in combined width + //=----------------------------------= + Swift.assert(z1.dropFirst(base[k2c...].count).allSatisfy({ $0.isZero })) + z1 = UnsafeMutableBufferPointer(rebasing: z1.prefix(base[k2c...].count)) + //=----------------------------------= + // pointee: initialization 2 + //=----------------------------------= + pointer = base.baseAddress! as UnsafeMutablePointer + pointer += NBK.initializeGetCount(pointer, to: UnsafeBufferPointer(z0)) + pointer += NBK.initializeGetCount(pointer, repeating: Base.Element.zero, count: k2c - z0.count) + pointer += NBK.initializeGetCount(pointer, to: UnsafeBufferPointer(z1)) + pointer += NBK.initializeGetCount(pointer, repeating: Base.Element.zero, count: base.count - (k2c + z1.count)) + Swift.assert(base.baseAddress!.distance(to: pointer) == base.count) + //=----------------------------------= + var slice = UnsafeMutableBufferPointer(rebasing: base.suffix(from: k1c)) + NBK.SUISS.increment(&slice, by: UnsafeBufferPointer(z0)) + NBK.SUISS.increment(&slice, by: UnsafeBufferPointer(z1)) + + z2: do { // reuse z0 and z1 to form z2. pointee is trivial + z0 = UnsafeMutableBufferPointer(start: z0.baseAddress!, count: x1.count + y1.count) + NBK.SUISS.initialize(&(z0), to: UnsafeBufferPointer(x1), times: UnsafeBufferPointer(y1)) + }; if xs == ys { + NBK.SUISS.decrement(&slice, by: UnsafeBufferPointer(z0)) + } else { + NBK.SUISS.increment(&slice, by: UnsafeBufferPointer(z0)) + } + } + } + + /// Initializes `base` to the square [Karatsuba][algorithm] product of `elements`. + /// + /// + /// high : z1 low : z0 + /// ┌───────────┴───────────┬───────────┴───────────┐ + /// │ x1 * x1 │ x0 * x0 │ + /// └───────────┬───────────┴───────────┬───────────┘ + /// add │ x0 * x0 │ + /// ├───────────────────────┤ + /// add │ x1 * x1 │ + /// ├───────────────────────┤ + /// sub │ (x1 - x0) * (x1 - x0) │ : z2 + /// └───────────────────────┘ + /// + /// + /// - Parameter base: A buffer of size `2 * elements.count`. + /// + /// - Note: The `base` must be uninitialized or `pointee` must be trivial. + /// + /// [algorithm]: https://en.wikipedia.org/wiki/karatsuba_algorithm + /// + /// ### Development + /// + /// - TODO: Add a test case where `elements` has redundant zeros. + /// + @inline(never) @inlinable public static func initializeByKaratsubaAlgorithm( + _ base: inout Base, toSquareProductOf elements: UnsafeBufferPointer) where Base == UnsafeMutableBufferPointer { + //=--------------------------------------= + Swift.assert(base.count == 2 * elements.count, NBK.callsiteOutOfBoundsInfo()) + //=--------------------------------------= + let (k0c): Int = elements.count + let (k1c): Int = k0c &>> 1 + let (k2c): Int = k1c &<< 1 + + let (x1s, x0s) = NBK.SUISS.partitionTrimmingRedundantZeros(elements, at: k1c) + let (z1c, z0c) = (2 * x1s.count, 2 * x0s.count) + + let count = 3 * (x1s.count + x0s.count) as Int + Swift.withUnsafeTemporaryAllocation(of: Base.Element.self, capacity: count) { buffer in + var pointer: UnsafeMutablePointer + //=----------------------------------= + // pointee: deferred deinitialization + //=----------------------------------= + defer { + buffer.baseAddress!.deinitialize(count: count) + } + //=----------------------------------= + // temporary memory regions + //=----------------------------------= + pointer = buffer.baseAddress! as UnsafeMutablePointer + var x0 = UnsafeMutableBufferPointer(start: pointer, count: x0s.count); pointer += x0.count + var x1 = UnsafeMutableBufferPointer(start: pointer, count: x1s.count); pointer += x1.count + var z0 = UnsafeMutableBufferPointer(start: pointer, count: z0c); pointer += z0.count + var z1 = UnsafeMutableBufferPointer(start: pointer, count: z1c); pointer += z1.count + Swift.assert(buffer.baseAddress!.distance(to: pointer) == count) + //=----------------------------------= + // pointee: initialization 1 + //=----------------------------------= + x0.baseAddress!.initialize(from: x0s.baseAddress!, count: x0s.count) + x1.baseAddress!.initialize(from: x1s.baseAddress!, count: x1s.count) + + NBK.SUISS.initialize(&z0, toSquareProductOf: UnsafeBufferPointer(x0)) + NBK.SUISS.initialize(&z1, toSquareProductOf: UnsafeBufferPointer(x1)) + + if NBK.SUISS.compare(x1, to: x0).isLessThanZero { + Swift.swap(&x0, &x1) + } + + NBK.SUISS.decrement( &x1, by: UnsafeBufferPointer(x0)) + //=----------------------------------= + // product must fit in combined width + //=----------------------------------= + Swift.assert(z1.dropFirst(base[k2c...].count).allSatisfy({ $0.isZero })) + z1 = UnsafeMutableBufferPointer(rebasing: z1.prefix(base[k2c...].count)) + //=----------------------------------= + // pointee: initialization 2 + //=----------------------------------= + pointer = base.baseAddress! as UnsafeMutablePointer + pointer += NBK.initializeGetCount(pointer, to: UnsafeBufferPointer(z0)) + pointer += NBK.initializeGetCount(pointer, repeating: Base.Element.zero, count: k2c - z0.count) + pointer += NBK.initializeGetCount(pointer, to: UnsafeBufferPointer(z1)) + pointer += NBK.initializeGetCount(pointer, repeating: Base.Element.zero, count: base.count - (k2c + z1.count)) + Swift.assert(base.baseAddress!.distance(to: pointer) == base.count) + //=----------------------------------= + var slice = UnsafeMutableBufferPointer(rebasing: base.suffix(from: k1c)) + NBK.SUISS.increment(&slice, by: UnsafeBufferPointer(z0)) + NBK.SUISS.increment(&slice, by: UnsafeBufferPointer(z1)) + + z2: do { // reuse z0 and z1 to form z2. pointee is trivial + z0 = UnsafeMutableBufferPointer(start: z0.baseAddress!, count: x1.count &<< 1) + NBK.SUISS.initialize(&(z0), toSquareProductOf: UnsafeBufferPointer(x1)) + NBK.SUISS.decrement(&slice, by:/*-----------*/ UnsafeBufferPointer(z0)) + } + } + } +} diff --git a/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication+Long.swift b/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication+Long.swift new file mode 100644 index 00000000..4657e444 --- /dev/null +++ b/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication+Long.swift @@ -0,0 +1,112 @@ +//=----------------------------------------------------------------------------= +// This source file is part of the Numberick open source project. +// +// Copyright (c) 2023 Oscar Byström Ericsson +// Licensed under Apache License, Version 2.0 +// +// See http://www.apache.org/licenses/LICENSE-2.0 for license information. +//=----------------------------------------------------------------------------= + +//*============================================================================* +// MARK: * NBK x Strict Unsigned Integer x Mul. x Long x Sub Sequence +//*============================================================================* +//=----------------------------------------------------------------------------= +// MARK: + where Base is Unsafe Buffer Pointer +//=----------------------------------------------------------------------------= + +extension NBK.StrictUnsignedInteger.SubSequence where Base: MutableCollection { + + //=------------------------------------------------------------------------= + // MARK: Initializers + //=------------------------------------------------------------------------= + + /// Initializes `base` to the [long][algorithm] product of `lhs` and `rhs` plus `addend`. + /// + /// - Parameter base: A buffer of size `lhs.count` + `rhs.count`. + /// + /// - Note: The `base` must be uninitialized or `pointee` must be trivial. + /// + /// [algorithm]: https://en.wikipedia.org/wiki/multiplication_algorithm + /// + @inline(never) @inlinable public static func initializeByLongAlgorithm( + _ base: inout Base, to lhs: UnsafeBufferPointer, times rhs: UnsafeBufferPointer, + plus addend: Base.Element = .zero) where Base == UnsafeMutableBufferPointer { + //=--------------------------------------= + Swift.assert(base.count == lhs.count + rhs.count, NBK.callsiteOutOfBoundsInfo()) + //=--------------------------------------= + var pointer = base.baseAddress! + //=--------------------------------------= + // pointee: initialization 1 + //=--------------------------------------= + var carry: Base.Element = addend + let first: Base.Element = rhs.first ?? Base.Element.zero + + for element in lhs { + var wide = element.multipliedFullWidth(by: first) + carry = Base.Element(bit: wide.low.addReportingOverflow(carry)) &+ wide.high + pointer.initialize(to: wide.low) // done, uninitialized or discarded pointee + pointer = pointer.successor() + } + + if !rhs.isEmpty { + pointer.initialize(to: carry) + pointer = pointer.successor() + } + + Swift.assert(base.baseAddress!.distance(to: pointer) == lhs.count + Int(bit: !rhs.isEmpty)) + //=--------------------------------------= + // pointee: initialization 2 + //=--------------------------------------= + for var index in rhs.indices.dropFirst() { + pointer.initialize(to: 00000) + pointer = pointer.successor() + NBK.SUISS.incrementInIntersection(&base, by: lhs, times: rhs[index], plus: Base.Element.zero, at: &index) + } + + Swift.assert(base.baseAddress!.distance(to: pointer) == base.count) + } + + /// Initializes `base` to the square [long][algorithm] product of `elements` plus `addend`. + /// + /// - Parameter base: A buffer of size `2 * elements`. + /// + /// - Note: The `base` must be uninitialized or `pointee` must be trivial. + /// + /// [algorithm]: https://en.wikipedia.org/wiki/multiplication_algorithm + /// + @inline(never) @inlinable public static func initializeByLongAlgorithm( + _ base: inout Base, toSquareProductOf elements: UnsafeBufferPointer, plus addend: Base.Element = .zero) + where Base == UnsafeMutableBufferPointer { + //=--------------------------------------= + Swift.assert(base.count == 2 * elements.count, NBK.callsiteOutOfBoundsInfo()) + //=--------------------------------------= + // pointee: initialization + //=--------------------------------------= + base.initialize(repeating: 0 as Base.Element) + //=--------------------------------------= + var index: Int + var carry: Base.Element = addend + //=--------------------------------------= + var baseIndex = elements.startIndex; while baseIndex < elements.endIndex { + let multiplier = elements[ baseIndex] + let productIndex = 2 * baseIndex + elements.formIndex(after: &baseIndex) + + index = productIndex + 1 // add non-diagonal products + + NBK.SUISS.incrementInIntersection( + &base, by: UnsafeBufferPointer(rebasing: elements[baseIndex...]), + times: multiplier, plus: 00000, at: &index) + + index = productIndex // partially double non-diagonal products + + NBK.SUISS.multiply(&base, by: 0002, add: &carry, from: &index, to: productIndex + 2) + + index = productIndex // add this iteration's diagonal product + + carry &+= Base.Element(bit: NBK.SUISS.incrementInIntersection( + &base, by: CollectionOfOne((multiplier)), + times: multiplier, plus: 00000, at: &index)) + } + } +} diff --git a/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication.swift b/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication.swift index e5212400..18eae628 100644 --- a/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication.swift +++ b/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Multiplication.swift @@ -58,3 +58,50 @@ extension NBK.StrictUnsignedInteger.SubSequence where Base: MutableCollection { } } } + +//=----------------------------------------------------------------------------= +// MARK: + where Base is Unsafe Buffer Pointer +//=----------------------------------------------------------------------------= + +extension NBK.StrictUnsignedInteger.SubSequence where Base: MutableCollection { + + //=------------------------------------------------------------------------= + // MARK: Initializers + //=------------------------------------------------------------------------= + + /// Initializes `base` to the product of `lhs` and `rhs`. + /// + /// - Parameters: + /// - base: A buffer of size `lhs.count` + `rhs.count`. + /// + /// - Note: The `base` memory must be uninitialized or the pointee must be a trivial type. + /// + @inlinable public static func initialize( + _ base: inout Base, to lhs: UnsafeBufferPointer, times rhs: UnsafeBufferPointer) + where Base == UnsafeMutableBufferPointer { + //=--------------------------------------= + if lhs.count < 20 || rhs.count < 20 { + return self.initializeByLongAlgorithm(&base, to: lhs, times: rhs, plus: Base.Element.zero) + } else { + return self.initializeByKaratsubaAlgorithm(&base, to: lhs, times: rhs) + } + } + + /// Initializes `base` to the square product of `elements`. + /// + /// - Parameters: + /// - base: A buffer of size `2 * elements.count`. + /// + /// - Note: The `base` memory must be uninitialized or the pointee must be a trivial type. + /// + @inlinable public static func initialize( + _ base: inout Base, toSquareProductOf elements: UnsafeBufferPointer) + where Base == UnsafeMutableBufferPointer { + //=--------------------------------------= + if elements.count < 20 { + return self.initializeByLongAlgorithm(&base, toSquareProductOf: elements, plus: Base.Element.zero) + } else { + return self.initializeByKaratsubaAlgorithm(&base, toSquareProductOf: elements) + } + } +} diff --git a/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Partition.swift b/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Partition.swift new file mode 100644 index 00000000..473a4afc --- /dev/null +++ b/Sources/NBKCoreKit/Private/NBKStrictUnsignedInteger+Partition.swift @@ -0,0 +1,31 @@ +//=----------------------------------------------------------------------------= +// This source file is part of the Numberick open source project. +// +// Copyright (c) 2023 Oscar Byström Ericsson +// Licensed under Apache License, Version 2.0 +// +// See http://www.apache.org/licenses/LICENSE-2.0 for license information. +//=----------------------------------------------------------------------------= + +//*============================================================================* +// MARK: * NBK x Strict Unsigned Integer x Partition x Sub Sequence +//*============================================================================* +//=----------------------------------------------------------------------------= +// MARK: + where Base is Unsafe Buffer Pointer +//=----------------------------------------------------------------------------= + +extension NBK.StrictUnsignedInteger.SubSequence { + + //=------------------------------------------------------------------------= + // MARK: Transformations + //=------------------------------------------------------------------------= + + /// Splits `base` at `index` then trims redundant zeros from each sequence. + @inlinable public static func partitionTrimmingRedundantZeros( + _ base: Base, at index: Base.Index) -> HL where Base == UnsafeBufferPointer { + let partition = Swift.min(base.count, index) + let low = Base(start: base.baseAddress!, /*------*/ count: NBK.dropLast(from: base[..