Skip to content

Commit

Permalink
Release 1.8.0
Browse files Browse the repository at this point in the history
  • Loading branch information
@leif-ibsen
leif-ibsen committed Mar 22, 2023
1 parent 95a21a6 commit 2b2f7f4
Show file tree
Hide file tree
Showing 15 changed files with 273 additions and 160 deletions.
70 changes: 35 additions & 35 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ BigInt requires Swift 5.0. It also requires that the Int and UInt types be 64 bi
In your projects Package.swift file add a dependency like<br/>

dependencies: [
.package(url: "https://github.com/leif-ibsen/BigInt", from: "1.7.0"),
.package(url: "https://github.com/leif-ibsen/BigInt", from: "1.8.0"),
]

<h2 id="ex"><b>Examples</b></h2>
Expand Down Expand Up @@ -90,49 +90,49 @@ Four large numbers 'a1000', 'b1000', 'c2000' and 'p1000' were used throughout th

<table width="90%">
<tr><th width="40%" align="left">Operation</th><th width="35%" align="right">Swift code</th><th width="25%" align="right">Time</th></tr>
<tr><td>As string</td><td align="right">c2000.asString()</td><td align="right">15 uSec</td></tr>
<tr><td>As signed bytes</td><td align="right">c2000.asSignedBytes()</td><td align="right">0.27 uSec</td></tr>
<tr><td>Bitwise and</td><td align="right">a1000 & b1000</td><td align="right">0.080 uSec</td></tr>
<tr><td>Bitwise or</td><td align="right">a1000 | b1000</td><td align="right">0.080 uSec</td></tr>
<tr><td>Bitwise xor</td><td align="right">a1000 ^ b1000</td><td align="right">0.080 uSec</td></tr>
<tr><td>Bitwise not</td><td align="right">~c2000</td><td align="right">0.084 uSec</td></tr>
<tr><td>Test bit</td><td align="right">c2000.testBit(701)</td><td align="right">0.0045 uSec</td></tr>
<tr><td>Flip bit</td><td align="right">c2000.flipBit(701)</td><td align="right">0.0072 uSec</td></tr>
<tr><td>Set bit</td><td align="right">c2000.setBit(701)</td><td align="right">0.0066 uSec</td></tr>
<tr><td>Clear bit</td><td align="right">c2000.clearBit(701)</td><td align="right">0.0063 uSec</td></tr>
<tr><td>As string</td><td align="right">c2000.asString()</td><td align="right">13 uSec</td></tr>
<tr><td>As signed bytes</td><td align="right">c2000.asSignedBytes()</td><td align="right">0.31 uSec</td></tr>
<tr><td>Bitwise and</td><td align="right">a1000 & b1000</td><td align="right">0.087 uSec</td></tr>
<tr><td>Bitwise or</td><td align="right">a1000 | b1000</td><td align="right">0.086 uSec</td></tr>
<tr><td>Bitwise xor</td><td align="right">a1000 ^ b1000</td><td align="right">0.087 uSec</td></tr>
<tr><td>Bitwise not</td><td align="right">~c2000</td><td align="right">0.098 uSec</td></tr>
<tr><td>Test bit</td><td align="right">c2000.testBit(701)</td><td align="right">0.018 uSec</td></tr>
<tr><td>Flip bit</td><td align="right">c2000.flipBit(701)</td><td align="right">0.019 uSec</td></tr>
<tr><td>Set bit</td><td align="right">c2000.setBit(701)</td><td align="right">0.019 uSec</td></tr>
<tr><td>Clear bit</td><td align="right">c2000.clearBit(701)</td><td align="right">0.019 uSec</td></tr>
<tr><td>Addition</td><td align="right">a1000 + b1000</td><td align="right">0.10 uSec</td></tr>
<tr><td>Subtraction</td><td align="right">a1000 - b1000</td><td align="right">0.11 uSec</td></tr>
<tr><td>Multiplication</td><td align="right">a1000 * b1000</td><td align="right">0.32 uSec</td></tr>
<tr><td>Division</td><td align="right">c2000 / a1000</td><td align="right">2.5 uSec</td></tr>
<tr><td>Modulus</td><td align="right">c2000.mod(a1000)</td><td align="right">2.5 uSec</td></tr>
<tr><td>Inverse modulus</td><td align="right">c2000.modInverse(p1000)</td><td align="right">0.39 mSec</td></tr>
<tr><td>Modular exponentiation</td><td align="right">a1000.expMod(b1000, c2000)</td><td align="right">3.8 mSec</td></tr>
<tr><td>Equal</td><td align="right">c2000 + 1 == c2000</td><td align="right">0.0016 uSec</td></tr>
<tr><td>Less than</td><td align="right">b1000 + 1 < b1000</td><td align="right">0.010 uSec</td></tr>
<tr><td>Shift 1 left</td><td align="right">c2000 <<= 1</td><td align="right">0.15 uSec</td></tr>
<tr><td>Shift 1 right</td><td align="right">c2000 >>= 1</td><td align="right">0.11 uSec</td></tr>
<tr><td>Shift 100 left</td><td align="right">c2000 <<= 100</td><td align="right">0.19 uSec</td></tr>
<tr><td>Shift 100 right</td><td align="right">c2000 >>= 100</td><td align="right">0.15 uSec</td></tr>
<tr><td>Is probably prime</td><td align="right">p1000.isProbablyPrime()</td><td align="right">6.2 mSec</td></tr>
<tr><td>Make probable 1000 bit prime</td><td align="right">BInt.probablePrime(1000)</td><td align="right">63 mSec</td></tr>
<tr><td>Next probable prime</td><td align="right">c2000.nextPrime()</td><td align="right">790 mSec</td></tr>
<tr><td>Inverse modulus</td><td align="right">c2000.modInverse(p1000)</td><td align="right">110 uSec</td></tr>
<tr><td>Modular exponentiation</td><td align="right">a1000.expMod(b1000, c2000)</td><td align="right">3.7 mSec</td></tr>
<tr><td>Equal</td><td align="right">c2000 + 1 == c2000</td><td align="right">0.017 uSec</td></tr>
<tr><td>Less than</td><td align="right">b1000 + 1 < b1000</td><td align="right">0.021 uSec</td></tr>
<tr><td>Shift 1 left</td><td align="right">c2000 <<= 1</td><td align="right">0.10 uSec</td></tr>
<tr><td>Shift 1 right</td><td align="right">c2000 >>= 1</td><td align="right">0.07 uSec</td></tr>
<tr><td>Shift 100 left</td><td align="right">c2000 <<= 100</td><td align="right">0.15 uSec</td></tr>
<tr><td>Shift 100 right</td><td align="right">c2000 >>= 100</td><td align="right">0.13 uSec</td></tr>
<tr><td>Is probably prime</td><td align="right">p1000.isProbablyPrime()</td><td align="right">6.1 mSec</td></tr>
<tr><td>Make probable 1000 bit prime</td><td align="right">BInt.probablePrime(1000)</td><td align="right">60 mSec</td></tr>
<tr><td>Next probable prime</td><td align="right">c2000.nextPrime()</td><td align="right">770 mSec</td></tr>
<tr><td>Primorial</td><td align="right">BInt.primorial(100000)</td><td align="right">16 mSec</td></tr>
<tr><td>Binomial</td><td align="right">BInt.binomial(100000, 10000)</td><td align="right">26 mSec</td></tr>
<tr><td>Factorial</td><td align="right">BInt.factorial(100000)</td><td align="right">69 mSec</td></tr>
<tr><td>Factorial</td><td align="right">BInt.factorial(100000)</td><td align="right">68 mSec</td></tr>
<tr><td>Fibonacci</td><td align="right">BInt.fibonacci(100000)</td><td align="right">0.74 mSec</td></tr>
<tr><td>Greatest common divisor</td><td align="right">a1000.gcd(b1000)</td><td align="right">31 uSec</td></tr>
<tr><td>Extended gcd</td><td align="right">a1000.gcdExtended(b1000)</td><td align="right">0.59 mSec</td></tr>
<tr><td>Least common multiple</td><td align="right">a1000.lcm(b1000)</td><td align="right">31 uSec</td></tr>
<tr><td>Make random number</td><td align="right">c2000.randomLessThan()</td><td align="right">1.3 uSec</td></tr>
<tr><td>Square</td><td align="right">c2000 ** 2</td><td align="right">0.80 uSec</td></tr>
<tr><td>Greatest common divisor</td><td align="right">a1000.gcd(b1000)</td><td align="right">32 uSec</td></tr>
<tr><td>Extended gcd</td><td align="right">a1000.gcdExtended(b1000)</td><td align="right">90 uSec</td></tr>
<tr><td>Least common multiple</td><td align="right">a1000.lcm(b1000)</td><td align="right">35 uSec</td></tr>
<tr><td>Make random number</td><td align="right">c2000.randomLessThan()</td><td align="right">1.2 uSec</td></tr>
<tr><td>Square</td><td align="right">c2000 ** 2</td><td align="right">0.83 uSec</td></tr>
<tr><td>Square root</td><td align="right">c2000.sqrt()</td><td align="right">14 uSec</td></tr>
<tr><td>Square root and remainder</td><td align="right">c2000.sqrtRemainder()</td><td align="right">14 uSec</td></tr>
<tr><td>Is perfect square</td><td align="right">(c2000 * c2000).isPerfectSquare()</td><td align="right">18 uSec</td></tr>
<tr><td>Square root modulo</td><td align="right">b1000.sqrtMod(p1000)</td><td align="right">1.7 mSec</td></tr>
<tr><td>Power</td><td align="right">c2000 ** 111</td><td align="right">2.3 mSec</td></tr>
<tr><td>Root</td><td align="right">c2000.root(111)</td><td align="right">17 uSec</td></tr>
<tr><td>Root</td><td align="right">c2000.root(111)</td><td align="right">16 uSec</td></tr>
<tr><td>Root and remainder</td><td align="right">c2000.rootRemainder(111)</td><td align="right">18 uSec</td></tr>
<tr><td>Is perfect root</td><td align="right">c2000.isPerfectRoot()</td><td align="right">17 mSec</td></tr>
<tr><td>Is perfect root</td><td align="right">c2000.isPerfectRoot()</td><td align="right">14 mSec</td></tr>
<tr><td>Jacobi symbol</td><td align="right">c2000.jacobiSymbol(p1000)</td><td align="right">0.15 mSec</td></tr>
<tr><td>Kronecker symbol</td><td align="right">c2000.kroneckerSymbol(p1000)</td><td align="right">0.15 mSec</td></tr>
</table>
Expand All @@ -155,9 +155,9 @@ Fractions are created by
<li>Specifying the numerator and denominator explicitly f.ex. BFraction(17, 4)</li>
<li>Specifying the decimal value explictly f.ex. BFraction(4.25)</li>
</ul>
Defining a fraction by giving its decimal value (like 4.25) migth lead to surprises,
Defining a fraction by giving its decimal value (like 4.25) might lead to surprises,
because not all decimal values can be represented exactly as a floating point number.
For example, one migth think that BFraction(0.1) would equal 1/10,
For example, one might think that BFraction(0.1) would equal 1/10,
but in fact it equals 3602879701896397 / 36028797018963968 = 0.1000000000000000055511151231257827021181583404541015625</br>
<h3><b>Operations</b></h3>
The operations available to fractions are:
Expand Down Expand Up @@ -194,7 +194,7 @@ There are references in the source code where appropriate.
<li>[HACKER] - Henry S. Warren, Jr.: Hacker's Delight. Second Edition, Addison-Wesley</li>
<li>[HANDBOOK] - Menezes, Oorschot, Vanstone: Handbook of Applied Cryptography. CRC Press 1996</li>
<li>[JEBELEAN] - Tudor Jebelean: An Algorithm for Exact Division. Journal of Symbolic Computation, volume 15, 1993</li>
<li>[KNUTH] - Donald E. Knuth: Seminumerical Algorithms. Addison-Wesley 1971</li>
<li>[KNUTH] - Donald E. Knuth: Seminumerical Algorithms, Third Edition</li>
</ul>

