Category talk:Wren-math

Source code

/* Module "math.wren" */

/* Math supplements the Num class with various other operations on numbers. */
class Math {
    // Constants.
    static e    { 2.71828182845904523536 } // base of natural logarithms
    static phi  { 1.6180339887498948482  } // golden ratio
    static ln2  { 0.69314718055994530942 } // natural logarithm of 2
    static ln10 { 2.30258509299404568402 } // natural logarithm of 10

    // Log function.
    static log10(x) { x.log/ln10 }  // Base 10 logarithm

    // Hyperbolic trig functions.
    static sinh(x) { (x.exp - (-x).exp)/2 } // sine
    static cosh(x) { (x.exp + (-x).exp)/2 } // cosine
    static tanh(x) { sinh(x)/cosh(x)      } // tangent

    // Inverse hyperbolic trig functions.
    static asinh(x) { (x + (x*x + 1).sqrt).log } // sine
    static acosh(x) { (x + (x*x - 1).sqrt).log } // cosine
    static atanh(x) { ((1+x)/(1-x)).log/2 }      // tangent

    // Angle conversions.
    static radians(d) { d * Num.pi / 180}
    static degrees(r) { r * 180 / Num.pi }

    // Returns the square root of 'x' squared + 'y' squared.
    static hypot(x, y) { (x*x + y*y).sqrt }

    // Returns the integer and fractional parts of 'x'. Both values have the same sign as 'x'.
    static modf(x) { [x.truncate, x.fraction] }

    // Returns the IEEE 754 floating-point remainder of 'x/y'.
    static rem(x, y) {
        if (x.isNan || y.isNan || x.isInfinity || y == 0) return nan
        if (!x.isInfinity && y.isInfinity) return x
        var nf = modf(x/y)
        if (nf[1] != 0.5) {
            return x - (x/y).round * y
        } else {
            var n = nf[0]
            if (n%2 == 1) n = (n > 0) ? n + 1 : n - 1
            return x - n * y
        }
    }

    // Round away from zero.
    static roundUp(x) { (x >= 0) ? x.ceil : x.floor }

    // Round to 'p' decimal places, maximum 14.
    // Mode parameter specifies the rounding mode:
    // < 0 towards zero, == 0 nearest, > 0 away from zero.
    static toPlaces(x, p, mode) {
        if (p < 0) p = 0
        if (p > 14) p = 14
        var pw = 10.pow(p)
        var nf = modf(x)
        x = nf[1] * pw
        x = (mode < 0) ? x.truncate : (mode == 0) ? x.round : roundUp(x)
        return nf[0] + x/pw
    }

    // Convenience version of above method which uses 0 for the 'mode' parameter.
    static toPlaces(x, p) { toPlaces(x, p, 0) }

    // Returns the linear interpolation of t between x and y, if t is in [0, 1] or
    // linear extrapolation otherwise.
    static lerp (x, y, t) { x + t * (y - x) }

    // The power of 'p' which equals or first exceeds a non-negative number 'x'.
    // 'p' must be an integer greater than 1.
    // Returns a two element list containing the power and the exponent.
    static nextPower(x, p) {
        if (!((x is Num) && x >= 0)) Fiber.abort("x must be a non-negative number.")
        if (!((p is Num) && p.isInteger && p > 1)) Fiber.abort("p must be an integer > 1.")
        if (x <= 1) return [1, 0]
        var pow = 1
        var count = 0
        while (x > pow) {
            pow = pow * p
            count = count + 1
        }
        return [pow, count]
    }

    // Returns the value of a polynomial when the unknown is 'x' using Horner's rule.
    // Polynomials are represented by an ordered list of coefficients
    // from the term with the highest degree down to the constant term.
    static evalPoly(coefs, x) {
        if (!(coefs is List) || coefs.count == 0 || !(coefs[0] is Num)) {
           Fiber.abort("coefs must be a non-empty ordered list of numbers.")
        }
        if (!(x is Num)) Fiber.abort("x must be a number.")
        var c = coefs.count
        if (c == 1) return coefs[0]
        var sum = coefs[0]
        for (i in 1...c) sum = sum * x + coefs[i]
        return sum
    }

