Category talk:Wren-long

From Rosetta Code

64-bit integer arithmetic[edit]

Many tasks on RC require the use of large numbers (usually integers).

Wren is unable to deal natively and precisely with integers outside a 53-bit range and so either has to limit such tasks to integers which it can deal with or use the BigInt class in the Wren-big module which can deal with integers of arbitrary size.

However, some tasks on RC only require integers within the 64-bit range and BigInt feels like overkill for such tasks.

This module aims to remedy that situation by providing ULong and Long classes for 64-bit unsigned and signed arithmetic, respectively, which are easy to use.

ULong[edit]

Internally, a 64-bit unsigned integer is represented by two unsigned 32-bit integral Nums, a low and a high part, though this is generally invisible to the user.

The design of ULong follows closely that of the BigInt class though the implementation is generally much simpler and hopefully more performant for numbers within its range.

The ULong.new constructor accepts either a string or an unsigned safe integer (i.e. within the normal 53-bit limitation). I couldn't see much point in accepting larger numbers (or floats) as this would inevitably mean that the ULong object would be imprecise from the outset. It is an error to pass a value which is out of range. Passing an empty string will also produce an error (rather than a zero object) as it's just as easy to pass 0 and an empty string may well be a mistake in any case.

It's also possible to create a ULong directly from a pair of unsigned 32-bit integers, low and high, or from a little-endian list of 8 bytes using the fromPair or fromBytes constructors respectively.

If you need to generate a ULong from a string in a non-decimal base or from another ULong, you can use the fromBaseString method for any base, fromHexString for base 16 or copy methods respectively. However, the first two - unlike the new constructor - simply 'wrap around' out of range values (even for bases 10 or 16) as doing anything else was felt to be too expensive.

Similarly and where appropriate, all operations on ULongs 'wrap around' out of range results not just because of the expense of doing anything else but also because this is often more convenient than throwing errors.

Compared to BigInt, most operations are relatively trivial to implement the exceptions being multiplication and division/remainder.

Where the operands are small enough, I've simply used Num multiplication or division.

However, in the case of multiplication of larger numbers, one can't just multiply all the 32-bit portions and add them together because of the 53-bit limitation. What I've done instead, is to split them into 16-bit portions to do the multiplication and then combined the results using bit arithmetic.

The most obvious way to perform division of larger numbers is to first obtain an estimate of the result by treating them as Nums and then refine this result iteratively until the exact result is obtained. However, this is easier said than done because the accuracy of the initial estimate will depend on the magnitude of the numbers and one has to guard against overflow when performing the iterations. Consequently, this is quite a slow operation.

Long[edit]

Internally, a 64-bit signed integer is represented by a ULong (its magnitude) and a sign (-1 for negative, 0 for zero and +1 for positive integers). This design has the following advantages:

1. Longs have double the range of ULongs and not one bit less as is usually the case for 64-bit signed integers where one bit is needed to represent the sign.

2. Apart from manipulation of the sign, most of the methods for the Long type can simply defer to the corresponding methods for the ULong type which makes the code for the Long class much shorter and simpler than it otherwise would have been.

3. There is no need to support bitwise operations and other stuff which is only really relevant to unsigned integers such as prime number routines. This again shortens and simplifies the code.

4. It is relatively easy to mix both Longs and ULongs in arithmetic expressions by performing the necessary conversions automatically.

However, a big disadvantage is that Longs use 50% more memory (3 Nums instead of two) and are slightly slower than ULongs which means that the latter should normally be preferred unless you know or suspect that you will need to deal with negative numbers.

A peculiarity of the design (though not necessarily a disadvantage) is that an overflow wraps around to zero rather than the opposite end of the range as is usual for 64-bit signed integers.

Source code[edit]

/* Module "long.wren" */

import "./trait" for Comparable

/*
    ULong represents a 64-bit unsigned integer together with arithmetic operations thereon.
    ULong objects are stored as two non-negative integers < 2^32 and are immutable.
*/
class ULong is Comparable {
    // Constants
    static zero     { lohi_( 0, 0) }
    static one      { lohi_( 1, 0) }
    static two      { lohi_( 2, 0) }
    static three    { lohi_( 3, 0) }
    static four     { lohi_( 4, 0) }
    static five     { lohi_( 5, 0) }
    static ten      { lohi_(10, 0) }

    // Returns the maximum 'short' ULong = 2^32-1 = 4294967295
    static maxShort { lohi_(4294967295, 0) }

    // Returns the maximum 'small' ULong = 2^53-1 = 9007199254740991
    static maxSmall { lohi_(4294967295, 2097151) }

    // Returns the largest ULong = 2^64-1 = 18446744073709551615
    static largest  { lohi_(4294967295, 4294967295) }

    // Returns the smallest ULong = 0
    static smallest { lohi_(0, 0) }

    // Private method to determine whether a number is a short unsigned integer or not.
    static isShort_(n) { (n is Num) && n.isInteger && n >= 0 && n <= 4294967295 }
    
    // Private method to determine whether a number is a small unsigned integer or not.
    static isSmall_(n) { (n is Num) && n.isInteger && n >= 0 && n <= 9007199254740991 }

    // Private method to determine whether a number is a non-negative Long or not.
    static isNonNegLong_(n) { (n is Long) && !n.isNegative }

    // Private method to calculate the base 2 logarithm of a number.
    static log2_(n) { n.log / 0.69314718055994530942 }

    // Private helper function to convert a string to lower case.
    static lower_(s) { s.codePoints.map { |c|
        return String.fromCodePoint((c >= 65 && c <= 90) ? c + 32 : c)
    }.join() }


    // Private helper method to convert a lower case base string to a small integer.
    static atoi_(s, base, digits) {
        var res = 0
        for (d in s) {
            var ix = digits.indexOf(d)
            res = res * base + ix
        }
        return res
    }