<h2 id="alg"><b>Algorithms</b></h2>
Expand All @@ -212,13 +212,13 @@ Some of the algorithms used in BigInt are described below.
<li>Basecase - Knuth algorithm D</li>
<li>Exact Division - Jebelean's exact division algorithm</li>
</ul>
<h3><b>Greatest Common Divisor</b></h3>
<h3><b>Greatest Common Divisor and Extended Greatest Common Divisor</b></h3>
Lehmer's algorithm [KNUTH] chapter 4.5.2, with binary GCD basecase.
<h3><b>Modular Exponentiation</b></h3>
Sliding window algorithm 14.85 from [HANDBOOK] using Barrett reduction for exponents with fewer than 2048 bits
and Montgomery reduction for larger exponents.
<h3><b>Inverse Modulus</b></h3>
Extended Euclid algorithm 2.1.4 from [CRANDALL].
Computed via extended GCD algorithm.
<h3><b>Square Root</b></h3>
Algorithm 1.12 (SqrtRem) from [BRENT] with algorithm 9.2.11 from [CRANDALL] as basecase.
<h3><b>Square Root Modulo a Prime Number</b></h3>
Expand Down
141 changes: 96 additions & 45 deletions Sources/BigInt/BigInt.swift
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,12 +34,6 @@ infix operator ** : ExponentiationPrecedence
/// The representation is minimal, there is no leading zero Limbs.
/// The exception is that the value 0 is represented as a single 64 bit zero Limb and sign *false*
public struct BInt: CustomStringConvertible, Comparable, Equatable, Hashable {

static let digits: [Character] = [
"0", "1", "2", "3", "4", "5", "6", "7", "8", "9",
"a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z",
"A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z"]


// MARK: - Constants

Expand Down Expand Up @@ -141,7 +135,6 @@ public struct BInt: CustomStringConvertible, Comparable, Equatable, Hashable {
if number.isEmpty {
return nil
}
var magnitude = [Limb(0)]

// Find the number of digits that fits in a single Limb for the given radix
// Process that number of digits at a time
Expand All @@ -158,26 +151,26 @@ public struct BInt: CustomStringConvertible, Comparable, Equatable, Hashable {
if g == 0 {
g = digits
}
var magnitude = [Limb(0)]
var i = 0
var l = Limb(0)
for c in number {
if let digit = BInt.digits.firstIndex(of: c) {
let d = digit < 36 ? digit : digit - 26
if d >= radix {
return nil
}
l *= Limb(radix)
l += Limb(d)
} else {
return nil
}
var from = number.startIndex
var to = from
while from != number.endIndex {
to = number.index(after: to)
i += 1
if i == g {
magnitude.multiply(pow)
guard let l = Limb(number[from ..< to], radix: radix) else {
return nil
}
if radix == 16 {
magnitude.shiftLeft(60)
} else {
magnitude.multiply(pow)
}
magnitude.add(l)
g = digits
l = 0
i = 0
from = to
}
}
self.init(magnitude, sign)
Expand Down Expand Up @@ -1212,9 +1205,7 @@ public struct BInt: CustomStringConvertible, Comparable, Equatable, Hashable {
}

/*
* [CRANDALL] - algorithm 2.1.4
*
* Return self modinverse m
* Return self modinverse m - using the extended gcd
*/
/// Inverse modulus - BInt parameter
///
Expand All @@ -1223,16 +1214,9 @@ public struct BInt: CustomStringConvertible, Comparable, Equatable, Hashable {
/// - Returns: If *self* and m are coprime, x such that (*self* * x) mod m = 1
public func modInverse(_ m: BInt) -> BInt {
precondition(m.isPositive, "Modulus must be positive")
var a = BInt.ONE
var g = self.mod(m)
var u = BInt.ZERO
var w = m
while w.isPositive {
let (q, r) = g.quotientAndRemainder(dividingBy: w)
(a, g, u, w) = (u, w, a - q * u, r)
}
let (g, a, _) = self.gcdExtended(m)
precondition(g.isOne, "Modulus and self are not coprime")
return a.isNegative ? a + m : a
return a.mod(m)
}

/// Inverse modulus - Int parameter
Expand Down Expand Up @@ -2311,24 +2295,91 @@ public struct BInt: CustomStringConvertible, Comparable, Equatable, Hashable {
}

/*
* [CRANDALL] - algorithm 2.1.4
* Lehmer's gcd algorithm
* [KNUTH] - chapter 4.5.2, algorithm L - exercise 18
*/
/// Extended greatest common divisor
///
/// - Parameter x: Operand
/// - Returns: Greatest common divisor *g* of *self* and *x*, and *a* and *b* such that *a* * *self* + *b* * *x* = *g*
public func gcdExtended(_ x: BInt) -> (g: BInt, a: BInt, b: BInt) {
var a = BInt.ONE
var b = BInt.ZERO
var g = self
var u = BInt.ZERO
var v = BInt.ONE
var w = x
while w.isNotZero {
let q = g / w
(a, b, g, u, v, w) = (u, v, w, a - q * u, b - q * v, g - q * w)
}
return g.isNegative ? (-g, -a, -b) : (g, a, b)
let selfabs = self.abs
let xabs = x.abs
if self.isZero {
return (xabs, BInt.ZERO, x.isNegative ? -BInt.ONE : BInt.ONE)
}
if x.isZero {
return (selfabs, self.isNegative ? -BInt.ONE : BInt.ONE, BInt.ZERO)
}
var u: BInt
var v: BInt
var u2 = BInt.ZERO
var v2 = BInt.ONE
if selfabs < xabs {
u = xabs
v = selfabs
} else {
u = selfabs
v = xabs
}
var u3 = u
var v3 = v
while v >= BInt.B62 {
let size = u.bitWidth - 62
var x = (u >> size).asInt()!
var y = (v >> size).asInt()!
var A = 1
var B = 0
var C = 0
var D = 1
while true {
let yC = y + C
let yD = y + D
if yC == 0 || yD == 0 {
break
}
let q = (x + A) / yC
if q != (x + B) / yD {
break
}
(A, B, x, C, D, y) = (C, D, y, A - q * C, B - q * D, x - q * y)
}
if B == 0 {
(u, v) = (v, u.mod(v))
let q = u3 / v3
(u2, v2) = (v2, u2 - q * v2)
(u3, v3) = (v3, u3 - q * v3)
} else {
(u, v) = (A * u + B * v, C * u + D * v)
(u2, v2) = (A * u2 + B * v2, C * u2 + D * v2)
(u3, v3) = (A * u3 + B * v3, C * u3 + D * v3)
}
}
while v3.isNotZero {
let q = u3 / v3
(u2, v2) = (v2, u2 - v2 * q)
(u3, v3) = (v3, u3 - v3 * q)
}

var u1: BInt
if selfabs < xabs {
// u3 = u2 * selfabs + u1 * xabs
u1 = (u3 - u2 * selfabs) / xabs
(u1, u2) = (u2, u1)
} else {
// u3 = u1 * selfabs + u2 * xabs
u1 = (u3 - u2 * xabs) / selfabs
}

// Fix the signs

if x.isNegative {
u2 = -u2
}
if self.isNegative {
u1 = -u1
}
return (u3, u1, u2)
}

/*
Expand Down
61 changes: 61 additions & 0 deletions Tests/BigIntTests/ConstructorTest.swift
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,4 +169,65 @@ class ConstructorTest: XCTestCase {
XCTAssertTrue((x - x1).abs.asDouble() / d <= 1.0e-15)
}
}

let strings = ["0",
"10",
"120",
"1230",
"12340",
"123450",
"1234560",
"12345670",
"123456780",
"1234567890",
"123456789a0",
"123456789ab0",
"123456789abc0",
"123456789abcd0",
"123456789abcde0",
"123456789abcdef0",
"123456789abcdefg0",
"123456789abcdefgh0",
"123456789abcdefghi0",
"123456789abcdefghij0",
"123456789abcdefghijk0",
"123456789abcdefghijkl0",
"123456789abcdefghijklm0",
"123456789abcdefghijklmn0",
"123456789abcdefghijklmno0",
"123456789abcdefghijklmnop0",
"123456789abcdefghijklmnopq0",
"123456789abcdefghijklmnopqr0",
"123456789abcdefghijklmnopqrs0",
"123456789abcdefghijklmnopqrst0",
"123456789abcdefghijklmnopqrstu0",
"123456789abcdefghijklmnopqrstuv0",
"123456789abcdefghijklmnopqrstuvw0",
"123456789abcdefghijklmnopqrstuvwx0",
"123456789abcdefghijklmnopqrstuvwxy0",
"123456789abcdefghijklmnopqrstuvwxyz0",
]

// Radix test
func test9() {
for r in 2 ... 36 {
for i in 0 ..< strings.count {
let x = BInt(strings[i], radix: r)
let X = BInt(strings[i].uppercased(), radix: r)
if i < r {
XCTAssertNotNil(x)
let z = x?.asString(radix: r, uppercase: false)
XCTAssertEqual(z, strings[i])
XCTAssertNotNil(X)
let Z = X?.asString(radix: r, uppercase: true)
XCTAssertEqual(Z, strings[i].uppercased())
XCTAssertEqual(BInt(z!, radix: r), BInt(Z!, radix: r))
} else {
XCTAssertNil(x)
XCTAssertNil(X)
}
}
}
}

}
1 change: 1 addition & 0 deletions Tests/BigIntTests/GcdExtendedTest.swift
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ class GcdExtendedTest: XCTestCase {
doTest(BInt.ONE, BInt.ONE)
for _ in 0 ..< 100 {
let a = BInt(bitWidth: 100)
doTest(a, BInt.ZERO)
doTest(a, BInt.ONE)
doTest(a, BInt.TWO)
for _ in 0 ..< 100 {
Expand Down
Loading

0 comments on commit 2b2f7f4

Please sign in to comment.