From ba414707e3bfc2c93ede586b1d5e5da59580b46f Mon Sep 17 00:00:00 2001 From: Michael Griebling Date: Thu, 13 Jul 2023 08:47:15 -0400 Subject: [PATCH] Fix bug in Burnikel-Ziegler division algorithm and update tests. --- Sources/BigInt/BigInt.swift | 4 +- Sources/BigInt/BurnikelZiegler.swift | 179 ++++++++++----------------- Tests/BigIntTests/DivModBZTest.swift | 48 ++++--- 3 files changed, 99 insertions(+), 132 deletions(-) diff --git a/Sources/BigInt/BigInt.swift b/Sources/BigInt/BigInt.swift index b8bfa34..ea4f036 100755 --- a/Sources/BigInt/BigInt.swift +++ b/Sources/BigInt/BigInt.swift @@ -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) } diff --git a/Sources/BigInt/BurnikelZiegler.swift b/Sources/BigInt/BurnikelZiegler.swift index 4fe46f2..5c25629 100644 --- a/Sources/BigInt/BurnikelZiegler.swift +++ b/Sources/BigInt/BurnikelZiegler.swift @@ -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) } } diff --git a/Tests/BigIntTests/DivModBZTest.swift b/Tests/BigIntTests/DivModBZTest.swift index bba6d60..100dffc 100644 --- a/Tests/BigIntTests/DivModBZTest.swift +++ b/Tests/BigIntTests/DivModBZTest.swift @@ -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) + } }