    // Returns the coefficients of a polynomial following differentiation.
    // Polynomials are represented as described in 'evalPoly' above.
    static diffPoly(coefs) {
        if (!(coefs is List) || coefs.count == 0 || !(coefs[0] is Num)) {
           Fiber.abort("coefs must be a non-empty ordered list of numbers.")
        }
        var c = coefs.count - 1
        if (c == 0) return [0]
        var deriv = List.filled(c, 0)
        for (i in 0...c) deriv[i] = (c - i) * coefs[i]
        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.
    static gamma(x) {
        var p = [
            0.99999999999980993,
          676.5203681218851,
        -1259.1392167224028,
          771.32342877765313,
         -176.61502916214059,
           12.507343278686905,
           -0.13857109526572012,
            9.9843695780195716e-6,
            1.5056327351493116e-7
        ]
        var t = x + 6.5
        var sum = p[0]
        for (i in 0..7) sum = sum + p[i+1]/(x + i)
        return 2.sqrt * Num.pi.sqrt * t.pow(x-0.5) * (-t).exp * sum
    }

    // The natural logarithm of the absolute value of gamma(x).
    static lgamma(x) {
        if (x == 1 || x == 2) return 0
        return Math.gamma(x).abs.log
    }

    // Returns the positive difference of two numbers.
    static dim(x, y) { (x >= y) ? x - y : 0 }

    // Static alternatives to instance methods in Num class.
    // Clearer when both arguments are complex expressions.    
    static min(x, y)  { (x < y) ? x : y }
    static max(x, y)  { (x > y) ? x : y }
    static atan(x, y) { x.atan(y) }
}

/* Int contains various routines which are only applicable to integers. */
class Int {
    // Truncated integer division (consistent with % operator).
    static quo(x, y) { (x/y).truncate }

    // Floored integer division (consistent with 'mod' method below).
    static div(x, y) { (x/y).floor }

    // Floored integer division modulus.
    static mod(x, y) { ((x % y) + y) % y }

    // Returns the integer square root of 'x' or null if 'x' is negative.
    static sqrt(x) { (x >= 0) ? x.sqrt.floor : null }

    // Returns the integer cube root of 'x'.
    static cbrt(x) { x.cbrt.truncate }

    // Returns the integer 'n'th root of 'x' or null if 'x' is negative and 'n' is even.
    static root(n, x) {
        if (!(n is Num) || !n.isInteger || n < 1) {
            Fiber.abort("n must be a positive integer.")
        }
        return (n == 1) ? x :
               (n == 2) ? sqrt(x) :
               (n == 3) ? cbrt(x) :
               (n % 2 == 1) ? x.sign * x.abs.pow(1/n).floor :
               (n >= 0) ? x.pow(1/n).floor : null
    }

    // Returns whether or not 'x' is a perfect square.
    static isSquare(x) {
        var s = sqrt(x)
        return s * s == x
    }

    // Returns whether or not 'x' is a perfect cube.
    static isCube(x) {
        var c = cbrt(x)
        return c * c * c == x
    }

    // Returns whether or not 'x' is a perfect 'n'th power.
    static isPower(n, x) {
        var r = root(n, x)
        return r.pow(n) == x
    }

    // Returns the greatest common divisor of 'x' and 'y'.
    static gcd(x, y) {
        while (y != 0) {
            var t = y
            y = x % y
            x = t
        }
        return x.abs
    }

    // Returns the least common multiple of 'x' and 'y'.
    static lcm(x, y) { 
        if (x == 0 && y == 0) return 0
        return (x*y).abs / gcd(x, y) 
    }

    // Returns the greatest common divisor of a list of integers 'a'.
    static gcd(a) {
        if (!(a is List) || a.count < 2) {
            Fiber.abort("Argument must be a list of at least two integers.")
        }
        var g = gcd(a[0], a[1])
        if (a.count == 2) return g
        return gcd(a[2..-1] + [g])
    }