    // Private helper method to convert a small integer to a string with a base between 2 and 36.
    static itoa_(n, base, digits) {
        if (n == 0) return "0"
        var res = ""
        while (n > 0) {
            res = res + "%(digits[n%base])"
            n = (n/base).floor
        }
        return res[-1..0]
    }

    // Private helper method to check whether a 20 digit string exceeds the largest value string.
    static exceedsMaxString_(s) {
        var m = "18446744073709551615"
        if (s == m) return false
        var mb = m.bytes
        var sb = s.bytes
        for (i in 0..19) {
            if (mb[i] < sb[i]) return true
            if (mb[i] > sb[i]) return false
        }
    }

    // Returns the greater of two ULongs.
    static max(a, b) {
        if (!(a is ULong)) a = new(a)
        if (!(b is ULong)) b = new(b)
        return (a > b) ? a : b
    }

    // Returns the lesser of two ULongs.
    static min(a, b) {
        if (!(a is ULong)) a = new(a)
        if (!(b is ULong)) b = new(b)
        return (a < b) ? a : b
    }

    // Returns the greatest common divisor of a and b.
    static gcd(a, b) {
        if (!(a is ULong)) a = new(a)
        if (!(b is ULong)) b = new(b)
        while (!b.isZero) {
            var t = b
            b = a % b
            a = t
        }
        return a
    }

    // Returns the least common multiple of a and b.
    static lcm(a, b) {
        if (!(a is ULong)) a = new(a)
        if (!(b is ULong)) b = new(b)
        if (a.isZero && b.isZero) return ULong.zero 
        return a / gcd(a, b) * b
    }

