Category talk:Wren-math: Difference between revisions
→Source code: Added squareFree, cubeFree and powerFree methods to Int class.
(→Source code: Fixed a problem with Int.modMul, added Int.nextPrime2 method.) |
(→Source code: Added squareFree, cubeFree and powerFree methods to Int class.) |
||
| (15 intermediate revisions by the same user not shown) | |||
Line 1:
===Source code===
<syntaxhighlight lang="
/* Math supplements the Num class with various other operations on numbers. */
Line 113:
return deriv
}
// Attempts to find a real root of a polynomial using Newton's method from
// an initial guess, a given tolerance, a maximum number of iterations
// and a given multiplicity (usually 1).
// If a root is found, it is returned otherwise 'null' is returned.
// If the root is near an integer checks if the integer is in fact the root.
static rootPoly(coefs, guess, tol, maxIter, mult) {
var deg = coefs.count - 1
if (deg == 0) return null
if (deg == 1) return -coefs[1]/coefs[0]
if (deg == 2 && coefs[1]*coefs[1] - 4*coefs[0]*coefs[2] < 0) return null
if (Math.evalPoly(coefs, 0) == 0) return 0
var eps = 0.001
var deriv = Math.diffPoly(coefs)
var x0 = guess
var iter = 1
while (true) {
var den = Math.evalPoly(deriv, x0)
if (den == 0) {
x0 = (x0 >= 0) ? x0 + eps : x0 - eps
} else {
var num = Math.evalPoly(coefs, x0)
var x1 = x0 - num/den * mult
if ((x1 - x0).abs <= tol) {
var r = x1.round
if ((r - x1).abs <= eps && Math.evalPoly(coefs, r) == 0) return r
return x1
}
x0 = x1
}
if (iter == maxIter) break
iter = iter + 1
}
x0 = x0.round
if (Math.evalPoly(coefs, x0) == 0) return x0
return null
}
// Convenience versions of 'rootPoly' which use defaults for some parameters.
static rootpoly(coefs, guess) { rootPoly(coefs, guess, 1e-15, 100, 1) }
static rootPoly(coefs) { rootPoly(coefs, 0.001, 1e-15, 100, 1) }
// Gamma function using Lanczos approximation.
Line 335 ⟶ 377:
static binomial(n, k) { multinomial(n, [k, n-k]) }
//
static maxSafePrime { 9007199254740881 }
// Private helper method for 'isPrime' method.
// Returns true if 'n' is prime, false otherwise but uses the
// Miller-Rabin test which should be faster for 'n' >= 2 ^ 32.
static
if (
if (n%2 == 0) return n == 2
if (n%3 == 0) return n == 3
if (n%5 == 0) return n == 5
if (n%7 == 0) return n == 7
var a
if (n <
a = [2, 7, 61]
} else if (n < 1122004669633) {
Line 423 ⟶ 428:
}
}
}
return true
}
// Determines whether 'n' is prime using a wheel with basis [2, 3, 5].
// Not accessing the wheel via a list improves performance by about 25%.
// Automatically changes to Miller-Rabin approach if 'n' >= 2 ^ 32.
static isPrime(n) {
if (!n.isInteger || n < 2) return false
if (n > 4294967295) return isPrimeMR_(n)
if (n%2 == 0) return n == 2
if (n%3 == 0) return n == 3
if (n%5 == 0) return n == 5
var d = 7
while (d*d <= n) {
if (n%d == 0) return false
d = d + 4
if (n%d == 0) return false
d = d + 2
if (n%d == 0) return false
d = d + 4
if (n%d == 0) return false
d = d + 2
if (n%d == 0) return false
d = d + 4
if (n%d == 0) return false
d = d + 6
if (n%d == 0) return false
d = d + 2
if (n%d == 0) return false
d = d + 6
if (n%d == 0) return false
}
return true
Line 429 ⟶ 466:
// Returns the next prime number greater than 'n'.
static nextPrime(n) {
if (n >= maxSafePrime) Fiber.abort("Argument is larger than maximum safe prime.")
if (n < 2) return 2
n = (n%2 == 0) ? n + 1 : n + 2
Line 437 ⟶ 475:
}
// Returns the previous prime number less than 'n' or null if there isn't one.
static prevPrime(n) {
if (n > Num.maxSafeInteger) Fiber.abort("Argument must be a 'safe' integer.")
if (n <
n = (n%2 == 0) ? n - 1 : n - 2
while (true) {
if (Int.
n = n
}
}
// Returns the 'n'th prime number. For the formula used, see:
// https://en.wikipedia.org/wiki/Prime_number_theorem#Approximations_for_the_nth_prime_number
static getPrime(n) {
if (n < 1 || !n.isInteger) Fiber.abort("Argument must be a positive integer.")
if (n < 8) return [2, 3, 5, 7, 11, 13, 17][n-1]
var l1 = n.log
var l2 = l1.log
var l3 = (l2 - 2) / l1
var l4 = (l2*l2 - 6*l2 + 11) / (2*l1*l1)
var p = ((l1 + l2 + l3 - l4 - 1) * n).floor
var g = p.log.round
var pc = Int.primeCount(p)
p = p + (n - pc) * g
if (p % 2 == 0) p = p - 1
pc = Int.primeCount(p)
if (pc < n) {
for (i in pc...n) p = Int.nextPrime(p)
} else if (Int.isPrime(p)) {
for (i in n...pc) p = Int.prevPrime(p)
} else {
for (i in n..pc) p = Int.prevPrime(p)
}
return p
}
Line 547 ⟶ 611:
var rtlmt = n.sqrt.floor
var mxndx = Int.quo(rtlmt - 1, 2)
var arrlen = (mxndx + 1).
var smalls = List.filled(arrlen, 0)
var roughs = List.filled(arrlen, 0)
Line 617 ⟶ 681:
}
return ans + 1
}
// Returns how many positive integers up to 'n' are relatively prime to 'n'.
static totient(n) {
if (n < 1) return 0
var tot = n
var i = 2
while (i*i <= n) {
if (n%i == 0) {
while (n%i == 0) n = (n/i).floor
tot = tot - (tot/i).floor
}
if (i == 2) i = 1
i = i + 2
}
if (n > 1) tot = tot - (tot/n).floor
return tot
}
Line 622 ⟶ 703:
static primeFactors(n) {
if (!n.isInteger || n < 2) return []
if (n > Num.maxSafeInteger) Fiber.abort("Argument must be a 'safe' integer.")
var inc = [4, 2, 4, 2, 4, 6, 2, 6]
var factors = []
Line 650 ⟶ 732:
return factors
}
// Returns a list of the distinct prime factors of 'n' in order.
static distinctPrimeFactors(n) {
var factors = primeFactors(n)
if (factors.count < 2) return factors
for (i in factors.count-1..1) {
if (factors[i] == factors[i-1]) factors.removeAt(i)
}
return factors
}
// Returns a list of the distinct prime factors of 'n' together with the number
// of times each such factor is repeated.
static primePowers(n) {
var factors = primeFactors(n)
if (factors.count == 0) return []
var prev = factors[0]
var res = [[prev, 1]]
for (f in factors.skip(1)) {
if (f == prev) {
res[-1][1] = res[-1][1] + 1
} else {
res.add([f, 1])
}
prev = f
}
return res
}
// Returns true if the prime factorization of 'n' does not contain any
// factors which are repeated 'm' (or more) times, or false otherwise.
static powerFree(m, n) { Int.primePowers(n).all { |pp| pp[1] < m } }
// Covenience versions of above method for squares and cubes.
static squareFree(n) { powerFree(2, n) }
static cubeFree(n) { powerFree(3, n) }
// Returns all the divisors of 'n' including 1 and 'n' itself.
static divisors(n) {
if (!n.isInteger || n < 1) return []
if (n > Num.maxSafeInteger) Fiber.abort("Argument must be a 'safe' integer.")
var divisors = []
var divisors2 = []
Line 725 ⟶ 844:
static divisorSum(n) {
if (!n.isInteger || n < 1) return 0
if (n > Num.maxSafeInteger) Fiber.abort("Argument must be a 'safe' integer.")
var total = 1
var power = 2
Line 751 ⟶ 871:
static divisorCount(n) {
if (!n.isInteger || n < 1) return 0
if (n > Num.maxSafeInteger) Fiber.abort("Argument must be a 'safe' integer.")
var count = 0
var prod = 1
Line 774 ⟶ 895:
// Private helper method which checks a number and base for validity.
static check_(n, b) {
if (!(n is Num && n.isInteger && n >= 0 && n <= Num.maxSafeInteger)) {
Fiber.abort("Number must be a non-negative 'safe' integer.")
}
if (!(b is Num && b.isInteger && b >= 2 && b < 64)) {
Line 816 ⟶ 937:
return [n, ap]
}
// Returns a new integer formed by reversing the base 10 digits of 'n'.
// Works with positive, negative and zero integers.
static reverse(n, check) {
if (check) check_(n, b)
var m = (n >= 0) ? n : -n
var r = 0
Line 837 ⟶ 950:
return r * n.sign
}
// Convenience versions of the above methods which never check for validity
// and/or use base 10 by default.
static digits(n, b) { digits(n, b, false) }
static digits(n) { digits(n, 10, false) }
static digitSum(n, b) { digitSum(n, b, false) }
static digitSum(n) { digitSum(n, 10, false) }
static digitalRoot(n, b) { digitalRoot(n, b, false) }
static digitalRoot(n) { digitalRoot(n, 10, false) }
static reverse(n) { reverse(n, false) }
// Returns the unique non-negative integer that is associated with a pair
Line 950 ⟶ 1,073:
// Converts a numeric list to a corresponding string list.
static toStrings(a) { a.map { |n| n.toString }.toList }
// Draws a horizontal bar chart on the terminal representing a list of numerical 'data'
// which must be non-negative. 'labels' is a list of the corresponding text for each bar.
// When printed and unless they are all the same length, 'labels' are right justified
// within their maximum length if 'rjust' is true, otherwise they are left justified.
// 'width' is the total width of the chart in characters to which the labels/data will
// automatically be scaled. 'symbol' is the character to be used to draw each bar,
// 'symbol2' is used to represent scaled non-zero data and 'symbol3' zero data which
// would otherwise be blank. The actual data is printed at the end of each bar.
static barChart (title, width, labels, data, rjust, symbol, symbol2, symbol3) {
var barCount = labels.count
if (data.count != barCount) Fiber.abort("Mismatch between labels and data.")
var maxLabelLen = max(labels.map { |l| l.count })
if (labels.any { |l| l.count < maxLabelLen }) {
if (rjust) {
labels = labels.map { |l| (" " * (maxLabelLen - l.count)) + l }.toList
} else {
labels = labels.map { |l| l + (" " * (maxLabelLen - l.count)) }.toList
}
}
var maxData = max(data)
var maxLen = width - maxLabelLen - maxData.toString.count - 2
var scaledData = data.map { |d| (d * maxLen / maxData).floor }.toList
System.print(title)
System.print("-" * Math.max(width, title.count))
for (i in 0...barCount) {
var bar = symbol * scaledData[i]
if (bar == "") bar = data[i] > 0 ? symbol2 : symbol3
System.print(labels[i] + " " + bar + " " + data[i].toString)
}
System.print("-" * Math.max(width, title.count))
}
// Convenience version of the above method which uses default symbols.
static barChart(title, width, labels, data, rjust) {
barChart(title, width, labels, data, rjust, "■", "◧", "□")
}
// As above method but uses right justification for labels.
static barChart(title, width, labels, data) {
barChart(title, width, labels, data, true, "■", "◧", "□")
}
// Plots a sequence of points 'pts' on x/y axes of lengths 'lenX'/'lenY' on the
// terminal automatically scaling them to those lengths. 'symbol' marks each point.
static plot(title, pts, lenX, lenY, symbol) {
var px = pts.map { |p| p[0] }
var py = pts.map { |p| p[1] }
var maxX = max(px).ceil
var minX = min(px).floor
var maxY = max(py).ceil
var minY = min(py).floor
var cy = max([maxY.toString.count, minY.toString.count, 5])
pts = pts.map { |p|
var q = [0, 0]
q[0] = Math.lerp(1, lenX, (p[0] - minX) / (maxX - minX)).round
q[1] = Math.lerp(1, lenY, (p[1] - minY) / (maxY - minY)).round
return q
}.toList
pts.sort { |p, q| p[1] != q[1] ? p[1] > q[1] : p[0] < q[0] }
System.print(title)
System.print("-" * (lenX + cy + 3))
var pc = 0
for (line in lenY..1) {
var s = ""
if (line == lenY) {
s = maxY.toString
} else if (line == 1) {
s = minY.toString
}
System.write("|")
if (s.count < cy) System.write(" " * (cy - s.count))
System.write(s + "|")
var xs = []
while (pc < pts.count && pts[pc][1] == line) {
xs.add(pts[pc][0])
pc = pc + 1
}
if (xs.isEmpty) {
System.print(" " * lenX + "|")
} else {
var lastX = 0
for (x in xs) {
if (x == lastX) continue
if (x - lastX > 1) System.write(" " * (x - lastX - 1))
System.write(symbol)
lastX = x
}
if (lenX > lastX) System.write(" " * (lenX - lastX))
System.print("|")
}
}
System.print("|" + "-" * (lenX + cy + 1) + "|")
var minL = minX.toString.count
var maxL = maxX.toString.count
System.write("| y/x" + " " * (cy - 4) + "|")
System.print(minX.toString + " " * (lenX - minL - maxL) + maxX.toString + "|")
System.print("-" * (lenX + cy + 3))
}
// As above method but uses a default symbol.
static plot(title, pts, lenX, lenY) { plot(title, pts, lenX, lenY, "▮") }
}
| |||