    // Returns the least common multiple of a list of integers 'a'.
    static lcm(a) {
        if (!(a is List) || a.count < 2) {
            Fiber.abort("Argument must be a list of at least two integers.")
        }
        var l = lcm(a[0], a[1])
        if (a.count == 2) return l
        return lcm(a[2..-1] + [l])
    }

    // Private helper method for modMul method.
    static checkedAdd_(a, b, c) {
        var room = c - 1 - a
        if (b <= room) {
            a = a + b
        } else {
            a = b - room - 1
        }
        return a
    }

    // Returns 'x' multiplied by 'n' modulo 'm'.
    static modMul(x, n, m) {
        if (m == 0) Fiber.abort("Cannot take modMul with modulus 0.")
        var sign = x.sign * n.sign
        var z = 0
        m = m.abs
        x = x.abs % m
        n = n.abs % m
        if (n > x) {
            var t = x
            x = n
            n = t
        }
        while (n > 0) {
            if (n % 2 == 1) z = checkedAdd_(z, x, m)
            x = checkedAdd_(x, x, m)
            n = div(n, 2)
        }
        return z * sign
    }

    // Returns the remainder when 'b' raised to the power 'e' is divided by 'm'.
    static modPow(b, e, m) {
        if (m == 0) Fiber.abort("Cannot take modPow with modulus 0.")
        if (m == 1 || m == -1) return 0
        var r = 1
        b = b % m
        if (e < 0) {
            e = -e
            b = modInv(b, m)
        }
        while (e > 0) {
            if (b == 0) return 0
            if (e % 2 == 1) r = modMul(r, b, m)
            e = div(e, 2)
            b = modMul(b, b, m)
        }
        return r
    }

    // Returns the multiplicative inverse of 'x'  modulo 'n'.
    static modInv(x, n) {
        var r = n
        var newR = x.abs
        var t = 0
        var newT = 1
        while (newR != 0) {
            var q = quo(r, newR)
            var lastT = t
            var lastR = r
            t = newT
            r = newR
            newT = lastT - q * newT
            newR = lastR - q * newR
        }
        if (r != 1) Fiber.abort("%(x) and %(n) are not co-prime.")
        if (t < 0) t = t + n
        return (x < 0) ? -t : t
    }

    // Returns the factorial of 'n'.
    static factorial(n) {
        if (!(n is Num && n >= 0 && n < 19)) {
            Fiber.abort("Argument must be a non-negative integer < 19.")
        }
        if (n < 2) return 1
        var fact = 1
        for (i in 2..n) fact = fact * i
        return fact
    }

    // Returns the multinomial coefficient of n over a list f where sum(f) == n.
    static multinomial(n, f) {
        if (!(n is Num && n >= 0 && n < 19)) {
            Fiber.abort("First argument must be a non-negative integer < 19.")
        }
        if (!(f is List)) Fiber.abort("Second argument must be a list.")
        var sum = f.reduce { |acc, i| acc + i }
        if (n != sum) {
            Fiber.abort("The elements of the list must sum to 'n'.")
        }
        var prod = 1
        for (e in f) {
            if (e < 0) Fiber.abort("The elements of the list must be non-negative integers.")
            if (e > 1) prod = prod * factorial(e)
        }
        return factorial(n)/prod
    }

    // Returns the binomial coefficent of n over k.
    static binomial(n, k) { multinomial(n, [k, n-k]) }

    // The highest prime less than 2^53.
    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 isPrimeMR_(n) {
        if (n > Num.maxSafeInteger) Fiber.abort("Argument must be a 'safe' integer.")
        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 < 4759123141) {
            a = [2, 7, 61]
        } else if (n < 1122004669633) {
            a = [2, 13, 23, 1662803]
        } else if (n < 2152302898747) {
            a = [2, 3, 5, 7, 11]
        } else if (n < 3474749660383) {
            a = [2, 3, 5, 7, 11, 13]
        } else if (n < 341550071728321) {
            a = [2, 3, 5, 7, 11, 13, 17]
        } else {
            a = [2, 3, 5, 7, 11, 13, 17, 19, 23]
        }
        var nPrev = n - 1
        var b = nPrev
        var r = 0
        while (b % 2 == 0) {
            b = div(b, 2)
            r = r + 1
        }
        for (i in 0...a.count) {
            if (n >= a[i]) {
                var x = modPow(a[i], b, n)
                if (x != 1 && x != nPrev) {
                    var d = r - 1
                    var next = false
                    while (d != 0) {
                        x = modMul(x, x, n)
                        if (x == 1) return false
                        if (x == nPrev) {
                            next = true
                            break
                        }
                        d = d - 1
                    }
                    if (!next) return false
                }
            }
        }
        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
    }