    // Returns the factorial of 'n'. Can only be used for n <= 20
    static factorial(n) {
        if (!(((n is Num) || (n is ULong)) && n >= 0 && n <= 20)) {
            Fiber.abort("Argument must be a non-negative integer no larger than 20.")
        }
        if (n < 2) return one
        var fact = one
        var i = two
        while (i <= n) {
            fact = fact * i
            i = i + 1
        }
        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 <= 20)) {
            Fiber.abort("First argument must be a non-negative integer <= 20.")
        }
        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 = one
        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]) }

    // Returns whether or not 'n' is an instance of ULong.
    static isInstance(n) { n is ULong }

    // Private method to determine if a ULong is a basic prime or not.
    static isBasicPrime_(n) {
        if (!(n is ULong)) n = new(n)
        if (n.isOne) return false
        if (n == two || n == three || n == five) return true
        if (n.isEven || n.isDivisibleBy(three) || n.isDivisibleBy(five)) return false
        if (n < 49) return true
        return null // unknown if prime or not
    }

    // Private method to apply the Miller-Rabin test.
    static millerRabinTest_(n, a) {
        var nPrev = n.dec
        var b = nPrev
        var r = 0
        while (b.isEven) {
            b = b >> 1
            r = r + 1
        }
        for (i in 0...a.count) {
            if (n >= a[i]) {
                var x = (a[i] is ULong) ? a[i] : new(a[i])
                x = x.modPow(b, n)
                if (!x.isOne && x != nPrev) {
                    var d = r - 1
                    var next = false
                    while (d != 0) {
                        x = x.square % n
                        if (x.isOne) return false
                        if (x == nPrev) {
                            next = true
                            break
                        }
                        d = d - 1
                    }
                    if (!next) return false
                }
            }
        }
        return true
    }

    // Private constructor which creates a ULong object from low and high components.
    construct lohi_(low, high) {
        _lo = low
        _hi = high
    }

    // Private constructor which creates a ULong object from a 'small' integer.
    construct fromSmall_(v) {
        var p = 4294967296  // 2 ^ 32
        _lo = (v % p)
        _hi = (v / p).floor
    }

    // Private method which creates a ULong object from a Num.
    // If 'v' is not small, will probably lose accuracy.
    static fromNum_(v) {
        if (v < 0 || v.isNan) return ULong.zero
        var m = 4294967296  // 2 ^ 32
        if (v >= 2.pow(64)) return ULong.lohi_(m - 1, m - 1)
        return ULong.lohi_(v % m, (v / m).floor)
    }

    // Private method which creates a ULong object from a base 10 numeric string.
    // Scientific notation is permitted.
    // Raises an error if the result is out of bounds.
    static fromString_(v) {
        v = v.trim()
        if (v.count == 0 || v[0] == "-") Fiber.abort("Invalid unsigned integer.")
        if (v[0] == "+") {
            v = v[1..-1]
            if (v.count == 0) Fiber.abort("Invalid unsigned integer.")
        }
        v = v.trimStart("0")
        if (v == "") v = "0"
        v = ULong.lower_(v)
        var split = v.split("e")
        if (split.count > 2) Fiber.abort("Invalid unsigned integer.")
        if (split.count == 2) {
            var exp = split[1]
            if (exp[0] == "+") exp = exp[1..-1]
            exp = Num.fromString(exp)
            if (!ULong.isSmall_(exp)) Fiber.abort("Exponent is not valid.")
            var text = split[0]
            var dp = text.indexOf(".")
            if (dp >= 0) {
                exp = exp - (text.count - dp - 1)
                text = text[0...dp] + text[dp+1..-1]
            }
            if (exp < 0) Fiber.abort("Exponent cannot be negative.")
            text = text + ("0" * exp)
            v = text
        }
        var len = v.count
        var isValid = len > 0 && v.all { |d| "0123456789".contains(d) }
        if (!isValid) Fiber.abort("Invalid unsigned integer.")
        if (len > 20 || (len == 20 && ULong.exceedsMaxString_(v))) {
            Fiber.abort("Integer is too big.")
        }
        if (len <= 16) {
            var n = Num.fromString(v)
            if (ULong.isSmall_(n)) return ULong.fromSmall_(n)
        }
        // process in 10 digit chunks
        var r = ULong.zero
        var pow10 = ULong.fromSmall_(10.pow(10))
        var i = 0
        while (i < len) {
            var chunkSize = ((len - i) < 10) ? len - i : 10
            var chunk = Num.fromString(v[i...i + chunkSize])
            if (chunkSize < 10) {
                var psize = ULong.fromSmall_(10.pow(chunkSize))
                r = r * psize + ULong.fromSmall_(chunk)
            } else {
                r = r * pow10 + ULong.fromSmall_(chunk)
            }
            i = i + 10
        }
        return r
    }

    // Creates a ULong object from an (unprefixed) numeric string in a given base (2 to 36).
    // Scientific notation is not permitted.
    // Wraps out of range values.
    static fromBaseString(v, base) {
        if (!(v is String)) Fiber.abort("Value must be a numeric string in the given base.")
        if (!((base is Num) && base.isInteger && base >= 2 && base <= 36)) {
            Fiber.abort("Base must be an integer between 2 and 36.")
        }
        v = v.trim()
        if (v.count == 0 || v[0] == "-") Fiber.abort("Invalid unsigned integer.")
        if (v[0] == "+") {
            v = v[1..-1]
            if (v.count == 0) Fiber.abort("Invalid unsigned integer.")
        }
        v = v.trimStart("0")
        if (v == "") v = "0"
        if (base > 10) v = ULong.lower_(v)
        var alphabet = "0123456789abcdefghijklmnopqrstuvwxyz"
        var digits = alphabet[0...base]
        var len = v.count
        var isValid = len > 0 && v.all { |d| digits.contains(d) }
        if (!isValid) Fiber.abort("Invalid unsigned integer.")
        // process in 10 digit chunks
        var r = ULong.zero
        var powb = ULong.fromSmall_(base.pow(10))
        var i = 0
        while (i < len) {
            var chunkSize = ((len - i) < 10) ? len - i : 10
            var chunk = ULong.atoi_(v[i...i + chunkSize], base, digits)
            if (chunkSize < 10) {
                var psize = ULong.fromSmall_(base.pow(chunkSize))
                r = r * psize + ULong.fromSmall_(chunk)
            } else {
                r = r * powb + ULong.fromSmall_(chunk)
            }
            i = i + 10
        }
        return r
    }    

    // Creates a ULong object from either a numeric base 10 string, an unsigned 'small' integer or a non-negative Long.
    static new(value) {
         if (!(value is String) && !isSmall_(value) && !isNonNegLong_(value)) {
              Fiber.abort("Value must be a base 10 numeric string, an unsigned small integer or a non-negative Long.")
         }
         return (value is String) ? fromString_(value) : (value is Num) ? fromSmall_(value) : value.toULong
    }

    // Creates a ULong object from an (unprefixed) hexadecimal string. Wraps out of range values. 
    static fromHexString(v) { fromBaseString(v, 16) }

    // Creates a ULong object from a pair (low and high) of unsigned 'short' integers.
    static fromPair(low, high) {
        if (!isShort_(low) || !isShort_(high)) {
            Fiber.abort("Low and high components must both be unsigned 32-bit integers.")
        }
        return lohi_(low, high)
    }

    // Creates a ULong object from a list of 8 unsigned bytes in little-endian format.
    static fromBytes(bytes) {
        if (bytes.count != 8) Fiber.abort("There must be exactly 8 bytes in the list.")
        for (b in bytes) {
            if (!(b is Num && b.isInteger && b >= 0 && b < 256)) {
                Fiber.abort("Each byte must be an integer between 0 and 255.")
            }
        }
        var low  = bytes[0] | bytes[1] << 8 | bytes[2] << 16 | bytes[3] << 24
        var high = bytes[4] | bytes[5] << 8 | bytes[6] << 16 | bytes[7] << 24
        return lohi_(low, high)
    }

    // Properties to return the low and high 32-bit portions of this instance.
    low  { _lo }
    high { _hi }

    // Public self-evident properties.
    isShort { _hi == 0 }
    isSmall { _hi <= 2097151 }
    isEven  { _lo % 2 == 0 }
    isOdd   { _lo % 2 == 1 }
    isOne   { _lo == 1 && _hi == 0 }
    isZero  { _lo == 0 && _hi == 0 }

    // Returns true if 'n' is a divisor of the current instance, false otherwise
    isDivisibleBy(n) { (this % n).isZero }

    // Returns true if the current instance is prime, false otherwise.
    // Fails due to overflow on numbers > 9223372036854775807 unless divisible by 2,3 or 5.
    isPrime {
        var isbp = ULong.isBasicPrime_(this)
        if (isbp != null) return isbp
        if (this > ULong.largest/2) Fiber.abort("Cannot test %(this) for primality.")
        return ULong.millerRabinTest_(this, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37])
    }

    // Returns the next prime number greater than the current instance.
    nextPrime {
        var n = isEven ? this + ULong.One : this + ULong.two
        while (true) {
            if (n.isPrime) return n
            n = n + ULong.two
        }
    }

    // Private method to calculate the index of this instance's most significant bit (0 to 63).
    msb_ {
        if (_hi == 0) return (_lo > 0) ? ULong.log2_(_lo).floor : 0
        return ULong.log2_(_hi).floor + 32
    }

    // Returns the bitwise complement of the current instance.
    ~ { ULong.lohi_(~_lo, ~_hi) }

    // Adds a ULong to the current instance. Wraps on overflow.
    + (n) {
        if (!(n is ULong)) n = ULong.new(n)
        if (this.isZero) return n.copy()
        if (n.isZero) return this.copy()
        var low  = _lo + n.low
        var high = _hi + n.high
        var m = 4294967296 // 2^32
        if (low >= m) {
            low = low - m
            high = high + 1
        }
        if (high >= m) {
            high = high - m
        }
        return ULong.lohi_(low, high)
    }

    // Subtracts a ULong from the current instance. Wraps on underflow.
    - (n) {
        if (!(n is ULong)) n = ULong.new(n)
        if (n.isZero) return this.copy()
        return this + (~n) + ULong.one
    }

    // Multiplies the current instance by a ULong. Wraps on overflow.
    * (n) {
        if (!(n is ULong)) n = ULong.new(n)
        if (this.isZero || n.isZero) return ULong.zero
        // if the sum of the msbs for both ULongs is less than 51 use Num multiplication
        if (this.msb_ + n.msb_ < 51) return ULong.fromSmall_(this.toNum * n.toNum)
        // otherwise split the operands into 16 bit pieces to do the multiplication
        var a3 = _hi >> 16
        var a2 = _hi & 0xffff
        var a1 = _lo >> 16
        var a0 = _lo & 0xffff

        var b3 = n.high >> 16
        var b2 = n.high & 0xffff
        var b1 = n.low >> 16
        var b0 = n.low & 0xffff

        var c0 = a0 * b0
        var c1 = (c0 >> 16) + a1 * b0
        c0 = c0 & 0xffff
        var c2 = c1 >> 16
        c1 = (c1 & 0xffff) + a0 * b1
        c2 = c2 + (c1 >> 16) + a2 * b0
        c1 = c1 & 0xffff
        var c3 = c2 >> 16
        c2 = (c2 & 0xffff) + a1 * b1
        c3 = c3 + (c2 >> 16)
        c2 = (c2 & 0xffff) + a0 * b2
        c3 = c3 + (c2 >> 16)
        c2 = c2 & 0xffff
        c3 = (c3 + a3 * b0 + a2 * b1 + a1 * b2 + a0 * b3) & 0xffff

        return ULong.lohi_((c1 << 16) | c0, (c3 << 16) | c2)
    }

    // Returns a list containing the quotient and the remainder after dividing
    // the current instance by a ULong.
    divMod(n) {
        if (!(n is ULong)) n = ULong.new(n)
        if (n.isZero) Fiber.abort("Cannot divide by zero.")
        // if both operands are 'small' use Num division.
        if (this.isSmall && n.isSmall) {
             var a = this.toNum
             var b = n.toNum
             return [ULong.fromSmall_((a/b).floor), ULong.fromSmall_(a%b)]
        }
        if (this.isZero) return [ULong.zero, ULong.zero]
        if (n.isOne) return [this.copy(), ULong.zero]
        if (n > this) return [ULong.zero, this.copy()]
        // use Num division to estimate the answer and refine it until the exact answer is found.
        var div = ULong.zero
        var rem = this.copy()
        // iterate until the remainder is less than the divisor
        while (rem >= n) {
            var est = (rem.toNum / n.toNum).floor
            if (est < 1) est = 1  // must be at least 1
            var emsb = ULong.log2_(est).ceil
            // calculate an adjustment to use based on the size of the estimate
            var adj = (emsb <= 53) ? 1 : 1 << (emsb-53)
            var div2 = ULong.fromNum_(est)
            var rem2 = div2 * n
            var rem3 = rem2 - div2 // to check whether rem2 has overflowed
            // reduce the estimated remainder until it is no greater than the actual remainder
            while (rem2 > rem || rem3 >= rem2) {
                est = est - adj
                div2 = ULong.fromNum_(est)
                rem2 = div2 * n
                rem3 = rem2 - div2
            }
            if (div2.isZero) div2 = ULong.one // must be at least one
            div = div + div2
            rem = rem - rem2
        }
        return [div, rem]
    }

    // Divides the current instance by a ULong.
    / (n) { divMod(n)[0] }

    // Returns the remainder after dividing the current instance by a ULong.
    % (n) { divMod(n)[1] }

    //Returns the bitwise 'and' of the current instance and another ULong.
    & (n) {
        if (!(n is ULong)) n = ULong.new(n)
        return ULong.lohi_(_lo & n.low, _hi & n.high)
    }

    // Returns the bitwise 'or' of the current instance and another ULong.
    | (n) {
        if (!(n is ULong)) n = ULong.new(n)
        return ULong.lohi_(_lo | n.low, _hi | n.high)
    }

    // Returns the bitwise 'xor' of the current instance and another ULong.
    ^ (n) {
        if (!(n is ULong)) n = ULong.new(n)
        return ULong.lohi_(_lo ^ n.low, _hi ^ n.high)
    }

    // Shifts the bits of the current instance 'n' places to the left. Wraps modulo 64.
    // Negative shifts are allowed.
    << (n) {
        if (n is ULong) n = n.toNum
        n = n & 63
        if (n == 0) return this.copy()
        if (n < 32) {
            return ULong.lohi_(_lo << n, (_hi << n) | (_lo >> (32 - n)))
        }
        return ULong.lohi_(0, _lo << (n - 32))
    }

    // Shifts the bits of the current instance 'n' places to the right. Wraps modulo 64.
    // Negative shifts are allowed.
    >> (n) {
        if (n is ULong) n = n.toNum
        n = n & 63
        if (n == 0) return this.copy()
        if (n < 32) {
            return ULong.lohi_(_lo >> n | (_hi << (32 - n)), _hi >> n)
        }
        return ULong.lohi_(_hi >> (n - 32), 0)
    }

    // The inherited 'clone' method just returns 'this' as ULong objects are immutable.
    // If you need an actual copy use this method instead.
    copy() { ULong.lohi_(_lo, _hi) }

    // Compares the current instance with a ULong. If they are equal returns 0.
    // If 'this' is greater, returns 1. Otherwise returns -1.
    // Also allows a comparison with positive infinity.
    compare(n) {
        if ((n is Num) && n.isInfinity && n > 0) return -1
        if (!(n is ULong)) n = ULong.new(n)
        if (_hi == n.high && _lo == n.low) return 0
        if (_hi > n.high) return 1
        if (_hi < n.high) return -1
        return (_lo < n.low) ? -1 : 1
    }

    // Returns the greater of this instance and another ULong instance.
    max(n) { ULong.max(this, n) }

    // Returns the smaller of this instance and another ULong instance.
    min(n) { ULong.min(this, n) }

    // Clamps this instance into the range [a, b].
    // If this instance is less than min, min is returned.
    // If it's more than max, max is returned. Otherwise, this instance is returned.
    clamp(a, b) {
        if (!(a is ULong)) a = ULong.new(a)
        if (!(b is ULong)) b = ULong.new(b)
        if (a > b) Fiber.abort("Range cannot be decreasing.")
        if (this < a) return a
        if (this > b) return b
        return this.copy()
    }

    // Squares the current instance. Wraps on overflow.
    square { this * this }

    // Cubes the current instance. Wraps on overflow.
    cube { this * this * this }

    // Returns true if the current instance is a perfect square, false otherwise.
    isSquare {
        var root = isqrt
        return root.square == this
    }

    // Returns true if the current instance is a perfect cube, false otherwise.
    isCube {
        var root = icbrt
        return root.cube == this
    }

    // Returns the integer n'th root of the current instance
    // i.e. the precise n'th root truncated towards zero if not an integer.
    iroot(n) {
        if (!((n is Num) && n.isInteger && n > 0)) {
            Fiber.abort("Argument must be a positive integer.")
        }
        if (n == 1) return this
        var t = copy()
        n = n - 1
        var s = t + 1
        var u = t
        while (u < s) {
            s = u
            u = ((u * n) + t / u.pow(n)) / (n + 1)
        }
        return s
    }

    // Returns the integer cube root of the current instance i.e. the largest integer 'x0'
    // such that x0.cube <= this.
    icbrt {
        if (isSmall) return ULong.fromSmall_(toNum.cbrt.floor)
        return iroot(3)
    }

    // Returns the integer square root of the current instance i.e. the largest integer 'x0'
    // such that x0.square <= this.
    isqrt {
        if (isSmall) return ULong.fromSmall_(toNum.sqrt.floor)
        // otherwise use Newton's method
        var x0 = this >> 1
        var x1 = (x0 + this/x0) >> 1
        while (x1 < x0) {
            x0 = x1
            x1 = (x0 + this/x0) >> 1
        }
        return x0
    }

    // Returns the current instance raised to the power of a 'small' ULong. Wraps on overflow.
    // If the exponent is less than 0, returns 0. O.pow(0) returns one.
    pow(n) {
        if (!(n is ULong)) n = ULong.new(n)
        if (n.isZero) return ULong.one
        if (n.isOne) return this.copy()
        if (this.isZero) return ULong.zero
        if (this.isOne) return ULong.one
        if (!n.isSmall) Fiber.abort("The exponent %(n) is too large.")
        if (this.isSmall) {
            var value = this.toNum.pow(n.toNum)
            if (ULong.isSmall_(value)) return ULong.fromSmall_(value)
        }
        var x = this
        var y = ULong.one
        var z = n
        while (true) {
            if (z.isOdd) {
                y = y * x
                z = z - 1
            }
            if (z.isZero) break
            z = z >> 1
            x = x.square
        }
        return y
    }

    // Returns the current instance multiplied by 'n' modulo 'mod'.
    modMul(n, mod) {
        var x = ULong.zero
        var y = this
        while (n > ULong.zero) {
            if ((n & 1) == 1) x = (x + y) % mod
            y = (y << 1) % mod
            n = n >> 1
        }
        return x
    }

    // Returns the current instance to the power 'exp' modulo 'mod'.
    modPow(exp, mod) {
        if (!(exp is ULong)) exp = ULong.new(exp)
        if (!(mod is ULong)) mod = ULong.new(mod)
        if (mod.isZero) Fiber.abort("Cannot take modPow with modulus 0.")
        var r = ULong.one
        var base = this % mod
        while (!exp.isZero) {
            if (base.isZero) return ULong.zero
            if (exp.isOdd) r = r.modMul(base, mod)
            exp = exp >> 1
            base = base.modMul(base, mod)
        }
        return r
    }

    // Increments the current instance by one.
    inc { this + ULong.one }

    // Decrements the current instance by one.
    dec { this - ULong.one }

    // Returns 0 if the current instance is zero or 1 otherwise.
    sign { isZero ? 0 : 1 }

    // Returns the number of digits required to represent the current instance in binary.
    bitLength { msb_ + 1 }

    // Returns true if the 'n'th bit of the current instance is set or false otherwise.
    testBit(n) {
        if (n.type != Num || !n.isInteger || n < 0 || n > 63) {
            Fiber.abort("Argument must be a non-negative integer less than 64.")
        }
        return (this >> n) & ULong.one != ULong.zero
    }

    // Converts the current instance to a Num where possible.
    // Will probably lose accuracy if the current instance is not 'small'.
    toNum { _hi * 4294967296 + _lo }

    // Converts the current instance to a 'small' integer where possible.
    // Otherwise returns null.
    toSmall { isSmall ? toNum : null }

    // Expresses the current instance as a pair of 'short' integers, low and high.
    toPair { [_lo, _hi] }

    // Expresses the current instance as a list of 8 unsigned bytes in little-endian format.
    toBytes {
        return [_lo & 0xff, _lo >> 8 & 0xff, _lo >> 16 & 0xff, _lo >> 24      ,
                _hi & 0xff, _hi >> 8 & 0xff, _hi >> 16 & 0xff, _hi >> 24]
    }

    // Converts the current instance to a Long.
    toLong { Long.sigma_(sign, this) }

    // Private worker method for toBaseString, toHexString and toString.
    toBaseString_(base) {
        if (isZero) return "0"
        // process in 6 digit chunks
        var pow6 = ULong.fromSmall_(base.pow(6))
        var alphabet = "0123456789abcdefghijklmnopqrstuvwxyz"
        var rem = this
        var res = ""
        while (true) {
            var div = rem / pow6
            var val = (rem - div * pow6).toNum >> 0
            var digits = ULong.itoa_(val, base, alphabet[0...base])
            rem = div
            if (rem.isZero) return digits + res
            if (digits.count < 6) digits = "0" * (6 - digits.count) + digits
            res = digits + res
        }
    }

    // Returns the string representation of the current instance in a given base (2 to 36).
    toBaseString(base) {
        if (!((base is Num) && base.isInteger && base >= 2 && base <= 36)) {
            Fiber.abort("Base must be an integer between 2 and 36.")
        }
        return toBaseString_(base)
    }

    // Returns the string representation of the current instance in base 16.
    toHexString { toBaseString_(16) }

    // Returns the string representation of the current instance in base 10.
    toString { toBaseString_(10) }

    /* Prime factorization methods. */

    // Private worker method for Pollard's Rho algorithm.
    static pollardRho_(m, seed, c) {
        var g = Fn.new { |x| (x*x + c) % m }
        var x = ULong.new(seed)
        var y = ULong.new(seed)
        var z = ULong.one
        var d = ULong.one
        var count = 0
        while (true) {
            x = g.call(x)
            y = g.call(g.call(y))
            d = (x >= y) ? (x - y) % m : (y - x) % m
            z = z * d
            count = count + 1
            if (count == 100) {
                d = ULong.gcd(z, m)
                if (d != ULong.one) break
                z = ULong.one
                count = 0
            }
        }
        if (d == m) return ULong.zero
        return d
    }

    // Returns a factor (ULong) of 'm' (a ULong or an unsigned 'small' integer) using the
    // Pollard's Rho algorithm. Both the 'seed' and 'c' can be set to integers.
    // Returns ULong.zero in the event of failure.
    static pollardRho(m, seed, c) {
        if (m < 2) return ULong.zero
        if (m is Num) m = ULong.new(m)
        if (m.isPrime) return m.copy()
        if (m.isSquare) return m.isqrt
        for (p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]) {
            if (m.isDivisibleBy(p)) return ULong.new(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 (ULong) 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 = ULong.new(37)
        var i = 0
        while (k * k <= n) {
            if (n.isDivisibleBy(k)) {
                factors.add(k.copy())
                n = n / k
            } else {
                k = k + inc[i]
                i = (i + 1) % 8
            }
        }
        if (n > ULong.one) factors.add(n)
        return factors
    }

    // Private worker method (recursive) to obtain the prime factors of a number (ULong).
    static primeFactors_(m, trialDivs) {
        if (m.isPrime) 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.isDivisibleBy(p)) {
                    factors.add(ULong.new(p))
                    n = n / p
                }
            }
        }
        while (n > ULong.one) {
            if (checkPrime && n.isPrime) {
                factors.add(n)
                break
            }
            if (n >= threshold) {
                var d = pollardRho_(n, seed, c)
                if (d != ULong.zero) {
                    factors.addAll(primeFactors_(d, false))
                    n = n / d
                    checkPrime = true
                } else if (c == 1) {
                    if (n.isSquare) {
                        n = n.isqrt
                        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 (ULong) of 'm' (a ULong or an unsigned small integer)
    // 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 = ULong.new(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 = ULong.new(n)
        var divs = []
        var divs2 = []
        var i = one
        var k = (n.isEven) ? one : two
        var sqrt = n.isqrt
        while (i <= sqrt) {
            if (n.isDivisibleBy(i)) {
                divs.add(i)
                var j = n / i
                if (j != i) divs2.add(j)
            }
            i = i + 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 == ULong.zero) return []
        var factors = ULong.primeFactors(n)
        var divs = [ULong.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]
    }
}

/*  ULongs contains various routines applicable to lists of unsigned 64-bit integers. */
class ULongs {
    static sum(a)  { a.reduce(ULong.zero) { |acc, x| acc + x } }
    static mean(a) { sum(a)/a.count }
    static prod(a) { a.reduce(ULong.one) { |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 } }
}

/*
    Long represents a 64-bit signed integer together with arithmetic operations thereon.
    Long objects, which are immutable, are stored as a magnitude (ULong) and a sign
    +1 for positive numbers, 0 for zero and -1 for negative numbers.
*/
class Long is Comparable {
    // Constants
    static minusOne { sigma_(-1, ULong.one)   }
    static zero     { sigma_( 0, ULong.zero)  }
    static one      { sigma_( 1, ULong.one)   }
    static two      { sigma_( 1, ULong.two)   }
    static three    { sigma_( 1, ULong.three) }
    static four     { sigma_( 1, ULong.four)  }
    static five     { sigma_( 1, ULong.five)  }
    static ten      { sigma_( 1, ULong.ten)   }

    // Returns the maximum 'short' Long = 2^32-1 = 4294967295
    static maxShort { sigma_(1, ULong.maxShort) }

    // Returns the minimum 'short' Long = -(2^32-1) = -4294967295
    static minShort { sigma_(-1, ULong.maxShort) }

    // Returns the maximum 'small' Long = 2^53-1 = 9007199254740991
    static maxSmall { sigma_(1, ULong.maxSmall) }

    // Returns the minimum 'small' Long = -(2^53-1) = -9007199254740991
    static minSmall { sigma_(-1, ULong.maxSmall) }

    // Returns the largest Long = 2^64-1 = 18446744073709551615
    static largest  { sigma_(1, ULong.largest) }

    // Returns the smallest Long = -(2^64-1) = -18446744073709551615
    static smallest { sigma_(-1, ULong.smallest) }

    // Private method to determine whether a number is a short integer or not.
    static isShort_(n) { (n is Num) && n.isInteger && n.abs <= 4294967295 }

    // Private method to determine whether a number is a small integer or not.
    static isSmall_(n) { (n is Num) && n.isInteger && n.abs <= 9007199254740991 }

    // Returns the greater of two Longs.
    static max(a, b) {
        if (!(a is Long)) a = new(a)
        if (!(b is Long)) b = new(b)
        return (a > b) ? a : b
    }

    // Returns the lesser of two Longs.
    static min(a, b) {
        if (!(a is Long)) a = new(a)
        if (!(b is Long)) b = new(b)
        return (a < b) ? a : b
    }

    // Returns the greatest common divisor of a and b. The result is never negative.
    static gcd(a, b) {
        if (!(a is Long)) a = new(a)
        if (!(b is Long)) b = new(b)
        while (!b.isZero) {
            var t = b
            b = a % b
            a = t
        }
        return a.abs
    }

    // Returns the least common multiple of a and b. The result is never negative.
    static lcm(a, b) {
        if (!(a is Long)) a = new(a)
        if (!(b is Long)) b = new(b)
        if (a.isZero && b.isZero) return Long.zero
        return (a * b).abs / gcd(a, b)
    }

    // Returns whether or not 'n' is an instance of Long.
    static isInstance(n) { n is Long }

    // Private constructor which creates a Long object from a sign and a ULong.
    construct sigma_(sig, mag) {
        _sig = sig
        _mag = mag
        if (mag.isZero && sig != 0) _sig = 0
    }

    // Private constructor which creates a Long object from a 'small' integer.
    construct fromSmall_(v) {
        _sig = v.sign
        _mag = ULong.fromSmall_(v.abs)
    }

    // Private method which creates a Long object from a Num.
    // If 'v' is not small, will probably lose accuracy.
    static fromNum_(v) { sigma_(v.sign, ULong.fromNum_(v.abs)) }

    // Private method which creates a Long object from a base 10 numeric string.
    // A leading sign is permitted and so is scientific notation.
    // Raises an error if the result is out of bounds.
    static fromString_(v) {
        v = v.trim()
        if (v.count == 0) Fiber.abort("Invalid integer.")
        var sig
        var mag
        if (v[0] == "-") {
            sig = -1
            mag = ULong.fromString_(v[1..-1])
        } else {
            sig = 1
            mag = ULong.fromString_(v)
        }
        if (mag == ULong.zero) sig = 0
        return sigma_(sig, mag)
    }

    // Creates a Long object from an (unprefixed) numeric string in a given base (2 to 36).
    // A leading sign is permitted but scientific notation is not.
    // Wraps out of range values.
    static fromBaseString(v, base) {
        v = v.trim()
        if (v.count == 0) Fiber.abort("Invalid integer.")
        var sig
        var mag
        if (v[0] == "-") {
            sig = -1
            mag = ULong.fromBaseString(v[1..-1], base)
        } else {
            sig = 1
            mag = ULong.fromBaseString(v, base)
        }
        if (mag == ULong.zero) sig = 0
        return sigma_(sig, mag)
    }

    // Creates a Long object from either a numeric base 10 string, a 'small' integer or a ULong.
    static new(value) {
        if (!(value is String) && !isSmall_(value) && !(value is ULong)) {
            Fiber.abort("Value must be a base 10 numeric string, a small integer or a ULong.")
        }
        return (value is String) ? fromString_(value) : (value is Num) ? fromSmall_(value) : value.toLong
    }

    // Creates a Long object from an (unprefixed) hexadecimal string. Wraps out of range values. 
    static fromHexString(v) { fromBaseString(v, 16) }

    // Properties to return the sign and magnitude of this instance.
    sign      { _sig }
    magnitude { _mag }
    mag       { _mag } // synonym for magnitude

    // Public self-evident properties.
    isShort    { _mag.isShort }
    isSmall    { _mag.isSmall }
    isEven     { _mag.isEven  }
    isOdd      { _mag.isOdd   }
    isOne      { _mag.isOne && _sig == 1 }
    isUnit     { _mag.isOne   }
    isZero     { _sig == 0    }
    isPositive { _sig == 1    }
    isNegative { _sig == -1   }

    // Returns true if 'n' is a divisor of the current instance, false otherwise
    isDivisibleBy(n) { (this % n).isZero }

    // negates 'this'
    - { Long.sigma_(-_sig, _mag) }

    // Adds a Long to the current instance. Wraps on overflow.
    + (n) {
        if (!(n is Long)) n = Long.new(n)
        var res
        if (_sig >= 0 && n.sign >= 0) {
            res = Long.sigma_(1, _mag + n.mag)
        } else if (_sig >= 0 && n.sign < 0) {
            res = (_mag >= n.mag) ? Long.sigma_(1, _mag - n.mag) : Long.sigma_(-1, n.mag - _mag)  
        } else if (_sig < 0 && n.sign >= 0) {
            res = (n.mag >= _mag) ? Long.sigma_(1, n.mag - _mag) : Long.sigma_(-1, _mag - n.mag)
        } else  { // _sig < 0 && n.sign < 0)
            res = Long.sigma_(-1, _mag + n.mag)
        }
        return res
    }

    // Subtracts a Long from the current instance. Wraps on underflow.
    - (n) { this + (-n) }

    // Multiplies the current instance by a Long. Wraps on overflow.
    * (n) {
        if (!(n is Long)) n = Long.new(n)
        return Long.sigma_(_sig * n.sign, _mag * n.mag)
    }

    // Returns a list containing the quotient and the remainder after dividing
    // the current instance by a Long. The remainder always has the same sign as 'this'.
    divMod(n) {
        if (!(n is Long)) n = Long.new(n)
        if (n.isZero) Fiber.abort("Cannot divide by zero.")
        var dm = _mag.divMod(n.mag)
        return [Long.sigma_(_sig / n.sign, dm[0]), Long.sigma_(_sig, dm[1])]
    }

    // Divides the current instance by a Long.
    / (n) { divMod(n)[0] }

    // Returns the remainder after dividing the current instance by a Long.
    % (n) { divMod(n)[1] }

    // Returns the absolute value of 'this'.
    abs { isNegative ? -this : this }

    // The inherited 'clone' method just returns 'this' as Long objects are immutable.
    // If you need an actual copy use this method instead.
    copy() { Long.sigma_(_sig, _mag) }

    // Compares the current instance with a Long. If they are equal returns 0.
    // If 'this' is greater, returns 1. Otherwise returns -1.
    // Also allows a comparison with positive or negative infinity.
    compare(n) {
        if ((n is Num) && n.isInfinity && n > 0) return -1
        if ((n is Num) && n.isInfinity && n < 0) return 1
        if (!(n is Long)) n = Long.new(n)
        if (_sig == n.sign && _mag == n.mag) return 0
        if (_sig >= 0 && n.sign >= 0) return (_mag > n.mag) ? 1 : -1
        if (_sig >= 0 && n.sign < 0) return 1
        if (_sig < 0 && n.sign >= 0) return -1
        return (n.mag > _mag) ? 1 : -1
    }

    // Returns the greater of this instance and another Long instance.
    max(n) { Long.max(this, n) }

    // Returns the smaller of this instance and another Long instance.
    min(n) { Long.min(this, n) }

    // Clamps this instance into the range [a, b].
    // If this instance is less than min, min is returned.
    // If it's more than max, max is returned. Otherwise, this instance is returned.
    clamp(a, b) {
        if (!(a is Long)) a = Long.new(a)
        if (!(b is Long)) b = Long.new(b)
        if (a > b) Fiber.abort("Range cannot be decreasing.")
        if (this < a) return a
        if (this > b) return b
        return this.copy()
    }

    // Squares the current instance. Wraps on overflow.
    square { this * this }

    // Cubes the current instance. Wraps on overflow.
    cube { this * this * this }

    // Returns true if the current instance is a perfect square, false otherwise.
    isSquare {
        if (isNegative) System.print("A negative real number cannot be a perfect square.")
        return _mag.isSquare
    }

    // Returns true if the current instance is a perfect cube, false otherwise.
    isCube { _mag.isCube }

    // Returns the integer n'th root of the current instance
    // i.e. the precise n'th root truncated towards zero if not an integer.
    iroot(n) {
        if (!((n is Num) && n.isInteger && n > 0)) {
            Fiber.abort("Argument must be a positive integer.")
        }
        if (isNegative && n.isEven) Fiber.abort("A negative real number cannot have an even root.")
        return Long.sigma_(_sig, _mag.iroot(n))
    }

    // Returns the integer cube root of the current instance
    // i.e. the precise cube root truncated towards zero if not an integer.
    icbrt { Long.sigma_(_sig, _mag.icbrt) }

    // Returns the integer square root of the current instance
    // i.e. the precise sqaure root truncated towards zero if not an integer.
    isqrt {
        if (isNegative) System.print("A negative real number cannot have a square root.")
        return Long.sigma_(_sig, _mag.isqrt)
    }

    // Returns the current instance raised to the power of a 'small' ULong. Wraps on overflow.
    // If the exponent is less than 0, returns 0. O.pow(0) returns one.
    pow(n) {
        var mag = _mag.pow(n)
        var sign 
        if (mag == ULong.zero) {
            sign = 0
        } else if (isNegative && n % 2 == 1) {
            sign = -1
        } else {
            sign = 1
        }
        return Long.sigma_(sign, mag)
    }

    // Returns the current instance multiplied by 'n' modulo 'mod'.
    modMul(n, mod) { Long.sigma_(_sig, _mag.modMul(n, mod)) }

    // Returns the current instance to the power 'exp' modulo 'mod'.
    modPow(exp, mod) {
        var mag = _mag.modPow(exp, mod)
        var sign
        if (mag == ULong.zero) {
            sign = 0
        } else if (isNegative && (exp % mod).isOdd) {
            sign = -1
        } else {
            sign = 1
        }
        return Long.sigma_(sign, mag)
    }

    // Increments the current instance by one.
    inc { this + Long.one }

    // Decrements the current instance by one.
    dec { this - Long.one }

    // Converts the current instance to a Num where possible.
    // Will probably lose accuracy if the current instance is not 'small'.
    toNum { _sig * _mag.toNum }

    // Converts the current instance to a 'small' integer where possible.
    // Otherwise returns null.
    toSmall { isSmall ? toNum : null }

    // Converts the current instance to a ULong where possible.
    // Otherwise returns null.
    toULong { !isNegative ? _mag : null }

    // Returns the string representation of the current instance in a given base (2 to 36).
    toBaseString(base) {
        var bs = _mag.toBaseString(base)
        return isNegative ? "-" + bs : bs
    }

    // Returns the string representation of the current instance in base 16.
    toHexString { toBaseString(16) }

    // Returns the string representation of the current instance in base 10.
    toString { toBaseString(10) }
}

/*  Longs contains various routines applicable to lists of signed 64-bit integers. */
class Longs {
    static sum(a)  { a.reduce(Long.zero) { |acc, x| acc + x } }
    static mean(a) { sum(a)/a.count }
    static prod(a) { a.reduce(Long.one) { |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 } }
}