Rosetta Code

Category talk:Wren-i64: Difference between revisions

Category Discussion

Category talk:Wren-i64 (view source)

Revision as of 20:24, 9 October 2022

7,071 bytes added ,  1 year ago
→‎Source code (Wren): Fixed some bugs and added prime factorization etc. stuff.
m (Added quotes to 'lang' attribute.)
(→‎Source code (Wren): Fixed some bugs and added prime factorization etc. stuff.)
Line 427:
}
return s
}
 
// 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 U64)) x = from(x)
if (!(n is U64)) n = from(n) else n = n.copy()
if (!(mod is U64)) mod = from(mod)
var z = zero
var y = x.copy()
while (n > zero) {
if (n.isOdd) z.add(y).rem(mod)
y.lsh(1).rem(mod)
n.rsh(1)
}
return z
}
 
Line 433 ⟶ 449:
static modPow(x, exp, mod) {
if (!(x is U64)) x = from(x)
if (!(exp is U64)) exp = from(exp) else exp = exp.copy()
if (!(mod is U64)) mod = from(mod)
if (mod.isZero) Fiber.abort("Cannot take modPow with modulus 0.")
var r = one
var base = x % mod
while (!exp > 0.isZero) {
if (base.isZero) return zero
if (exp.isOdd) r.set(modmulmodMul(r, base, mod))
exp.rsh(1)
base.set(modMul(base, base, mod))
}
return r
}
 
// 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 U64)) x = from(x)
if (!(n is U64)) n = from(n)
if (!(mod is U64)) mod = from(mod)
var z = zero
var y = x.copy()
while (n > 0) {
if (n.isOdd) x.set((x + y) % mod)
y.lsh(1).rem(mod)
n.rsh(1)
}
return x
}
 
Line 497:
}
 
// Returns the factorial of 'n'.
// Returns the factorial of 'n'.
static factorial(n) {
Line 533 ⟶ 534:
static binomial(n, k) { multinomial(n, [k, n-k]) }
 
// DeterminesReturns whether or not 'x' is prime using a wheel withan basisinstance [2,of 3]U64.
static isPrimeisInstance(x) { x is U64 }
 
if (!(x is U64)) x = from(x)
// Private method to determine if (xa U64 is a <basic 2)prime returnor falsenot.
static isBasicPrime_(n) {
if (x%2 == 0) return x == 2
if (x%3 == 0n.isOne) return x == 3false
varif d(n == two || n == three || n == U64.five) return true
if (n.isEven || n.isDivisible(three) || n.isDivisible(five)) return false
while (d*d <= x) {
if (x%dn ==< 049) return falsetrue
return null // unknown d.add(2)if prime or not
}
if (x%d == 0) return false
 
d.add(4)
// 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.square.rem(n)
if (x.isOne) return false
if (x == nPrev) {
next = true
break
}
d = d - 1
}
if (!next) return false
}
}
}
return true
}
 
// Returns true if 'x' is prime, false otherwise.
// Fails due to overflow on numbers > 9223372036854775807 unless divisible by 2,3 or 5.
static isPrime(x) {
if (!(x is U64)) x = from(x)
var isbp = isBasicPrime_(x)
if (isbp != null) return isbp
if (x > largest/two) Fiber.abort("Cannot test %(x) for primality.")
return millerRabinTest_(x, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37])
}
 
Line 552 ⟶ 591:
static nextPrime(x) {
if (!(x is U64)) x = from(x)
xvar n = x.isEven ? x + one : x + two
while (true) {
if (isPrime(xn)) return xn
x = x + n.add(two)
}
}
 
// ReturnsPrivate whetherworker ormethod notfor Pollard'x' is an instances ofRho U64algorithm.
static isInstancepollardRho_(xm, seed, c) { x is U64 }
var g = Fn.new { |x| x.copy().square.add(c).rem(m) }
var x = from(seed)
var y = from(seed)
var z = one
var d = one
var count = 0
while (true) {
x = g.call(x)
y = 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)
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 large numbers with a small number of prime factors.
static divisors2(n) {
if (n <= 0) return []
var factors = primeFactors(n)
var divs = [one]
for (p in factors) {
for (i in 0...divs.count) divs.add(divs[i]*p)
}
divs.sort()
var c = divs.count
if (c > 1) {
for (i in c-1..1) {
if (divs[i-1] == divs[i]) divs.removeAt(i)
}
}
return divs
}
 
// As 'properDivisors' but uses 'divisors2' method.
static properDivisors2(n) {
var d = divisors2(n)
var c = d.count
return (c <= 1) ? [] : d[0..-2]
}
 
// Creates a new U64 object with a value of zero.
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/331492"