    // 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
        while (true) {
            if (Int.isPrime(n)) return n
            n = n + 2
        }
    }

    // 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 < 3) return null
        if (n == 3) return 2
        n = (n%2 == 0) ? n - 1 : n - 2
        while (true) {
            if (Int.isPrime(n)) return n
            n = n - 2
        }
    }

    // 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
    }

    // Sieves for primes up to and including 'n'.
    // If primesOnly is true returns a list of the primes found.
    // If primesOnly is false returns a bool list 'c' of size (n + 1) where:
    // c[i] is false if 'i' is prime or true if 'i' is composite.
    static primeSieve(n, primesOnly) {
        if (n < 2) return []
        var primes = [2]
        var k = ((n-3)/2).floor + 1
        var marked = List.filled(k, true)
        var limit = ((n.sqrt.floor - 3)/2).floor + 1
        limit = limit.max(0)
        for (i in 0...limit) {
            if (marked[i]) {
                var p = 2*i + 3
                var s = ((p*p - 3)/2).floor
                var j = s
                while (j < k) {
                    marked[j] = false
                    j = j + p
                }
            }
        }
        for (i in 0...k) {
            if (marked[i]) primes.add(2*i + 3)
        }
        if (primesOnly) return primes
        var c = List.filled(n+1, true)
        for (p in primes) c[p] = false
        return c
    }

    // Convenience version of above method which uses true for the primesOnly parameter.
    static primeSieve(n) { primeSieve(n, true) }

    // Sieves for primes up to and including 'limit' using a segmented approach
    // and returns a list of the primes found in order.
    // Second argument needs to be the cache size (L1) of your machine in bytes.
    // Translated from public domain C++ code at:
    // https://gist.github.com/kimwalisch/3dc39786fab8d5b34fee
    static segmentedSieve(limit, cacheSize) {
        cacheSize = (cacheSize/8).floor // 8 bytes per list element
        var sqroot = limit.sqrt.floor
        var segSize = sqroot.max(cacheSize)
        if (limit < 2) return []
        var allPrimes = [2]
        var isPrime = List.filled(sqroot + 1, true)
        var primes  = []
        var multiples = []
        var i = 3
        var n = 3
        var s = 3
        var low = 0
        while (low <= limit) {
            var sieve = List.filled(segSize, true)
            var high = limit.min(low + segSize - 1)
            while (i * i <= high) {
                if (isPrime[i]) {
                    var j = i * i
                    while (j <= sqroot) {
                        isPrime[j] = false
                        j = j + i
                    }
                }
                i = i + 2
            }
            while (s * s <= high) {
                if (isPrime[s]) {
                    primes.add(s)
                    multiples.add(s * s - low)
                }
                s = s + 2
            }
            for (ii in 0...primes.count) {
                var j = multiples[ii]
                var k = primes[ii] * 2
                while (j < segSize) {
                    sieve[j] = false
                    j = j + k
                }
                multiples[ii] = j - segSize
            }
            while (n <= high) {
                if (sieve[n - low]) allPrimes.add(n)
                n = n + 2
            }
            low = low + segSize
        }
        return allPrimes
    }

    // Computes and returns how many primes there are up to and including 'n'.
    // See https://rosettacode.org/wiki/Legendre_prime_counting_function for details.
    static primeCount(n) {
        if (n < 2) return 0
        if (n < 9) return ((n + 1)/2).floor
        var masks = [1, 2, 4, 8, 16, 32, 64, 128]
        var half = Fn.new { |n| (n - 1) >> 1 }
        var rtlmt = n.sqrt.floor
        var mxndx = Int.quo(rtlmt - 1, 2)
        var arrlen = (mxndx + 1).max(2)
        var smalls = List.filled(arrlen, 0)
        var roughs = List.filled(arrlen, 0)
        var larges = List.filled(arrlen, 0)
        for (i in 0...arrlen) {
            smalls[i] = i
            roughs[i] = i + i + 1
            larges[i] = Int.quo(Int.quo(n, i + i + 1) - 1 , 2)
        }
        var cullbuflen = Int.quo(mxndx + 8, 8)
        var cullbuf = List.filled(cullbuflen, 0)
        var nbps = 0
        var rilmt = arrlen
        var i = 1
        while (true) {
            var sqri = (i + i) * (i + 1)
            if (sqri > mxndx) break
            if ((cullbuf[i >> 3] & masks[i & 7]) != 0) {
                i = i + 1
                continue
            }
            cullbuf[i >> 3] = cullbuf[i >> 3] | masks[i & 7]
            var bp = i + i + 1
            var c = sqri
            while (c < arrlen) {
                cullbuf[c >> 3] = cullbuf[c >> 3] | masks[c & 7]
                c = c + bp
            }
            var nri = 0
            for (ori in 0...rilmt) {
                var r = roughs[ori]
                var rci = r >> 1
                if ((cullbuf[rci >> 3] & masks[rci & 7]) != 0) continue
                var d = r * bp
                var t = (d <= rtlmt) ? larges[smalls[d >> 1] - nbps] :
                                       smalls[half.call(Int.quo(n, d))]
                larges[nri] = larges[ori] - t + nbps
                roughs[nri] = r
                nri = nri + 1
            }
            var si = mxndx
            var pm = ((rtlmt/bp).floor - 1) | 1
            while (pm >= bp) {
                var c = smalls[pm >> 1]
                var e = (pm * bp) >> 1
                while (si >= e) {
                    smalls[si] = smalls[si] - c + nbps
                    si = si - 1
                }
                pm = pm - 2
            }
            rilmt = nri
            nbps = nbps + 1
            i = i + 1
        }
        var ans = larges[0] + ((rilmt + 2 * (nbps - 1)) * (rilmt - 1) / 2).floor
        for (ri in 1...rilmt) ans = ans - larges[ri]
        var ri = 1
        while (true) {
            var p = roughs[ri]
            var m = Int.quo(n, p)
            var ei = smalls[half.call(Int.quo(m, p))] - nbps
            if (ei <= ri) break
            ans = ans - (ei - ri) * (nbps + ri - 1)
            for (sri in (ri + 1..ei)) {
                ans = ans + smalls[half.call(Int.quo(m, roughs[sri]))]
            }
            ri = ri + 1
        }
        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
    }

    // Returns the prime factors of 'n' in order using a wheel with basis [2, 3, 5].
    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 = []
        while (n%2 == 0) {
            factors.add(2)
            n = (n/2).truncate
        }
        while (n%3 == 0) {
            factors.add(3)
            n = (n/3).truncate
        }
        while (n%5 == 0) {
            factors.add(5)
            n = (n/5).truncate
        }
        var k = 7
        var i = 0
        while (k * k <= n) {
            if (n%k == 0) {
                factors.add(k)
                n = (n/k).truncate
            } else {
                k = k + inc[i]
                i = (i + 1) % 8
            }
        }
        if (n > 1) factors.add(n)
        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 = []
        var i = 1
        var k = (n%2 == 0) ? 1 : 2
        while (i <= n.sqrt) {
            if (n%i == 0) {
                divisors.add(i)
                var j = (n/i).floor
                if (j != i) divisors2.add(j)
            }
            i = i + k
        }
        if (!divisors2.isEmpty) divisors = divisors + divisors2[-1..0]
        return divisors
    }

    // 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) {
        var pf = primeFactors(n)
        if (pf.count == 0) return (n == 1) ? [1] : pf
        var arr = []
        if (pf.count == 1) {
            arr.add([pf[0], 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, count])
                    prevPrime = pf[i]
                    count = 1
                }
            }
            arr.add([prevPrime, count])
        }
        var divisors = []
        var generateDivs
        generateDivs = Fn.new { |currIndex, currDivisor|
            if (currIndex == arr.count) {
                divisors.add(currDivisor)
                return
            }
            for (i in 0..arr[currIndex][1]) {
                generateDivs.call(currIndex+1, currDivisor)
                currDivisor = currDivisor * arr[currIndex][0]
            }
        }
        generateDivs.call(0, 1)
        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.isInteger || n < 1) return 0
        if (n > Num.maxSafeInteger) Fiber.abort("Argument must be a 'safe' integer.")
        var total = 1
        var power = 2
        while (n % 2 == 0) {
            total = total + power
            power = 2 * power
            n = div(n, 2)
        }
        var i = 3
        while (i * i <= n) {
            var sum = 1
            power = i
            while (n % i == 0) {
                sum = sum + power
                power = i * power
                n = div(n, i)
            }
            total = total * sum
            i = i + 2
        }
        if (n > 1) total = total * (n + 1)
        return total
    }

    // Returns the number of divisors of 'n' including 1 and 'n' itself.
    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
        while (n % 2 == 0) {
            count = count + 1
            n = div(n, 2)
        }
        prod = prod * (1 + count)
        var i = 3
        while (i * i <= n) {
            count = 0
            while (n % i == 0) {
                count = count + 1
                n = div(n, i)
            }
            prod = prod * (1 + count)
            i = i + 2
         }
         if (n > 2) prod = prod * 2
         return prod
    }

    // 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)) {
            Fiber.abort("Base must be an integer between 2 and 63.")
        }
    }

    // Returns a list of an integer n's digits in base b. Optionally checks n and b are valid.
    static digits(n, b, check) {
        if (check) check_(n, b)
        if (n == 0) return [0]
        var digs = []
        while (n > 0) {
            digs.add(n%b)
            n = (n/b).floor
        }
        return digs[-1..0]
    }

    // Returns the sum of an integer n's digits in base b. Optionally checks n and b are valid.
    static digitSum(n, b, check) {
        if (check) check_(n, b)
        var sum = 0
        while (n > 0) {
            sum = sum + (n%b)
            n = (n/b).floor
        }
        return sum
    } 

    // Returns the digital root and additive persistence of an integer n in base b.
    // Optionally checks n and b are valid.
    static digitalRoot(n, b, check) {
        if (check) check_(n, b)
        var ap = 0
        while (n > b - 1) {
            n = digitSum(n, b)
            ap = ap + 1
        }
        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
        while (m > 0) {
            r = 10*r + m%10
            m = div(m, 10)
        }
        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 
    // of non-negative integers 'x' and 'y' according to Cantor's pairing function. 
    static cantorPair(x, y) {
        if (x.type != Num || !x.isInteger || x < 0) {
            Fiber.abort("Arguments must be non-negative integers.")
        }
        if (y.type != Num || !y.isInteger || y < 0) {
            Fiber.abort("Arguments must be non-negative integers.")
        }
        return (x*x + 3*x + 2*x*y + y + y*y) / 2
    }

    // Returns the pair of non-negative integers that are associated with a single 
    // non-negative integer 'z' according to Cantor's pairing function. 
    static cantorUnpair(z) {
        if (z.type != Num || !z.isInteger || z < 0) {
            Fiber.abort("Argument must be a non-negative integer.")
        }
        var i = (((1 + 8*z).sqrt-1)/2).floor      
        return [z - i*(1+i)/2, i*(3+i)/2 - z]
    }
}

