Category talk:Wren-i64: Difference between revisions
→Source code (Wren): Bug fix
(→Source code (Wren): Added U64.nextPrime method.) |
(→Source code (Wren): Bug fix) |
||
| (17 intermediate revisions by the same user not shown) | |||
Line 1:
==64-bit integer arithmetic==
Although this is already supported in an easy to use way by the Wren-long module, there are some RC tasks which that module (written entirely in Wren) struggles to complete in an acceptable time or to a satisfactory standard
I have therefore written a wrapper for the C99 fixed length types, int64_t and
These types are designed to behave the same as the C types regards rounding, underflow and overflow.
Line 11:
As well as methods for basic arithmetic, I have also included some additional 'convenience' or utility methods which can be coded in a few lines of Wren code using the former.
This module is written to be as fast as possible within the confines of the Wren VM and therefore follows the model of Wren-GMP in recycling memory to minimize the number of variables and garbage collections needed in a typical application. This means that there is an inevitable trade-off between maximizing efficiency/performance and easy of use which is an important consideration for a simple language such as Wren. In an attempt to retain as much of the latter as possible:
1. I have included versions of the main functions which operate on and always return a reference to the current object which reduces verbosity and enables method chaining. These methods also sense the type of the argument eliminating the need to have separate versions for Wren-i64 types and built-in numbers (Num).
2. The two classes: I64 and U64 both inherit from the Comparable trait which
3. I have also included operator overloading for the basic arithmetic operations with the difference that these always produce a new object rather than mutating the current one. Consequently, these should be used sparingly in scripts which need to maximize performance but can be used more freely in other scripts to provide a more 'natural' programming experience.
Line 29 ⟶ 28:
==Source code (Wren)==
<
import "./trait" for Comparable
Line 48 ⟶ 47:
static four { from(4) }
static five { from(5) }
static six { from(6) }
static seven { from(7) }
static eight { from(8) }
static nine { from(9) }
static ten { from(10) }
Line 59 ⟶ 62:
// The arguments can be either I64 objects or Nums.
// If 'n' is less than 0, returns zero. O.pow(0) returns one.
static pow(x, n) {
if
if
if (n < 0) return zero
if (n.isZero) return one
Line 73 ⟶ 75:
value = value.pow(n.toNum)
if (value <= 9007199254740991) return neg ? from(-value) : from(value)
x = x.copy()
var y = one
var z = n.copy()
while (true) {
if (z.isOdd) {
Line 84 ⟶ 85:
if (z.isZero) break
z.rsh(1)
x.
}
return y
Line 94 ⟶ 95:
// Throws an error if 'x' is negative and 'n' is even.
static root(x, n) {
if
if (!((n is Num) && n.isInteger && n > 0)) {
Fiber.abort("'n' must be a positive integer.")
Line 102 ⟶ 103:
var neg = t < 0
if (neg) {
if (n
Fiber.abort("Cannot take the %(n)th root of a negative number.")
} else {
Line 116 ⟶ 117:
}
return (neg) ? -s : s
}
// Returns an I64 object equal in value to 'x' multiplied by 'n' modulo 'mod'.
// The arguments can either be I64 objects or Nums but 'mod' must be non-zero.
static modMul(x, n, mod) {
x = from(x)
n = from(n)
mod = from(mod)
if (mod.isZero) Fiber.abort("Cannot take modMul with modulus 0.")
var sign = x.sign * n.sign
var mag = U64.modMul(x.abs, n.abs, mod.abs).toI64
return mag * sign
}
Line 121 ⟶ 134:
// The arguments can either be I64 objects or Nums but 'mod' must be non-zero.
static modPow(x, exp, mod) {
if
if
if (mod.isZero) Fiber.abort("Cannot take modPow with modulus 0.")
var r = one
Line 133 ⟶ 146:
while (exp > 0) {
if (base.isZero) return zero
if (exp.isOdd) r.set(modMul(r
exp.rsh(1)
base.
}
return r
Line 143 ⟶ 156:
// The arguments can either be I64 objects or Nums but must be co-prime.
static modInv(x, n) {
if
if
var r = from(n)
var newR = from(x).abs
Line 166 ⟶ 179:
// Returns the smaller of two I64 objects.
static min(x, y) {
if (x is Num) x = from(x)
if (y is Num) y = from(y)
if (x < y) return x
return y
Line 172 ⟶ 187:
// Returns the greater of two I64 objects.
static max(x, y) {
if (x is Num) x = from(x)
if (y is Num) y = from(y)
if (x > y) return x
return y
Line 178 ⟶ 195:
// Returns the positive difference of two I64 objects.
static dim(x, y) {
if (x is Num) x = from(x)
if (y is Num) y = from(y)
if (x >= y) return x - y
return zero
Line 184 ⟶ 203:
// Returns the greatest common divisor of two I64 objects.
static gcd(x, y) {
x = from(x)
y = from(y)
while (y != zero) {
var t = y
y
x
}
return x.abs
Line 194 ⟶ 215:
// Returns the least common multiple of two I64 objects.
static lcm(x, y) {
x = from(x)
y = from(y)
if (x.isZero && y.isZero) return zero
return (x*y).abs / gcd(x, y)
Line 323 ⟶ 346:
foreign cmpU64(op) // compare this to a U64 object
/* Miscellaneous methods. */
Line 329 ⟶ 352:
max(op) { (op > this) ? op : this } // returns the maximum of this and another I64 object
clamp(min, max) {
if (
if (this < min) set(min) else if (this > max) set(max)
return this.copy()
}
copy() { I64.from(this) }
isOdd { this % I64.two
isEven { this % I64.two == I64.zero } // true if 'this' is even
isZero { this == I64.zero } // true if 'this' is zero
Line 342 ⟶ 366:
isSafe { this >= I64.minSafe && this <= I64.maxSafe } // true if 'this' is safe
isDivisible(d) { this % d == I64.zero }
isRoot(n) {
var r = I64.root(this, n)
return I64.pow(r, n) == this
}
isSquare { isRoot(2) }
isCube { isRoot(3) }
sign {
if (this < I64.zero) return -1
if (this > I64.zero) return 1
Line 373 ⟶ 397:
static four { from(4) }
static five { from(5) }
static six { from(6) }
static seven { from(7) }
static eight { from(8) }
static nine { from(9) }
static ten { from(10) }
Line 379 ⟶ 407:
static maxSafe { from(9007199254740991) }
static maxU32 { from(4294967295) }
// Returns the largest prime less than 2^64.
static largestPrime { U64.fromStr("18446744073709551557") }
// Returns a U64 object equal in value to 'x' raised to the power of 'n'.
// The arguments can be either U64 objects or Nums.
// If 'n' is less than 0, returns zero. O.pow(0) returns one.
static pow(x, n) {
if
if
if (n.isZero) return one
if (n.isOne) return x.copy()
Line 393 ⟶ 423:
var value = x.toNum.pow(n.toNum)
if (value <= 9007199254740991) return from(value)
x = x.copy()
var y = one
var z = n.copy()
while (true) {
if (z.isOdd) {
Line 404 ⟶ 433:
if (z.isZero) break
z.rsh(1)
x.
}
return y
Line 413 ⟶ 442:
// 'x' can either be a U64 object or a Num but 'n' must be a positive integer.
static root(x, n) {
if
if (!((n is Num) && n.isInteger && n > 0)) {
Fiber.abort("'n' must be a positive integer.")
Line 427 ⟶ 456:
}
return s
}
// Private helper method for modMul method to avoid overflow.
static checkedAdd_(a, b, c) {
var room = c - U64.one - a
if (b <= room) {
a.add(b)
} else {
a.set(b - room - U64.one)
}
}
// Returns a U64 object equal in value to 'x' multiplied by 'n' modulo 'mod'.
// The arguments can either be U64 objects or Nums.
static modMul(x, n, mod) {
if (x is Num) x = from(x)
n = from(n)
mod = from(mod)
if (mod.isZero) Fiber.abort("Cannot take modMul with modulus 0.")
var z = zero
var y = x % mod
n.rem(mod)
if (n > y) {
var t = y
y = n
n = t
}
while (n > zero) {
if (n.isOdd) checkedAdd_(z, y, mod)
checkedAdd_(y, y, mod)
n.rsh(1)
}
return z
}
Line 432 ⟶ 494:
// The arguments can either be U64 objects or Nums but 'mod' must be non-zero.
static modPow(x, exp, mod) {
if
if
if (mod.isZero) Fiber.abort("Cannot take modPow with modulus 0.")
var r = one
var base = x % mod
while (!exp
if (base.isZero) return zero
if (exp.isOdd) r.set(
exp.rsh(1)
base.set(modMul(base, base, mod))
Line 447 ⟶ 509:
}
// Returns
// The arguments can either be U64 objects or Nums but must be co-prime.
static
if
if
return I64.modInv(x.toI64, n.toI64).toU64
}
// Returns the smaller of two U64 objects.
static min(x, y) {
if (x is Num) x = from(x)
if (y is Num) y = from(y)
if (x < y) return x
return y
Line 471 ⟶ 527:
// Returns the greater of two U64 objects.
static max(x, y) {
if (x is Num) x = from(x)
if (y is Num) y = from(y)
if (x > y) return x
return y
Line 477 ⟶ 535:
// Returns the positive difference of two U64 objects.
static dim(x, y) {
if (x is Num) x = from(x)
if (y is Num) y = from(y)
if (x >= y) return x - y
return zero
Line 483 ⟶ 543:
// Returns the greatest common divisor of two U64 objects.
static gcd(x, y) {
x = from(x)
y = from(y)
while (y != zero) {
var t = y
y
x
}
return x
Line 493 ⟶ 555:
// Returns the least common multiple of two U64 objects.
static lcm(x, y) {
x = from(x)
y = from(y)
if (x.isZero && y.isZero) return zero
return x * y / gcd(x, y)
}
// Returns the factorial of 'n'.
// Returns the factorial of 'n'.
static factorial(n) {
Line 504 ⟶ 569:
if (n < 2) return one
var fact = one
var i =
while (i <= n) {
fact.mul(i)
i
}
return fact
Line 527 ⟶ 592:
if (e > 1) prod.mul(factorial(e))
}
return factorial(n)
}
Line 533 ⟶ 598:
static binomial(n, k) { multinomial(n, [k, n-k]) }
//
static
// Private method to apply the Miller-Rabin test.
static millerRabinTest_(n, a) {
var nPrev = n.copy().dec
var b = nPrev.copy()
var r = 0
while (b.isEven) {
b.rsh(1)
r = r + 1
}
for (i in 0...a.count) {
if (n >= a[i]) {
var x = (a[i] is U64) ? a[i].copy() : from(a[i])
x.set(modPow(x, b, n))
if (!x.isOne && x != nPrev) {
var d = r - 1
var next = false
while (d != 0) {
x.set(modMul(x, x, n))
if (x.isOne) return false
if (x == nPrev) {
next = true
break
}
d = d - 1
}
if (!next) return false
}
}
}
return true
}
// Private helper method for 'isPrime' method.
// Determines whether 'x' is prime using a wheel with basis [2, 3, 5].
// Should be faster than Miller-Rabin if 'x' is below 2^37.
static isPrimeWheel_(x) {
if (x < 2) return false
if (
var
while (d*d <= n) {
if ((n%d).isZero) return false
d.add(4)
if ((n%d).isZero) return false
d.add(2)
if (
d.add(4)
if ((n%d).isZero) return false
d.add(2)
if ((n%d).isZero) return false
d.add(4)
if ((n%d).isZero) return false
d.add(6)
if ((n%d).isZero) return false
d.add(2)
if ((n%d).isZero) return false
d.add(6)
if ((n%d).isZero) return false
}
return true
}
// Returns true if 'x' is prime, false otherwise.
static isPrime(x) {
if (x is Num) x = from(x)
if (x < 137438953472) return isPrimeWheel_(x) // use wheel if below 2^37
if (x.isEven || x.isDivisible(three) || x.isDivisible(five) ||
x.isDivisible(seven)) return false
var a
if (x < 1122004669633) {
a = [2, 13, 23, 1662803]
} else if (x < 2152302898747) {
a = [2, 3, 5, 7, 11]
} else if (x < 3474749660383) {
a = [2, 3, 5, 7, 11, 13]
} else if (x < 341550071728321) {
a = [2, 3, 5, 7, 11, 13, 17]
} else if (x < fromStr("3825123056546413051")) {
a = [2, 3, 5, 7, 11, 13, 17, 19, 23]
} else {
a = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
}
return millerRabinTest_(x, a)
}
// Returns the next prime number greater than 'x'.
static nextPrime(x) {
if
var n = x.isEven ? x + one : x + two
while (true) {
if (isPrime(
}
}
// Returns
static
if (x is Num) x = from(x)
if (x < three) return null
if (x == three) return U64.two
var n = x.isEven ? x - one : x - two
while (true) {
if (isPrime(n)) return n
n.sub(two)
}
}
// Private worker method for Pollard's Rho algorithm.
static pollardRho_(m, seed, c) {
var g = Fn.new { |x| (x * x).add(c).rem(m) }
var x = from(seed)
var y = from(seed)
var z = one
var d = one
var count = 0
while (true) {
x.set(g.call(x))
y.set(g.call(g.call(y)))
if (x >= y) d.sub(x, y).rem(m) else d.sub(y, x).rem(m)
z.mul(d)
count = count + 1
if (count == 100) {
d.set(gcd(z, m))
if (d != one) break
z.set(one)
count = 0
}
}
if (d == m) return zero
return d
}
// Returns a factor (U64) of 'm' (a U64 object or an integral Num) using the
// Pollard's Rho algorithm. Both the 'seed' and 'c' can be set to integers.
// Returns U64.zero in the event of failure.
static pollardRho(m, seed, c) {
if (m < 2) return U64.zero
if (m is Num) m = from(m)
if (isPrime(m)) return m.copy()
if (m.isSquare) return m.copy().sqrt
for (p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]) {
if (m.isDivisible(p)) return from(p)
}
return pollardRho_(m, seed, c)
}
// Convenience version of the above method which uses a seed of 2 and a value for c of 1.
static pollardRho(m) { pollardRho(m, 2, 1) }
// Private method for factorizing smaller numbers (U64) using a wheel with basis [2, 3, 5].
static primeFactorsWheel_(m) {
var n = m.copy()
var inc = [4, 2, 4, 2, 4, 6, 2, 6]
var factors = []
var k = from(37)
var i = 0
while (k * k <= n) {
if (n.isDivisible(k)) {
factors.add(k.copy())
n.div(k)
} else {
k.add(inc[i])
i = (i + 1) % 8
}
}
if (n > one) factors.add(n)
return factors
}
// Private worker method (recursive) to obtain the prime factors of a number (U64).
static primeFactors_(m, trialDivs) {
if (isPrime(m)) return [m.copy()]
var n = m.copy()
var factors = []
var seed = 2
var c = 1
var checkPrime = true
var threshold = 1e11 // from which using PR may be advantageous
if (trialDivs) {
for (p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]) {
while (n.isDivisible(p)) {
factors.add(from(p))
n.div(p)
}
}
}
while (n > one) {
if (checkPrime && isPrime(n)) {
factors.add(n.copy())
break
}
if (n >= threshold) {
var d = pollardRho_(n, seed, c)
if (d != zero) {
factors.addAll(primeFactors_(d, false))
n.div(d)
checkPrime = true
} else if (c == 1) {
if (n.isSquare) {
n.sqrt
var pf = primeFactors_(n, false)
factors.addAll(pf)
factors.addAll(pf)
break
} else {
c = 2
checkPrime = false
}
} else if (c < 101) {
c = c + 1
} else if (seed < 101) {
seed = seed + 1
} else {
factors.addAll(primeFactorsWheel_(n))
break
}
} else {
factors.addAll(primeFactorsWheel_(n))
break
}
}
factors.sort()
return factors
}
// Returns a list of the primes factors (U64) of 'm' (a U64 object or an integral Num)
// using the wheel based factorization and/or Pollard's Rho algorithm as appropriate.
static primeFactors(m) {
if (m < 2) return []
if (m is Num) m = from(m)
return primeFactors_(m, true)
}
// Returns all the divisors of 'n' including 1 and 'n' itself.
static divisors(n) {
if (n < 1) return []
if (n is Num) n = from(n)
var divs = []
var divs2 = []
var i = one
var k = n.isEven ? one : two
var sqrt = n.copy().sqrt
while (i <= sqrt) {
if (n.isDivisible(i)) {
divs.add(i.copy())
var j = n / i
if (j != i) divs2.add(j)
}
i.add(k)
}
if (!divs2.isEmpty) divs = divs + divs2[-1..0]
return divs
}
// Returns all the divisors of 'n' excluding 'n'.
static properDivisors(n) {
var d = divisors(n)
var c = d.count
return (c <= 1) ? [] : d[0..-2]
}
// As 'divisors' method but uses a different algorithm.
// Better for larger numbers.
static divisors2(n) {
if (n is Num) n = from(n)
var pf = primeFactors(n)
if (pf.count == 0) return (n == 1) ? [one] : pf
var arr = []
if (pf.count == 1) {
arr.add([pf[0].copy(), 1])
} else {
var prevPrime = pf[0]
var count = 1
for (i in 1...pf.count) {
if (pf[i] == prevPrime) {
count = count + 1
} else {
arr.add([prevPrime.copy(), count])
prevPrime = pf[i]
count = 1
}
}
arr.add([prevPrime.copy(), count])
}
var divisors = []
var generateDivs
generateDivs = Fn.new { |currIndex, currDivisor|
if (currIndex == arr.count) {
divisors.add(currDivisor.copy())
return
}
for (i in 0..arr[currIndex][1]) {
generateDivs.call(currIndex+1, currDivisor)
currDivisor = currDivisor * arr[currIndex][0]
}
}
generateDivs.call(0, one)
return divisors.sort()
}
// As 'properDivisors' but uses 'divisors2' method.
static properDivisors2(n) {
var d = divisors2(n)
var c = d.count
return (c <= 1) ? [] : d[0..-2]
}
// Returns the sum of all the divisors of 'n' including 1 and 'n' itself.
static divisorSum(n) {
if (n < 1) return zero
n = from(n)
var total = one
var power = two
while (n.isDivisible(two)) {
total.add(power)
power.mul(2)
n.div(2)
}
var i = three
while (i * i <= n) {
var sum = one
power.set(i)
while (n.isDivisible(i)) {
sum.add(power)
power.mul(i)
n.div(i)
}
total.mul(sum)
i.add(2)
}
if (n > 1) total.mul(n + 1)
return total
}
// Returns the number of divisors of 'n' including 1 and 'n' itself.
static divisorCount(n) {
if (n < 1) return 0
n = from(n)
var count = 0
var prod = 1
while (n.isDivisible(two)) {
count = count + 1
n.div(2)
}
prod = prod * (1 + count)
var i = three
while (i * i <= n) {
count = 0
while (n.isDivisible(i)) {
count = count + 1
n.div(i)
}
prod = prod * (1 + count)
i.add(2)
}
if (n > 2) prod = prod * 2
return prod
}
// Creates a new U64 object with a value of zero.
Line 683 ⟶ 1,084:
max(op) { (op > this) ? op : this } // returns the maximum of this and another U64 object
clamp(min, max) {
if (
if (this < min) set(min) else if (this > max) set(max)
return this.copy()
}
copy() { U64.from(this) }
isOdd { this % U64.two == 1 } // true if 'this' is odd
isEven { this % U64.two == 0 } // true if 'this' is even
isZero { this == 0 } // true if 'this' is zero
isOne { this == 1 } // true if 'this' is one
isSafe { this <= U64.maxSafe
isDivisible(d) { this % d == U64.zero } // true if 'this' is divisible by d
isRoot(n) {
var r = U64.root(this, n)
return U64.pow(r, n) == this
}
isSquare { isRoot(2) }
isCube { isRoot(3) }
sign {
if (this > U64.zero) return 1
return 0
}
// Return a list of the current instance's base 10 digits
digits { toString.map { |d| Num.fromString(d) }.toList }
digitSum {
var sum = 0
for (d in toString.bytes) sum = sum + d - 48
return sum
}
}
Line 729 ⟶ 1,137:
static min(sz) { sz.reduce { |acc, x| (x > acc) ? x : acc } }
static max(sz) { sz.reduce { |acc, x| (x < acc) ? x : acc } }
}</
==Source code (C)==
<
#include <stdio.h>
Line 1,643 ⟶ 2,051:
free(script);
return 0;
}</
| |||