[BigInt] Implement division (Knuth).

Author: xwu

Sources/BigIntModule/BigInt._Significand.swift

@@ -661,14 +661,102 @@ extension BigInt._Significand {
 
   // @inlinable
   @discardableResult
-  internal func quotientAndRemainder(
-    dividingBy other: Self
-  ) -> (quotient: Self, remainder: Self) {
-    if other.count == 1 {
-      var x = self
-      let remainder = x.divide(by: other._low)
-      return (quotient: x, remainder: remainder)
+  internal mutating func divide(by other: Self) -> /* remainder: */ Self {
+    func shift(_ lhs: inout Self, left rhs: Int) {
+      var carry = 0 as UInt
+      guard rhs != 0 else { return }
+      for i in 0..<lhs.count {
+        let temporary = lhs[i]
+        lhs[i] = temporary &<< rhs | carry
+        carry = temporary &>> (UInt.bitWidth &- rhs)
+      }
+      lhs.append(carry)
+    }
+    
+    func shift(_ lhs: inout Self, right rhs: Int) {
+      var carry = 0 as UInt
+      guard rhs != 0 else { return }
+      for i in (0..<lhs.count).reversed() {
+        let temporary = lhs[i]
+        lhs[i] = temporary &>> rhs | carry
+        carry = temporary &<< (UInt.bitWidth &- rhs)
+      }
+    }
+    // Based on Knuth's Algorithm D.
+    // For details see <https://skanthak.homepage.t-online.de/division.html>.
+    
+    // We'll remove any extraneous leading zero words while determining by how
+    // much to shift our operands.
+    var other = other
+    var n = other.count
+    guard let i = other.lastIndex(where: { $0 != 0 }) else {
+      fatalError("Divide by zero")
     }
-    fatalError("Not implemented")
+    if n > i &+ 1 {
+      other.removeLast(n &- (i &+ 1))
+      n = i &+ 1
+    }
+    guard n > 1 else { return divide(by: other._low) }
+    let clz = other[n &- 1].leadingZeroBitCount
+    
+    var m = count - n
+    guard let j = lastIndex(where: { $0 != 0 }) else { return Self() }
+    if m > j &+ 1 {
+      removeLast(m &- (j &+ 1))
+      m = j &+ 1
+    }
+    precondition(m >= 0)
+    
+    // 1. Normalize operands.
+    shift(&other, left: clz)
+    shift(&self, left: clz)
+    let word = other[n &- 1]
+    
+    // 2. Initialize loop.
+    var result = Self(repeating: 0, count: m &+ 1)
+    for idx in (n...(n &+ m)).reversed() {
+      let (high, low) = (self[idx], self[idx &- 1])
+      
+      // 3. Calculate trial quotient and remainder.
+      var overflow = false
+      var (q̂, r̂): (UInt, UInt)
+      if word > high {
+        (q̂, r̂) = word.dividingFullWidth((high, low))
+      } else {
+        (q̂, r̂) = word < high
+          ? word.dividingFullWidth((high &- word, low))
+          : low.quotientAndRemainder(dividingBy: word)
+        overflow = true
+      }
+      while overflow || q̂.multipliedFullWidth(by: other[n &- 2]) > (r̂, low) {
+        q̂ &-= 1
+        r̂ &+= word
+        if r̂ < word { break }
+        if q̂ == UInt.max { overflow = false }
+      }
+      
+      // 4. Multiply and subtract.
+      let slice = self[(idx &- n)...idx]
+      var x = Self(slice), y = other
+      y.multiply(by: q̂)
+      overflow = x.subtract(y)
+      
+      // 5. Test remainder.
+      if overflow {
+        // 6. Add back.
+        q̂ -= 1
+        x.add(y)
+      }
+      
+      replaceSubrange((idx &- n)...idx, with: x[0...n])
+      result[idx &- n] = q̂
+    } // 7. Loop.
+    
+    // 8. Unnormalize.
+    removeLast(m &+ 1)
+    shift(&self, right: clz)
+    
+    swap(&self, &result)
+    return result
   }
 }

Sources/BigIntModule/BigInt.swift

@@ -390,48 +390,70 @@ extension BigInt: BinaryInteger {
   @inlinable
   public static var isSigned: Bool { true }
 
-  // @inlinable
+  @inlinable
   public static func / (lhs: BigInt, rhs: BigInt) -> BigInt {
-    guard rhs._combination != 0 else { fatalError("Division by zero") }
-    guard abs(rhs) <= abs(lhs) else { return 0 }
-    
-    //TODO: Implement division. For now, this is a skeleton placeholder.
-    if lhs._exponent == rhs._exponent {
-      let combination = lhs._signum * rhs._signum
-      let (significand, _) =
-        lhs._significand.quotientAndRemainder(dividingBy: rhs._significand)
-      var result = BigInt(_combination: combination, significand: significand)
-      result._normalize()
-      return result
-    }
-    fatalError("Not implemented")
+    var result = lhs
+    result /= rhs
+    return result
   }
 
   // @inlinable
   public static func /= (lhs: inout BigInt, rhs: BigInt) {
-    lhs = lhs / rhs
-  }
-
-  // @inlinable
-  public static func % (lhs: BigInt, rhs: BigInt) -> BigInt {
     guard rhs._combination != 0 else { fatalError("Division by zero") }
-    guard abs(rhs) <= abs(lhs) else { return lhs }
+    guard lhs._combination != 0 && abs(lhs) >= abs(rhs) else {
+      lhs = 0
+      return
+    }
+    
+    var rhs = rhs
+    let exponents = (lhs._exponent, rhs._exponent)
+    if exponents.0 < exponents.1 {
+      rhs._significand.insert(
+        contentsOf: repeatElement(0, count: exponents.1 &- exponents.0), at: 0)
+    } else if exponents.0 > exponents.1 {
+      lhs._significand.insert(
+        contentsOf: repeatElement(0, count: exponents.0 &- exponents.1), at: 0)
+    }
+    lhs._combination = lhs._signum * rhs._signum
 
-    //TODO: Implement remainder. For now, this is a skeleton placeholder.
-    if lhs._exponent == rhs._exponent {
-      let combination = lhs._combination * rhs._signum
-      let (_, significand) =
-        lhs._significand.quotientAndRemainder(dividingBy: rhs._significand)
-      var result = BigInt(_combination: combination, significand: significand)
-      result._normalize()
-      return result
+    if lhs._significand != rhs._significand {
+      lhs._significand.divide(by: rhs._significand)
+    } else {
+      lhs._significand = _Significand(1)
     }
-    fatalError("Not implemented")
+    lhs._normalize()
+  }
+
+  @inlinable
+  public static func % (lhs: BigInt, rhs: BigInt) -> BigInt {
+    var result = lhs
+    result %= rhs
+    return result
   }
 
   // @inlinable
   public static func %= (lhs: inout BigInt, rhs: BigInt) {
-    lhs = lhs % rhs
+    guard rhs._combination != 0 else { fatalError("Division by zero") }
+    guard lhs._combination != 0 && abs(lhs) >= abs(rhs) else { return }
+    
+    var rhs = rhs
+    let exponents = (lhs._exponent, rhs._exponent)
+    if exponents.0 < exponents.1 {
+      rhs._significand.insert(
+        contentsOf: repeatElement(0, count: exponents.1 &- exponents.0), at: 0)
+    } else if exponents.0 > exponents.1 {
+      lhs._significand.insert(
+        contentsOf: repeatElement(0, count: exponents.0 &- exponents.1), at: 0)
+    }
+    lhs._combination = lhs._signum
+    
+    if lhs._significand != rhs._significand {
+      var result = lhs._significand.divide(by: rhs._significand)
+      swap(&result, &lhs._significand)
+    } else {
+      lhs._significand = _Significand()
+    }
+    lhs._normalize()
   }
 
   // @inlinable

Tests/BigIntTests/BigIntTests.swift

@@ -165,6 +165,33 @@ final class BigIntModuleTests: XCTestCase {
     }
   }
 
+  func testDivision() {
+    let x = BigInt("327441402998268901582239630362983195839")!
+    let y = BigInt("279240677930711642307518656231691197860")!
+    XCTAssertEqual(x / y, 1)
+
+    for _ in 0..<50 {
+      let a = _randomWords(count: 6)
+      let b = _randomWords(count: 4)
+      
+      XCTAssertEqual(a.0 / b.0, BigInt(a.1 / b.1))
+      XCTAssertEqual(a.0 % b.0, BigInt(a.1 % b.1))
+      XCTAssertEqual(
+        (a.0 << 128) / (b.0 << 42),
+        BigInt((a.1 << 128) / (b.1 << 42)))
+      XCTAssertEqual(
+        (a.0 << 128) % (b.0 << 42),
+        BigInt((a.1 << 128) % (b.1 << 42)))
+      XCTAssertEqual(
+        (a.0 << 42) / (b.0 << 128),
+        BigInt((a.1 << 42) / (b.1 << 128)))
+      XCTAssertEqual(
+        (a.0 << 42) % (b.0 << 128),
+        BigInt((a.1 << 42) % (b.1 << 128)))
+      XCTAssertEqual((a.0 << 128) / (b.0 << 128), a.0 / b.0)
+    }
+  }
+  
   func testBitwiseOperators() {
     for _ in 0..<100 {
       let a = _randomWords(count: 8)