/*
    Nums contains various routines applicable to lists or ranges of numbers
    many of which are useful for statistical purposes.
*/
class Nums {
    // Methods to calculate sum, various means, product and maximum/minimum element of 'a'.
    // The sum and product of an empty list are considered to be 0 and 1 respectively.
    static sum(a)  { a.reduce(0) { |acc, x| acc + x } }
    static mean(a) { sum(a)/a.count }
    static geometricMean(a) { a.reduce { |prod, x| prod * x}.pow(1/a.count) }
    static harmonicMean(a) { a.count / a.reduce { |acc, x| acc + 1/x } }
    static quadraticMean(a) { (a.reduce(0) { |acc, x| acc + x*x }/a.count).sqrt }
    static prod(a) { a.reduce(1) { |acc, x| acc * x } }
    static max(a)  { a.reduce { |acc, x| (x > acc) ? x : acc } }
    static min(a)  { a.reduce { |acc, x| (x < acc) ? x : acc } }

    // As above methods but applying a function 'f' to each element of 'a' 
    // before performing the operation. 
    // 'f' should take a single Num parameter and return a Num.
    static sum(a, f)  { a.reduce(0) { |acc, x| acc + f.call(x) } }
    static mean(a, f) { sum(a, f)/a.count }
    static geometricMean(a, f) { a.reduce { |prod, x| prod * f.call(x)}.pow(1/a.count) }
    static harmonicMean(a, f) { a.count / a.reduce { |acc, x| acc + 1/f.call(x) } }
    static quadraticMean(a, f) { (a.reduce(0) { |acc, x| acc + f.call(x).pow(2) }/a.count).sqrt }
    static prod(a, f) { a.reduce(1) { |acc, x| acc * f.call(x) } }
    static max(a, f)  { a.reduce { |acc, x|
         var fx = f.call(x)
         return (fx > acc) ? fx : acc
    } }
    static min(a, f)  { a.reduce { |acc, x| 
         var fx = f.call(x)
         return (fx < acc) ? fx : acc 
    } }

