Skip to content

Commit

Permalink
Fix bug in Burnikel-Ziegler division algorithm and update tests.
Browse files Browse the repository at this point in the history
  • Loading branch information
@mgriebling
mgriebling committed Jul 13, 2023
1 parent cbc7ae1 commit ba41470
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 132 deletions.
4 changes: 2 additions & 2 deletions Sources/BigInt/BigInt.swift
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1139,8 +1139,8 @@ extension BInt {
public func quotientAndRemainder(dividingBy x: BInt) -> (quotient: BInt, remainder: BInt) {
var quotient = BInt.zero
var remainder = BInt.zero
if x.limbs.count > Limbs.BZ_DIV_LIMIT && self.limbs.count > x.limbs.count + Limbs.BZ_DIV_LIMIT {
(quotient.limbs, remainder.limbs) = self.limbs.bzDivMod(x.limbs)
if x.limbs.count > BInt.BZ_DIV_LIMIT && self.limbs.count > x.limbs.count + BInt.BZ_DIV_LIMIT {
(quotient, remainder) = self.bzDivMod(x)
} else {
(quotient.limbs, remainder.limbs) = self.limbs.divMod(x.limbs)
}
Expand Down
179 changes: 68 additions & 111 deletions Sources/BigInt/BurnikelZiegler.swift
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,146 +5,103 @@
// Created by Leif Ibsen on 05/12/2021.
//


extension Array where Element == Limb {

extension BInt {

// Divisor limb limit for Burnikel-Ziegler division
static let BZ_DIV_LIMIT = 60
static var BZ_DIV_LIMIT = 60

/*
* [BURNIKEL] - algorithm 3
*/
func bzDivMod(_ v: Limbs) -> (quotient: Limbs, remainder: Limbs) {
var quotient: Limbs = [0]
var remainder: Limbs = []
var A = self
var B = v
let s = B.count
let m = 1 << (64 - (s / Limbs.BZ_DIV_LIMIT).leadingZeroBitCount)
func bzDivMod(_ v: BInt) -> (q: BInt, r: BInt) {
var A = self.abs
var B = v.abs
let s = B.limbs.count
let m = 1 << (64 - (s / BInt.BZ_DIV_LIMIT).leadingZeroBitCount)
let j = (s + m - 1) / m
let n = j * m
let n64 = n << 6
let sigma = Swift.max(0, n64 - B.bitWidth)
A.shiftLeft(sigma)
B.shiftLeft(sigma)
A <<= sigma
B <<= sigma
let t = Swift.max(2, (A.bitWidth + n64) / n64)
var Z = Limbs(repeating: 0, count: 2 * n)
var from = (t - 1) * n
var zi = n
for ai in from ..< A.count {
Z[zi] = A[ai]
zi += 1
}
from -= n
zi = 0
for ai in from ..< from + n {
Z[zi] = A[ai]
zi += 1
}
var Q = BInt.zero
var R = BInt.zero
var Z = A.blockA(n, t - 1) << n64 | A.blockA(n, t - 2)
for i in (0 ... t - 2).reversed() {
var (Q, R) = Div2n1n(n, Z, B)
R.normalize()
quotient.add(Q, from)
var Qi: BInt
(Qi, R) = Z.bzDiv2n1n(n, B)
Q = Q << n64 | Qi
if i > 0 {
from -= n
for zi in 0 ..< R.count {
Z[n + zi] = R[zi]
}
for zi in R.count ..< n {
Z[n + zi] = 0
}
zi = 0
for ai in from ..< from + n {
Z[zi] = A[ai]
zi += 1
}
} else {
remainder = R
remainder.shiftRight(sigma)
Z = R << n64 | A.blockA(n, i - 1)
}
}
return (quotient, remainder)
return (Q, R >> sigma)

}

func blockA(_ n: Int, _ i: Int) -> BInt {
let mc = self.limbs.count
assert(i * n <= mc)
if (i + 1) * n < mc {
return BInt(Limbs(self.limbs[i * n ..< (i + 1) * n]))
} else {
return BInt(Limbs(self.limbs[i * n ..< mc]))
}
}

/*
* [BURNIKEL] - algorithm 1
*/
func Div2n1n(_ n: Int, _ A: Limbs, _ B: Limbs) -> (Limbs, Limbs) {
if B.count & 1 == 1 || B.count < Limbs.BZ_DIV_LIMIT {
func bzDiv2n1n(_ n: Int, _ B: BInt) -> (q: BInt, r: BInt) {
if n & 1 == 1 || B.limbs.count < BInt.BZ_DIV_LIMIT {

// Basecase

var a = A
a.normalize()
var b = B
b.normalize()
return a.divMod(b)
}
var A1: Limbs
var A2: Limbs
var A3: Limbs
var A4: Limbs
let n12 = n >> 1
let n32 = 3 * n12
if A.count > n32 {
A1 = Limbs(A[n32 ..< A.count])
A2 = Limbs(A[n ..< n32])
A3 = Limbs(A[n12 ..< n])
A4 = Limbs(A[0 ..< n12])
} else if A.count > n {
A1 = [0]
A2 = Limbs(A[n ..< A.count])
A3 = Limbs(A[n12 ..< n])
A4 = Limbs(A[0 ..< n12])
} else if A.count > n12 {
A1 = [0]
A2 = [0]
A3 = Limbs(A[n12 ..< A.count])
A4 = Limbs(A[0 ..< n12])
} else {
A1 = [0]
A2 = [0]
A3 = [0]
A4 = Limbs(A[0 ..< A.count])
// Base case

let (q, r) = self.limbs.divMod(B.limbs)
return (BInt(q), BInt(r))
}
let (Q1, R1) = Div3n2n(n12, A1, A2, A3, B)
let R11 = Limbs(R1.count > n12 ? R1[n12 ..< R1.count] : [0])
let R12 = Limbs(R1.count > n12 ? R1[0 ..< n12] : R1[0 ..< R1.count])
let (Q2, R) = Div3n2n(n12, R11, R12, A4, B)
var Q = Q1.shiftedLeft(n12 << 6)
Q.add(Q2, 0)
return (Q, R)
let n2 = n >> 1
let n32 = n << 5
let (A123, A4) = self.split(n2)
let (Q1, R1) = A123.bzDiv3n2n(n2, B)
let (R11, R12) = R1.split(n)
let (Q2, R) = ((R11 << n32 | R12) << n32 | A4).bzDiv3n2n(n2, B)
return (Q1 << n32 + Q2, R)
}

/*
* [BURNIKEL] - algorithm 2
*/
func Div3n2n(_ n: Int, _ A1: Limbs, _ A2: Limbs, _ A3: Limbs, _ B: Limbs) -> (Limbs, Limbs) {
let B1 = Limbs(B.count > n ? B[n ..< B.count] : [0])
let B2 = Limbs(B.count > n ? B[0 ..< n] : B[0 ..< B.count])
var Q: Limbs
var R1: Limbs
if A1.compare(B1) < 0 {
var A = A1.shiftedLeft(n << 6)
A.add(A2, 0)
(Q, R1) = Div2n1n(n, A, B1)
func bzDiv3n2n(_ n: Int, _ B: BInt) -> (q: BInt, r: BInt) {
let n64 = n << 6
let (B1, B2) = B.split(n)
let (A12, A3) = self.split(n)
let (A1, _) = A12.split(n)
var Q: BInt
var R1: BInt
if A1 < B1 {
(Q, R1) = A12.bzDiv2n1n(n, B1)
} else {
R1 = A1
_ = R1.subtract(B1, 0)
R1.shiftLeft(n << 6)
R1.add(B1, 0)
Q = Limbs(repeating: 0xffffffffffffffff, count: n)
R1 = A12 - (B1 << n64) + B1
Q = BInt.one << n64 - BInt.one
}
let D = Q * B2
var R = R1 << n64 | A3
while R < D {
R += B
Q -= BInt.one
}
var D = Q
D.multiply(B2)
var R = R1.shiftedLeft(n << 6)
R.add(A3)
while R.compare(D) < 0 {
R.add(B)
_ = Q.subtract([1], 0)
return (Q, R - D)
}

func split(_ n: Int) -> (BInt, BInt) {
let mc = self.limbs.count
if mc > n {
return (BInt(Limbs(self.limbs[n ..< mc])), BInt(Limbs(self.limbs[0 ..< n])))
} else {
return (BInt.zero, self)
}
_ = R.subtract(D, 0)
return (Q, R)
}

}
48 changes: 29 additions & 19 deletions Tests/BigIntTests/DivModBZTest.swift
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,47 +14,57 @@ import XCTest

class DivModBZTest: XCTestCase {

override func setUpWithError() throws {
// Put setup code here. This method is called before the invocation of each test method in the class.
}

override func tearDownWithError() throws {
// Put teardown code here. This method is called after the invocation of each test method in the class.
}
let b1 = BInt("4782202788054612029528392986600059097414971724022365008513345109918378950942662970278927686112707894586824720981524256319306585052676834087480834429433264797425893247623688331021633208954847354805799943341309825989013743806187109581043148680813778321530496715601563282624414040398143207622036272190408590790537203475256105564071579263867875240985573356522656108542128577321057879052328865035355873615679363655889925711574420153832091752422843046918811427400662135559303516853703976812686385750376227787949580582081831261725701003498206512329872677233489510953469375683037038373999696771585788905639115522613405495707184524158219208223766442059014593330657009722153962376853423770486138578089775621301167811299166407361746606697808186757966914671246073712904200588408923186387737887675292886953797066980967406053530122853539036965490224784924649007954898678503314655546475504501686187354866964374552614120640782949622452027788962138602665933147687696322089504278791624651519312327831756553779377194524673395819281486668576384019590720179413349582970319393884388810494546040342087536563628332152073181614300721769371426238517540520845214665313301183551962591849558938499025348780376716477073930634436840084468255937443451690315999349137664638968972614199015304906547819056227171224947070739716300953775743441307920501863532234466545645695774331885044978250148663467372130392099894852145190998232878772486650513010816769902892518719250066947215706536216248696237476619548294280997466823267797934620158926088915926052113157320399149350072309099187966099088189643497513687928960839375289097070057698545052161834612441440300447910931798424330069494272576133664129571446005941039050449109761662280133298161365994610326139093159039351518411873726699434670218627794470163865731797258607504515467706568302430378820423176019440925413294968918946753234527062907457039327403799514654681482620563976230609510390120177018153378953914211667001350128690497063627989286291690626363431460266484993616898632404483747215192933193062784093302610061072941730610833278036184381530158256539585838588799897714285126286720260761560779987177142728853496833739616878241477395730921177557832273916603145310867588335248530358959095181725126608131193726511914438842826945650698812108368961639111997274610216346844001820682880739430341057410931474128548720115358991355261589298277633257899790971197558786178367745055254357923445045308508288059601997158296253020130655356530996489598877891324128600243475505599610596315196496323542179888792209304804935537197649858327858787700693018393896295993380379016269866179084654381355939982578243521489624917379652576570974212281051937534650720271599262478255593167843763594436777566708623653325250511141415367412037101277051776800861048051216123083433710443499239892012310615341329522047064721158043794002611244331208099097510130494454107820850957462322893961198397453714743859503326095545071195512880003809279")
let b2 = BInt("3069183075406921794879924958191326267980083732617257090967156575624265074992632621684985264011687141177463521361517911553021245099194483198628542584486348563466378125452177821735146139462477045607889665505906919596910519383246040059165099934889611077055111196533329299014397179377993006736606277358858761491911450868432133382295467939853925371563407906143030316759860088351677624708599821169550833590036545299179057356555493184183038259794167997381510402990078867498040076511934068262113082029079685194092585200557212952087870250306201570762147399263274328246438930899678576215096246269986660643586905677861972798413323900862279393294098586674450767550167631066398651798197094996096391423300063966130361456553414900428770512722081658104476364858583339960519510499984305492480977699199899564863730624237563221068841591071268000190139910882859705395243310564643593059495182047032409956901498079495621606539723949410789406025669661263490770568454509838497464363714458645954467110174994932482618754701038967174720237992450773625322184776742653465690932102840101388417452307232873785183975126764865419953046185303713655389674631027509028941278698713684303566707778368015911611508474550527223037689855")

func doTest1(_ dividend: BInt, _ divisor: BInt) {
var q1: Limbs = []
var r1: Limbs = []
var q2: Limbs = []
var r2: Limbs = []
var q2: BInt
var r2: BInt
(q1, r1) = dividend.limbs.divMod(divisor.limbs)
(q2, r2) = dividend.limbs.bzDivMod(divisor.limbs)
XCTAssertEqual(q1, q2)
XCTAssertEqual(r1, r2)
(q2, r2) = dividend.bzDivMod(divisor)
XCTAssertEqual(q1, q2.limbs)
XCTAssertEqual(r1, r2.limbs)
}

func doTest2(_ dividend: BInt, _ divisor: BInt) {
let (q, r) = dividend.quotientAndRemainder(dividingBy: divisor)
XCTAssertEqual(dividend, divisor * q + r)
}

func test1() {
for _ in 0 ..< 1000 {
let x = BInt(bitWidth: 2 * (Limbs.BZ_DIV_LIMIT + 1) * 64)
let y = BInt(bitWidth: (Limbs.BZ_DIV_LIMIT + 1) * 64)
let x = BInt(bitWidth: 2 * (BInt.BZ_DIV_LIMIT + 1) * 64)
let y = BInt(bitWidth: (BInt.BZ_DIV_LIMIT + 1) * 64)
doTest1(x, y)
doTest2(x, y)
doTest2(x, -y)
doTest2(-x, y)
doTest2(-x, -y)
}
doTest1(BInt(Limbs(repeating: UInt64.max, count: 2 * (Limbs.BZ_DIV_LIMIT + 1))), BInt(Limbs(repeating: UInt64.max, count: Limbs.BZ_DIV_LIMIT + 1)))
doTest1(BInt(Limbs(repeating: UInt64.max, count: 2 * (BInt.BZ_DIV_LIMIT + 1))), BInt(Limbs(repeating: UInt64.max, count: BInt.BZ_DIV_LIMIT + 1)))
}

func test2() {
var n = 2
for _ in 0 ..< 8 {
let y1 = BInt.one << ((Limbs.BZ_DIV_LIMIT + 1) * 64 * n)
let y2 = BInt(bitWidth: (Limbs.BZ_DIV_LIMIT + 1) * 64 * n)
let y1 = BInt.one << ((BInt.BZ_DIV_LIMIT + 1) * 64 * n)
let y2 = BInt(bitWidth: (BInt.BZ_DIV_LIMIT + 1) * 64 * n)
n *= 2
let x1 = BInt.one << ((Limbs.BZ_DIV_LIMIT + 1) * 64 * n)
let x2 = BInt(bitWidth: (Limbs.BZ_DIV_LIMIT + 1) * 64 * n)
let x1 = BInt.one << ((BInt.BZ_DIV_LIMIT + 1) * 64 * n)
let x2 = BInt(bitWidth: (BInt.BZ_DIV_LIMIT + 1) * 64 * n)
doTest1(x1, y1)
doTest1(x1, y2)
doTest1(x2, y1)
doTest1(x2, y2)
}
}

// This test used to fail in release 1.13.0
func test3() {
doTest2(b1, b1 * b1)
doTest2(b2, b2 * b2)
}

}

0 comments on commit ba41470

Please sign in to comment.