    // Returns the median of a sorted list 'a'.
    static median(a) {
        var c = a.count
        if (c == 0) {
            Fiber.abort("An empty list cannot have a median")
        } else if (c%2 == 1) {
            return a[(c/2).floor]
        } else {
            var d = (c/2).floor
            return (a[d] + a[d-1])/2
        }
    }

    // Returns a list whose first element is a list of the mode(s) of 'a'
    // and whose second element is the number of times the mode(s) occur.
    static modes(a) {
        var m = {}
        for (e in a) m[e] = (!m[e]) ? 1 : m[e] + 1
        var max = 0
        for (e in a) if (m[e] > max) max = m[e]
        var res = []
        for (k in m.keys) if (m[k] == max) res.add(k)
        return [max, res]
    }

    // Returns the sample variance of 'a'.
    static variance(a) {
        var m = mean(a)
        var c = a.count
        return (a.reduce(0) { |acc, x| acc + x*x } - m*m*c) / (c-1)
    }

    // Returns the population variance of 'a'.
    static popVariance(a) {
        var m = mean(a)
        return (a.reduce(0) { |acc, x| acc + x*x }) / a.count - m*m
    }

    // Returns the sample standard deviation of 'a'.
    static stdDev(a) { variance(a).sqrt }

    // Returns the population standard deviation of 'a'.
    static popStdDev(a) { popVariance(a).sqrt }

    // Returns the mean deviation of 'a'.
    static meanDev(a) {
        var m = mean(a)
        return a.reduce { |acc, x| acc + (x - m).abs } / a.count
    }

    // Converts a string list to a corresponding numeric list.
    static fromStrings(a) { a.map { |s| Num.fromString(s) }.toList }

    // 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, "▮") }
}

/* Boolean supplements the Bool class with bitwise operations on boolean values. */
class Boolean {
    // Private helper method to convert a boolean to an integer.
    static btoi_(b) { b ? 1 : 0 }

    // Private helper method to convert an integer to a boolean.
    static itob_(i) { i != 0 }

    // Private helper method to check its arguments are both booleans.
    static check_(b1, b2) {
        if (!((b1 is Bool) && (b2 is Bool))) Fiber.abort("Both arguments must be booleans.")
    }

    // Returns the logical 'and' of its boolean arguments.
    static and(b1, b2) {
        check_(b1, b2)
        return itob_(btoi_(b1) & btoi_(b2))
    }

    // Returns the logical 'or' of its boolean arguments.
    static or(b1, b2) {
        check_(b1, b2)
        return itob_(btoi_(b1) | btoi_(b2))
    }

    // Returns the logical 'xor' of its boolean arguments.
    static xor(b1, b2) {
        check_(b1, b2)
        return itob_(btoi_(b1) ^ btoi_(b2))
    }
}