#!/usr/bin/ruby

# Tests for some Number methods.

# Timings:
#   07 March 2023: 2.020s (with MPU)
#   07 March 2023: 2.418s (without MPU)

assert_eq(8 >> 32, 0)
assert_eq(112 >> 32, 0)
assert_eq(208 >> 64, 0)
assert_eq(496 >> 32, 0)
assert_eq(640 >> 32, 0)

assert(!18721.is_carmichael)

assert_eq(5040.remove(2), 315)
assert_eq(5040.neg.remove(2), -315)
assert_eq(5040.neg.remove(-2), -315)
assert_eq(5040.remove(-2), 315)

assert_eq(5040.remove(32), 5040)
assert_eq(5040.remove(24), 210)
assert_eq(5040.remove(29), 5040)
assert_eq(5040.remove(7), 720)
assert_eq(5040.remove(5), 1008)
assert_eq(5040.remove(-5), -1008)
assert_eq(5040.neg.remove(-3), -560)

assert_eq(12!.remove(-3), -1971200)
assert_eq(12!.neg.remove(-3), 1971200)

assert_eq(5040.remove(0), 5040)
assert_eq(5040.remove(1), 5040)
assert_eq(5040.remove(-1), 5040)

assert_eq(1.powerful(-12, 12), 1.powerful(-12, 12))

assert(irand(0, isqrt(2**64 - 1)).is_between(0, 2**32))
assert(irand(0, isqrt(2**128 - 1)).is_between(0, 2**64))

assert_eq(isqrt(2**64 - 1).factor, factor(2**32 - 1))
assert_eq(isqrt(2**128 - 1).factor, factor(2**64 - 1))

assert(43 `divides` 4171)
assert((-43) `divides` 4171)

assert(41 `divides` 4171 -> not)
assert((-41) `divides` 4171 -> not)

assert(4171 `is_div` 43)
assert(4171 `is_div` (-43))

assert(4171 `is_div` 41 -> not)
assert(4171 `is_div` (-41) -> not)

assert((2**128 + 1)*43 `is_div` -43)
assert((2**128 + 1)*43 `is_div` 43)

assert((2**128 + 1)*43 `is_div` -41 -> not)
assert((2**128 + 1)*43 `is_div` 41 -> not)

assert(43 `divides` ((2**128 + 1)*43))
assert((-43) `divides` ((2**128 + 1)*43))

assert(43 `divides` ((2**128 + 1)*41) -> not)
assert((-43) `divides` ((2**128 + 1)*41) -> not)

assert((2**128 + 1) `is_div` 5704689200685129054721)
assert((2**128 + 1) `is_div` (-5704689200685129054721))

assert((2**128 + 1) `is_div` 5704689200685129054723 -> not)
assert(5704689200685129054723 `divides` (2**128 + 1) -> not)

assert(5704689200685129054721 `divides` (2**128 + 1))
assert((-5704689200685129054721) `divides` (2**128 + 1))

assert(42.is_div(6))
assert(!42.is_div(8))

assert(2**32 - 1 -> is_div(2**16 + 1))
assert(2**32 - 1 -> is_div(2**16 + 3) -> not)
assert(2**64 - 1 -> is_div(2**32 + 1))
assert(2**64 - 1 -> is_div(2**32 + 3) -> not)
assert(2**64 - 1 -> is_div(2**32 - 1))
assert(2**128 - 1 -> is_div(2**64 + 1))

assert(2**128 + 1 -> is_div(59649589127497217))
assert(2**128 + 1 -> is_div(59649589127497215) -> not)
assert(2**128 + 1 -> is_div(5704689200685129054721))
assert(2**128 + 1 -> is_div(5704689200685129054723) -> not)

assert_eq(factorial_valuation(3**25 + 1, 2), 847288609427)
assert_eq(factorial_valuation(3**25 - 1, 3), 423644304696)
assert_eq(factorial_valuation(3**25, 127), 6724512769)
assert_eq(factorial_valuation(2**128 + 1, 2**64 + 1), 18446744073709551615)
assert_eq(factorial_valuation(2**128 + 1, 43), 8101961117165201511032728748375433594)

assert_eq(10.by{ .is_div(3) }, [0, 3, 6, 9, 12, 15, 18, 21, 24, 27])
assert_eq({|v| 10.by{ v.is_div(_) } }(20.primorial), 20.primorial.divisors(14))
assert_eq({|v| 10.by{ v.is_div(_) } }(100.primorial), 100.primorial.divisors(14))

assert_eq(10.by{.is_powerful},    2.powerful(10.nth_powerful))
assert_eq(10.by{.is_powerful(3)}, 3.powerful(10.nth_powerful(3)))

assert(all_prime(5))
assert([2].all_prime)
assert([5].all_prime)
assert([2,3,5,7,11].all_prime)

assert_eq(10.by{.is_safe_prime}, %n[5, 7, 11, 23, 47, 59, 83, 107, 167, 179])

assert_eq(8 >> 8, 0)
assert_eq(0 >> 0, 0)
assert_eq(1234 >> 8, 4)

assert_eq(13.smooth_count(1e16), 511985)

assert_eq(squarefree_sum(1e6), 303961062910)
assert_eq(squarefree_sum(1e7), 30396557311887)

assert_eq(squarefree_count(1e9), 607927124)
assert_eq(squarefree_count(1e10), 6079270942)

assert_eq(prime_count(1e5-1), 9592)
assert_eq(prime_count(1e6-1), 78498)

assert_eq(prime_count(1e5.next_prime-1), 9592)
assert_eq(prime_count(1e6.next_prime-1), 78498)

assert_eq(prime_count(1e5.next_prime), 9593)
assert_eq(prime_count(1e6.next_prime), 78499)

assert_eq(
    30.of { .inc.nth_nonpowerfree(1) },
    30.by { .is_nonpowerfree(1) }
)

assert_eq(
    30.of { .inc.nth_nonpowerfree(2) },
    30.by { .is_nonpowerfree(2) }
)

assert_eq(
    30.of { .inc.nth_nonpowerfree(3) },
    30.by { .is_nonpowerfree(3) }
)

assert_eq(
    100.of { .prime_power_sum },
    100.of {|n|
        sum(1..n.ilog2, {|k|
            sum_primes(2, n.iroot(k), k)
        })
    }
)

assert_eq(
    semiprimes(100),
    1..100 -> grep { .is_semiprime },
)

assert_eq(
    semiprimes(50, 100),
    50..100 -> grep { .is_semiprime },
)

assert_eq(
    prime_powers(100),
    1..100 -> grep { .is_prime_power },
)

assert_eq(
    prime_powers(50, 100),
    50..100 -> grep { .is_prime_power },
)

# RangeNumber k-almost primes

assert_eq(
    gather { 10..100 -> each_almost_prime(2, {|k| take(k) }) },
    2.almost_primes(10, 100)
)

assert_eq(
    gather { 10..100 -> each_squarefree_almost_prime(2, {|k| take(k) }) },
    2.squarefree_almost_primes(10, 100)
)

assert_eq(
    gather { 10..100 -> each_almost_prime(3, {|k| take(k) }) },
    3.almost_primes(10, 100)
)

assert_eq(
    gather { 10..100 -> each_omega_prime(3, {|k| take(k) }) },
    3.omega_primes(10, 100)
)

assert_eq(
    gather { 10..1e4 -> each_carmichael(3, {|k| take(k) }) },
    3.carmichael(10, 1e4)
)

assert_eq(
    gather { 10..1e4 -> each_lucas_carmichael(3, {|k| take(k) }) },
    3.lucas_carmichael(10, 1e4)
)

assert_eq(
    gather { 1e6..1e7 -> each_carmichael(5, {|k| take(k) }) },
    5.carmichael(1e6, 1e7)
)

assert_eq(
    gather { 1e6..1e7 -> each_lucas_carmichael(5, {|k| take(k) }) },
    5.lucas_carmichael(1e6, 1e7)
)

assert_eq(
    gather { 10..100 -> each_squarefree_almost_prime(3, {|k| take(k) }) },
    3.squarefree_almost_primes(10, 100)
)

assert_eq(
    gather { 10..100 `by` 3 -> each_almost_prime(3, {|k| take(k) }) },
    10..100 `by` 3 -> grep { .is_almost_prime(3) }
)

assert_eq(
    gather { 10..100 `by` 3 -> each_almost_prime(3, {|k| take(k) }) },
    10..100 `by` 3 -> grep { .is_almost_prime(3) }
)

assert_eq(
    gather { 3..2000 `by` 2 -> each_carmichael(3, {|k| take(k) }) },
    3..2000 `by` 2 -> grep { .is_carmichael && .is_almost_prime(3) }
)

assert_eq(
    gather { 3..1000 `by` 2 -> each_lucas_carmichael(3, {|k| take(k) }) },
    3..1000 `by` 2 -> grep { .is_lucas_carmichael && .is_almost_prime(3) }
)

assert_eq(
    gather { 10..200 `by` 3 -> each_squarefree_almost_prime(3, {|k| take(k) }) },
    10..200 `by` 3 -> grep { .is_almost_prime(3) && .is_squarefree }
)

assert_eq(
    gather { 10..100 `by` 3 -> each_almost_prime(2, {|k| take(k) }) },
    10..100 `by` 3 -> grep { .is_almost_prime(2) }
)

# RangeNumber k-powerful numbers

assert_eq(
    gather { 10..1000 -> each_powerful(2, {|k| take(k) }) },
    2.powerful(10, 1000)
)

assert_eq(
    gather { 10..1000 -> each_powerful(3, {|k| take(k) }) },
    3.powerful(10, 1000)
)

assert_eq(
    gather { 10..1000 `by` 3 -> each_powerful(3, {|k| take(k) }) },
    10..1000 `by` 3 -> grep { .is_powerful(3) }
)

assert_eq(
    gather { 10..1000 `by` 3 -> each_powerful(2, {|k| take(k) }) },
    10..1000 `by` 3 -> grep { .is_powerful(2) }
)

# RangeNumber squarefree numbers

assert_eq(
    gather { 10..100 -> each_squarefree {|k| take(k) } },
    squarefree(10, 100)
)

assert_eq(
    gather { 10..100 `by` 3 -> each_squarefree {|k| take(k) } },
    10..100 `by` 3 -> grep { .is_squarefree }
)

# RangeNumber prime numbers

assert_eq(
    gather { 10..100 -> each_prime {|k| take(k) } },
    primes(10, 100)
)

assert_eq(
    gather { 1..100 -> each_prime {|k| take(k) } },
    primes(100)
)

assert_eq(
    gather { 10..100 `by` 3 -> each_prime {|k| take(k) } },
    10..100 `by` 3 -> grep { .is_prime }
)

# RangeNumber composite numbers

assert_eq(
    gather { 10..100 -> each_composite {|k| take(k) } },
    composites(10, 100)
)

assert_eq(
    gather { 1..100 -> each_composite {|k| take(k) } },
    composites(100)
)

assert_eq(
    gather { 10..100 `by` 3 -> each_composite {|k| take(k) } },
    10..100 `by` 3 -> grep { .is_composite }
)

# RangeNumber powerfree numbers

assert_eq(
    gather { 10..100 -> each_powerfree(3, {|k| take(k) }) },
    3.powerfree(10, 100)
)

assert_eq(
    gather { 1..100 -> each_powerfree(2, {|k| take(k) }) },
    2.powerfree(100)
)

assert_eq(
    gather { 10..100 `by` 3 -> each_powerfree(3, {|k| take(k) }) },
    10..100 `by` 3 -> grep { .is_cubefree }
)

# RangeNumber non-powerfree numbers

assert_eq(
    gather { 10..100 -> each_nonpowerfree(3, {|k| take(k) }) },
    3.nonpowerfree(10, 100)
)

assert_eq(
    gather { 1..100 -> each_nonpowerfree(2, {|k| take(k) }) },
    2.nonpowerfree(100)
)

assert_eq(
    gather { 10..100 `by` 3 -> each_nonpowerfree(3, {|k| take(k) }) },
    10..100 `by` 3 -> grep { !.is_cubefree }
)

# RangeNumber semiprime numbers

assert_eq(
    gather { 10..100 -> each_semiprime {|k| take(k) } },
    semiprimes(10, 100)
)

assert_eq(
    gather { 10..100 `by` 3 -> each_semiprime {|k| take(k) } },
    10..100 `by` 3 -> grep { .is_semiprime }
)

# RangeNumber prime power numbers

assert_eq(
    gather { 50..100 -> each_prime_power {|k| take(k) } },
    prime_powers(50, 100)
)

assert_eq(
    gather { 10..100 `by` 3 -> each_prime_power {|k| take(k) } },
    10..100 `by` 3 -> grep { .is_prime_power }
)

# Counting methods

for k in (1..4) {

    var n = irand(50, 100)

    assert_eq(k.omega_prime_count(n), 1..n -> count { .is_omega_prime(k) })
    assert_eq(k.almost_prime_count(n), 1..n -> count { .is_almost_prime(k) })
    assert_eq(k.squarefree_almost_prime_count(n), 1..n -> count { .is_squarefree_almost_prime(k) })

    n = irand(100, 1000)

    assert_eq(k.omega_prime_count(n), k.omega_primes(n).len)
    assert_eq(k.almost_prime_count(n), k.almost_primes(n).len)
    assert_eq(k.squarefree_almost_prime_count(n), k.squarefree_almost_primes(n).len)
}

assert_eq(2.powerful(50, 100).len, 2.powerful_count(50, 100))
assert_eq(2.powerful(10, 106).len, 2.powerful_count(10, 106))
assert_eq(2.powerful(50, 106).len, 2.powerful_count(50, 106))
assert_eq(3.powerful(50, 106).len, 3.powerful_count(50, 106))
assert_eq(3.powerful(49, 105).len, 3.powerful_count(49, 105))

assert_eq(
    (faulhaber_sum(2**100 * 10000, 1) - 100.nonpowerfree(2**100 * 10000).sum),
    100.powerfree_sum(2**100 * 10000)
)

assert_eq(2.nonpowerfree(100), 0..100 -> grep { .is_nonpowerfree(2) })
assert_eq(3.nonpowerfree(100), 0..100 -> grep { .is_nonpowerfree(3) })
assert_eq(4.nonpowerfree(100), 0..100 -> grep { .is_nonpowerfree(4) })

assert_eq(2.nonpowerfree(65, 100, 2), %n[68, 72, 75, 76, 80, 81, 84, 88, 90, 92, 96, 98, 99, 100])
assert_eq(3.nonpowerfree(64, 150, 3), %n[64, 72, 80, 81, 88, 96, 104, 108, 112, 120, 125, 128, 135, 136, 144])

assert_eq(2.nonpowerfree_count(65, 100, 2), 14)
assert_eq(3.nonpowerfree_count(64, 150, 2), 15)

assert_eq(2.nonpowerfree_sum(65, 100, 2), 1199)
assert_eq(3.nonpowerfree_sum(64, 150, 2), 1593)

for k in (2..5) {
    var lo = 1e6.irand
    var hi = 1e8.irand
    assert_eq(k.powerful(lo, hi).sum, k.powerful_sum(lo, hi))
}

assert_eq(10.powerful_sum(1e20), 10.powerful(1e20).sum)
assert_eq(11.powerful_sum(1e20), 11.powerful(1e20).sum)

assert_eq(10.powerful_count(1e20), 10.powerful(1e20).len)
assert_eq(11.powerful_count(1e20), 11.powerful(1e20).len)

# Almost prime divisors

assert_eq(5040.almost_prime_divisors.flat.sort, 5040.divisors)
assert_eq(5040.omega_prime_divisors.flat.sort, 5040.divisors)

do {

    var small_k_omega_primes = [
      %n[1],
      %n[2,3,4,5,7,8,9,11,13,16,17,19,23,25,27,29,31,32,37,41,43,47,49,53,59,61,64,67,71,73,79,81,83,89,97,101,103,107,109,113],
      %n[6,10,12,14,15,18,20,21,22,24,26,28,33,34,35,36,38,39,40,44,45,46,48,50,51,52,54,55,56,57,58,62,63,65,68,69,72,74,75,76],
      %n[30,42,60,66,70,78,84,90,102,105,110,114,120,126,130,132,138,140,150,154,156,165,168,170,174,180,182,186,190,195,198,204,220,222,228,230,231,234,238,240],
      %n[210,330,390,420,462,510,546,570,630,660,690,714,770,780,798,840,858,870,910,924,930,966,990,1020,1050,1092,1110,1122,1140,1155,1170,1190,1218,1230,1254,1260,1290,1302,1320,1326],
      %n[2310,2730,3570,3990,4290,4620,4830,5460,5610,6006,6090,6270,6510,6630,6930,7140,7410,7590,7770,7854,7980,8190,8580,8610,8778,8970,9030,9240,9282,9570,9660,9690,9870,10010,10230,10374,10626,10710,10920,11130],
    ]

    small_k_omega_primes.each_kv {|k,v|
        assert_eq(k.omega_primes(v[-1]), v)
        assert_eq(gather { k.omega_primes_each(v[0], v[-1], { take(_) }) }, v)
    }
}

do {

    var (x1, x2) = iquadratic_formula(161022, 770228, -1589768114)

    assert_eq(x1, 97)
    assert_eq(x2, -102)

    var tests = [
        [6, 11, -35],
        [2, -4, -2],
        [-4, -7, 12],
        [20, -15, -10],
        [1, -1, -3],
        [-1, 6, 18],
    ]

    for a,b,c in (tests) {
        var roots = [quadratic_formula(a,b,c)]
        var rootsQ = [quadratic_formulaQ(a,b,c)]

        assert_eq(roots.len, 2)
        assert_eq(rootsQ.len, 2)

        assert(roots.all {|x|
            a*x**2 + b*x + c -> round(-30) =~= 0
        })

        assert(rootsQ.all {|x|
            assert_eq(a*x**2 + b*x + c, 0)
        })
    }
}

do {

    var tests = [
        [1, -6, 11, -6],
        [1, -23, 142, -120],
        [5, 2, -5, -3],
        [1, 4, 6, 5],
        [1, 5, 2, -8],
        [1, 4, 7, 6],
        [-36, 8, -82, 2850986],
        [15, -22, 8, -7520940423059310542039581],
    ]

    for a,b,c,d in (tests) {
        var roots = [cubic_formula(a,b,c,d)]
        assert_eq(roots.len, 3)
        assert(roots.all {|x|
            a*x**3 + b*x**2 + c*x + d -> round(-30) =~= 0
        })
    }
}

do {
    var m = (2**16 + 1)

    # OEIS: A020500
    var x1  = %n[0, 0, 2, 3, 2, 5, 1, 7, 2, 3, 1, 11, 1, 13, 1, 1, 2, 17, 1, 19, 1, 1, 1, 23, 1, 5, 1, 3, 1, 29, 1, 31, 2, 1, 1, 1, 1, 37, 1, 1, 1, 41, 1, 43, 1, 1, 1, 47, 1, 7, 1, 1, 1, 53, 1, 1, 1, 1, 1, 59, 1, 61, 1, 1, 2, 1, 1, 67, 1, 1, 1, 71, 1, 73, 1, 1, 1, 1, 1, 79, 1, 3, 1, 83, 1, 1, 1, 1, 1, 89, 1, 1, 1, 1, 1, 1, 1, 97, 1, 1]

    # OEIS: A020513
    var xn1 = %n[0, -2, 0, 1, 2, 1, 3, 1, 2, 1, 5, 1, 1, 1, 7, 1, 2, 1, 3, 1, 1, 1, 11, 1, 1, 1, 13, 1, 1, 1, 1, 1, 2, 1, 17, 1, 1, 1, 19, 1, 1, 1, 1, 1, 1, 1, 23, 1, 1, 1, 5, 1, 1, 1, 3, 1, 1, 1, 29, 1, 1, 1, 31, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 37, 1, 1, 1, 1, 1, 1, 1, 41, 1, 1, 1, 43, 1, 1, 1, 1, 1, 1, 1, 47, 1, 1, 1, 7, 1]

    assert_eq(100.of {|n| cyclotomic(n, 1) }, x1)
    assert_eq(100.of {|n| cyclotomicmod(n, 1, 29) }, x1.map{_ % 29})

    assert_eq(100.of {|n| cyclotomic(n, -1) }, xn1)
    assert_eq(100.of {|n| cyclotomicmod(n, -1, 29) }, xn1.map {_ % 29})

    for n in (-1 .. 10), x in (-10 .. 10) {
        assert_eq(Str(cyclotomic(n, x) % m), Str(cyclotomicmod(n, x, m)))
    }

    assert_eq(cyclotomicmod(2**127 - 1, 2, 2**127 - 1), 1)
    assert_eq(cyclotomicmod(2**127 - 1, 2*3*5, 2**127 - 1), 1)
    assert_eq(cyclotomicmod(2**127 - 1, 2**127 - 1, 2**127 - 1), 1)

    assert_eq(cyclotomicmod(2**64 + 1, 5040, 503*863), 287421)
    assert_eq(cyclotomicmod(2**64 + 1, -5040, 503*863), cyclotomicmod(2**64 + 1, (-5040)%(503*863), 503*863))

    assert_eq(cyclotomicmod(7!, 5040, 2**128 + 1), 235744758223530422383203401596174491694)
    assert_eq(cyclotomicmod(30!, 5040, 2**128 + 1), 40675970320518606495224484019728682382)

    assert_eq(cyclotomicmod(2**128 - 1, 5040, 2**128 + 1), 258216496104743231120502200247009922076)
    assert_eq(cyclotomicmod(5040*5704689200685129054721, 2**256 + 1, 2**64 + 1), 1)
    assert_eq(cyclotomicmod(5040*5704689200685129054721, next_prime(2**64), prev_prime(2**64)), 385968771661381427)
}

do {
    assert(%n[294409, 167979421, 1152091655881, 62411762908817281, 1516087654274358001].all{.is_imprimitive_carmichael})
    assert(42310088783100741554666880481.is_imprimitive_carmichael)
    assert(21593590390253023722267234622513201.is_imprimitive_carmichael)
    assert(16412975107923138847512341751620644377601.is_imprimitive_carmichael)
    assert(325533792014488126487416882038879701391121.is_imprimitive_carmichael)
    assert(1605045791181700950034233564955898780122791301414374937801.is_imprimitive_carmichael)
}

do {
    with (%n[341, 294409, 9972894583, 1264022137981459, 14054662152215842621]) {|a|
        assert(a.all {.is_super_psp})
        assert(a.none{.is_over_psp})
    }
    assert(192463418472849397730107809253922101.is_super_psp)
    assert(1347320741392600160214289343906212762456021.is_super_psp)
    assert(70865138168006643427403953978871929070133095474701.is_super_psp)
    assert(3363391752747838578311772729701478698952546288306688208857.is_super_psp)
    assert(132153369641266990823936945628293401491197666138621036175881960329.is_super_psp)
    assert(9105096650335639994239038954861714246150666715328403635257215036295306537.is_super_psp)
}

do {
    for n in (%n[2047, 13421773, 14073748835533, 1376414970248942474729, 48663264978548104646392577273, 294413417279041274238472403168164964689, 98117433931341406381352476618801951316878459720486433149, 1252977736815195675988249271013258909221812482895905512953752551821]) {
        assert(n.is_psp)
        assert(n.is_over_psp)
        assert(n.is_strong_psp)
        assert(n.is_super_psp)
        assert(n.is_euler_psp)
        assert(!n.is_pell_psp)
        assert(!n.is_bfw_psp)
        assert(!n.is_bpsw_prp)
        assert(!n.is_fib_psp)
        assert(!n.is_carmichael)
        assert(!n.is_chebyshev_psp)
        assert(!n.is_strong_fib_psp)
        assert(!n.is_lucas_psp)
        assert(!n.is_lucasU_psp(1, -1))
        assert(!n.is_lucasV_psp(1, -1))
        assert(!n.is_lucasU_psp)
        assert(!n.is_lucasV_psp)
        assert(!n.is_strong_lucas_psp)
        assert(!n.is_stronger_lucas_psp)
        assert(!n.is_strongish_lucas_psp)
        assert(!n.is_lucas_carmichael)
        assert(!n.is_imprimitive_carmichael)
    }
}

do {
    assert_eq(20.of {.lnsuperprimorial} ~Z=~= 20.of{.superprimorial.ln} -> uniq, [true])
    assert_eq(20.of {.lnsuperfactorial} ~Z=~= 20.of{.superfactorial.ln} -> uniq, [true])
    assert_eq(20.of {.lnhyperfactorial} ~Z=~= 20.of{.hyperfactorial.ln} -> uniq, [true])
}

do {
    var p = (2**127 - 1)
    var t = nil

    for k in (-100 .. 100) {

        if (1.rand < 0.3) {
            t = next_prime_power(p+k)
            assert(t.is_prime_power)
            assert(t > p+k)

            t = prev_prime_power(p+k)
            assert(t.is_prime_power)
            assert(t < p+k)
        }

        if (1.rand < 0.05) {
            t = next_perfect_power(p+k)
            assert(t.is_perfect_power)
            assert(t > p+k)

            t = prev_perfect_power(p+k)
            assert(t.is_perfect_power)
            assert(t < p+k)
        }
    }
}

do {
    var p = 4171
    var t = nil

    for k in (-100 .. 100) {

        if (1.rand < 0.3) {
            t = next_prime_power(p+k)
            assert(t.is_prime_power)
            assert(t > p+k)

            t = prev_prime_power(p+k)
            assert(t.is_prime_power)
            assert(t < p+k)
        }

        if (1.rand < 0.1) {
            t = next_perfect_power(p+k)
            assert(t.is_perfect_power)
            assert(t > p+k)

            t = prev_perfect_power(p+k)
            assert(t.is_perfect_power)
            assert(t < p+k)
        }
    }
}

do {
    var n = 0
    var arr = []
    20.of {
        arr << (n = next_prime_power(n))
    }
    assert_eq(arr, {.is_prime_power}.first(20))
}

assert_eq(prime_count(619, 619), 1)
assert_eq(prime_count(620, 620), 0)
assert_eq(prime_count(620, 621), 0)
assert_eq(prime_count(619, 621), 1)
assert_eq(prime_count(620, 621), 0)
assert_eq(prime_count(620, 630), 0)
assert_eq(prime_count(620, 631), 1)
assert_eq(prime_count(619, 631), 2)

assert_eq(prime_power_sum(626, 630), prime_powers(626, 630).sum)
assert_eq(prime_power_sum(620, 630), prime_powers(620, 630).sum)
assert_eq(prime_power_sum(620, 631), prime_powers(620, 631).sum)
assert_eq(prime_power_sum(619, 631), prime_powers(619, 631).sum)

assert_eq(prime_power_sum(43**2), prime_powers(43**2).sum)
assert_eq(prime_power_sum(43**3), prime_powers(43**3).sum)

assert_eq(prime_count(2**127 - 1, 2**127 - 1), 1)
assert_eq(3.almost_prime_count(82654, 7158821), 1736362)
assert_eq(4.almost_prime_count(83925, 4376019), 872251)

assert_eq(omega_prime_count(15, 880788066951082290), 10)
assert_eq(omega_prime_count(15, 1922334673882737115), 365)

with (14.omega_primes(20483443417467030)) {|a|
    assert_eq(a.last, 20483443417467030)
    assert_eq(a.len, 15)
}

with (15.omega_primes(1922334673882737115)) {|a|
    assert_eq(a.last, 1920731913647276730)
    assert_eq(a.len, 365)
}

assert_eq(next_squarefree(2**31 + 3), 2**31 + 6)
assert_eq(next_squarefree(2**128 + 1), 2**128 + 2)

assert_eq(prev_squarefree(2**31 + 3), 2**31 + 1)
assert_eq(prev_squarefree(2**128 + 3), 2**128 + 2)

assert_eq(next_powerfree(2**31 + 3, 3), 2**31 + 4)
assert_eq(next_powerfree(2**128 + 1, 3), 2**128 + 2)

assert_eq(prev_powerfree(2**31 + 3, 3), 2**31 + 2)
assert_eq(prev_powerfree(2**128 + 3, 3), 2**128 + 2)

assert_eq(next_semiprime(23), 25)
assert_eq(next_semiprime(4171), 4174)
assert_eq(next_semiprime(4171-1), 4171)
assert_eq(next_semiprime(2**128 + 1), 2**128 + 3)

assert_eq(next_squarefree_semiprime(23), 26)
assert_eq(next_squarefree_semiprime(4171), 4174)
assert_eq(next_squarefree_semiprime(4171-1), 4171)
assert_eq(next_squarefree_semiprime(2**128 + 1), 2**128 + 3)

assert_eq(prev_semiprime(4171), 4169)
assert_eq(prev_semiprime(4171+1), 4171)
assert_eq(prev_semiprime(2**128 + 1), 2**128 - 2)

assert_eq(prev_almost_prime(2**128 + 1, 3), 2**128 - 3)
assert_eq(prev_almost_prime(2**128 + 1, 125), 2**128 - 2**122)
assert_eq(next_almost_prime(2**128 + 1, 126), (75*128**18 + 120*128**17)/15)

#assert_eq(next_omega_prime(16.pn_primorial, 16), 36278497172720993190)  # quite slow

for k in (1..3) {

    var n = 30

    if (k > 1) {
        assert_eq(
            n.inc.of { .next_powerfree(k).prev_powerfree(k) }.grep.uniq,
            k.powerfree(n),
        )
    }

    assert_eq(
        n.inc.of { .next_powerful(k).prev_powerful(k) }.grep.uniq,
        k.powerful(n),
    )

    assert_eq(
        n.inc.of { .next_nonpowerfree(k).prev_nonpowerfree(k) }.grep.uniq,
        k.nonpowerfree(n),
    )

    assert_eq(
        n.inc.of { .next_almost_prime(k).prev_almost_prime(k) }.grep.uniq,
        k.almost_primes(n),
    )

    assert_eq(
        n.inc.of { .next_squarefree_almost_prime(k).prev_squarefree_almost_prime(k) }.grep.uniq,
        k.squarefree_almost_primes(n),
    )

    assert_eq(
        n.inc.of { .next_omega_prime(k).prev_omega_prime(k) }.grep.uniq,
        k.omega_primes(n),
    )
}

assert_eq(next_omega_prime(8.pn_primorial*17, 8), 164949330)
assert_eq(prev_omega_prime(8.pn_primorial*17, 8), 164862390)

assert_eq(next_omega_prime(1e8, 8), 100025310)
assert_eq(prev_omega_prime(1e8, 8), 99981420)

assert_eq(next_squarefree_almost_prime(8.pn_primorial/13 * 23, 8), 17687670)
assert_eq(prev_squarefree_almost_prime(8.pn_primorial/13 * 23, 8), 16546530)

assert_eq(next_squarefree_almost_prime(1e8, 8), 100025310)
assert_eq(prev_squarefree_almost_prime(1e8, 8), 99953490)

var seen = Set()

for n,k in ([[4171, 1], [4171, 2], [4171, 3], [1e4, 4], [1e8, 7], [1e8,6], [1e8,5], [1e12,7], [1e12,6], [1e12,5]]) {

    if (n <= 1e6) {
        assert_eq(prev_almost_prime(next_almost_prime(n, k), k), prev_almost_prime(n+1, k))
        assert_eq(next_almost_prime(prev_almost_prime(n, k), k), next_almost_prime(n-1, k))

        assert_eq(prev_squarefree_almost_prime(next_squarefree_almost_prime(n, k), k), prev_squarefree_almost_prime(n+1, k))
        assert_eq(next_squarefree_almost_prime(prev_squarefree_almost_prime(n, k), k), next_squarefree_almost_prime(n-1, k))

        assert_eq(prev_omega_prime(next_omega_prime(n, k), k), prev_omega_prime(n+1, k))
        assert_eq(next_omega_prime(prev_omega_prime(n, k), k), next_omega_prime(n-1, k))
    }

    if (!seen.has(n)) {     # methods that take only one argument (n)

        seen << n

        assert_eq(next_semiprime(n), ->(n) {
            ++n; while (!n.is_semiprime) { ++n }; n
        }(n))

        assert_eq(prev_semiprime(n), ->(n) {
            --n; while (!n.is_semiprime) { --n }; n
        }(n))

        assert_eq(next_squarefree_semiprime(n), ->(n) {
            ++n; while (!n.is_squarefree_semiprime) { ++n }; n
        }(n))

        assert_eq(prev_squarefree_semiprime(n), ->(n) {
            --n; while (!n.is_squarefree_semiprime) { --n }; n
        }(n))

        assert_eq(next_squarefree(n), ->(n) {
            ++n; while (!n.is_squarefree) { ++n }; n
        }(n))

        assert_eq(prev_squarefree(n), ->(n) {
            --n; while (!n.is_squarefree) { --n }; n
        }(n))

        assert_eq(next_prime(n), ->(n) {
            ++n; while (!n.is_prime) { ++n }; n
        }(n))

        assert_eq(prev_prime(n), ->(n) {
            --n; while (!n.is_prime) { --n }; n
        }(n))

        assert_eq(next_prime_power(n), ->(n) {
            ++n; while (!n.is_prime_power) { ++n }; n
        }(n))

        assert_eq(prev_prime_power(n), ->(n) {
            --n; while (!n.is_prime_power) { --n }; n
        }(n))
    }

    if (n <= 1e6) {
        with (next_omega_prime(n,k)) {|w|
            assert_eq([w], k.omega_primes(n+1, w))
        }

        with (prev_omega_prime(n,k)) {|w|
            assert_eq([w], k.omega_primes(w, n-1))
        }

        with (next_squarefree_almost_prime(n,k)) {|w|
            assert_eq([w], k.squarefree_almost_primes(n+1, w))
        }

        with (prev_squarefree_almost_prime(n,k)) {|w|
            assert_eq([w], k.squarefree_almost_primes(w, n-1))
        }
    }

    with (next_almost_prime(n,k)) {|w|
        assert_eq([w], k.almost_primes(n+1, w))
    }

    with (prev_almost_prime(n,k)) {|w|
        assert_eq([w], k.almost_primes(w, n-1))
    }

    assert_eq(next_omega_prime(n,k), ->(n) {
        ++n; while (!n.is_omega_prime(k)) { ++n }; n
    }(n))

    assert_eq(prev_omega_prime(n,k), ->(n) {
        --n; while (!n.is_omega_prime(k)) { --n }; n
    }(n))

    assert_eq(next_almost_prime(n,k), ->(n) {
        ++n; while (!n.is_almost_prime(k)) { ++n }; n
    }(n))

    assert_eq(prev_almost_prime(n,k), ->(n) {
        --n; while (!n.is_almost_prime(k)) { --n }; n
    }(n))

    assert_eq(next_squarefree_almost_prime(n,k), ->(n) {
        ++n; while (!n.is_squarefree_almost_prime(k)) { ++n }; n
    }(n))

    assert_eq(prev_squarefree_almost_prime(n,k), ->(n) {
        --n; while (!n.is_squarefree_almost_prime(k)) { --n }; n
    }(n))

    if (k > 1) {
        assert_eq(next_powerfree(n,k), ->(n) {
            ++n; while (!n.is_powerfree(k)) { ++n }; n
        }(n))

        assert_eq(prev_powerfree(n,k), ->(n) {
            --n; while (!n.is_powerfree(k)) { --n }; n
        }(n))

        with (next_powerfree(n,k)) {|w|
            assert_eq([w], k.powerfree(n+1, w))
        }

        with (prev_powerfree(n,k)) {|w|
            assert_eq([w], k.powerfree(w, n-1))
        }
    }

    if (n <= 1e6) {
        assert_eq(next_powerful(n,k), ->(n) {
            ++n; while (!n.is_powerful(k)) { ++n }; n
        }(n))

        assert_eq(prev_powerful(n,k), ->(n) {
            --n; while (!n.is_powerful(k)) { --n }; n
        }(n))

        assert_eq(next_nonpowerfree(n,k), ->(n) {
            ++n; while (!n.is_nonpowerfree(k)) { ++n }; n
        }(n))

        assert_eq(prev_nonpowerfree(n,k), ->(n) {
            --n; while (!n.is_nonpowerfree(k)) { --n }; n
        }(n))
    }

    with (next_powerful(n,k)) {|w|
        assert_eq([w], k.powerful(n+1, w))
    }

    with (prev_powerful(n,k)) {|w|
        assert_eq([w], k.powerful(w, n-1))
    }

    with (next_nonpowerfree(n,k)) {|w|
        assert_eq([w], k.nonpowerfree(n+1, w))
    }

    with (prev_nonpowerfree(n,k)) {|w|
        assert_eq([w], k.nonpowerfree(w, n-1))
    }
}

assert_eq(next_perfect_power(4171), do {
    var n = 4171; while (!n.is_perfect_power) { ++n }; n
})

assert_eq(prev_perfect_power(4171), do {
    var n = 4171; while (!n.is_perfect_power) { --n }; n
})

assert_eq(
    next_perfect_power(1<<3320).prev_perfect_power.prev_perfect_power.next_perfect_power,
    1 << 3320,
)

assert_eq(
    10.by { .is_squarefree_almost_prime(2) },
    10.of { .inc.nth_squarefree_almost_prime(2) },
)

assert_eq(
    10.by { .is_squarefree_almost_prime(3) },
    10.of { .inc.nth_squarefree_almost_prime(3) },
)

assert(is_squarefree_almost_prime(2**64 + 1, 2))
assert(!is_squarefree_almost_prime(2**64 + 1, 3))
assert(is_squarefree_almost_prime((2**64 + 1) * 2, 3))
assert(!is_squarefree_almost_prime((2**64 + 1) * 4, 3))
assert(is_squarefree_almost_prime((2**64 + 1) * 503 * 863, 4))
assert(is_squarefree_almost_prime((2**64 + 1) * 503 * 863 * 1000000007, 5))
assert(!is_squarefree_almost_prime((2**64 + 1) * 503 * 863 * 1000000007**2, 5))
assert(!is_squarefree_almost_prime((2**64 + 1) * 503 * 863 * 1000000007**2, 6))

assert(30.by { .is_cyclic }, 30.by { .is_coprime(.phi) })
assert(8177568910636879136524885826320973235.is_cyclic)
assert(17277031840122951876810012573270045985.is_cyclic)
assert(2059832906607460252767290568443059994787898033540634712711845135488141590979778401392385.is_cyclic)
assert(2050222770372550554323267720953963601363820698627252818938445687085323309254089582862054255135745.is_cyclic)

with (29830753961909190919091909190919091909187936016513) {|n|
    for k in (nil, 0, 1) {
        var f = n.special_factor(k)
        assert(f.len >= 2)
        assert_eq(f.prod, n)
    }
}

assert(!is_cyclic(5040))
assert(!is_cyclic(2**256 - 1))
assert(!is_cyclic(2**256 - 11))
assert(!is_cyclic(100!))

# OEIS: A376332
assert(1..12 -> map {|k|       489118019 - k! }.all{.is_semiprime})
assert(1..13 -> map {|k|      6233987183 - k! }.all{.is_semiprime})
assert(1..14 -> map {|k|     87199150463 - k! }.all{.is_semiprime})
assert(1..15 -> map {|k|   1310334397523 - k! }.all{.is_semiprime})
assert(1..16 -> map {|k|  20937254735843 - k! }.all{.is_semiprime})
assert(1..17 -> map {|k| 355693511854763 - k! }.all{.is_semiprime})

assert(33110689978317405605970.phi.is_power(14))
assert(1262838284613442604484690.phi.is_power(15))

assert_eq(sum_primes(1e6), 37550402023)
assert_eq(sum_primes(1, 1e6, 2), 24693298341834533)
assert_eq(sum_primes(2, 1e6, 3), 18393235410255348281725)
assert_eq(sum_primes(1e5, 1e6, 4), 14652141801059963090664358346)
assert_eq(sum_primes(1e4, 1e6, 5), 12174298052688955639420343987695627)

assert_eq(
    30.of { sum_primes(_) }
    30.of { .primes.sum }
)

assert_eq(
    30.of { sum_primes(1, _, 2) }
    30.of { .primes.sum{.sqr} }
)

assert_eq(
    30.of { sum_primes(1, _, 0) }
    30.of { .prime_count }
)

assert_eq(sum_primes(1e8-1000, 1e8, 1), primes(1e8-1000, 1e8).sum)
assert_eq(sum_primes(1e8-1000, 1e8, 2), primes(1e8-1000, 1e8).sum{.sqr})
assert_eq(sum_primes(1e8-1000, 1e8, 3), primes(1e8-1000, 1e8).sum{.ipow(3)})

assert_eq(prev_composite(2**127 + 2), 2**127 + 1)
assert_eq(prev_composite(2**127 + 1), 2**127)
assert_eq(prev_composite(2**127), 2**127 - 2)

assert_eq(prime_count(2**65, 2**65 + 10000), 229)
assert_eq(prime_power_count(2**65, 2**65 + 10000), 230)

assert_eq(
    10.by { .diff_of_squares.len == 3 },
    %n[45, 48, 63, 64, 72, 75, 80, 81, 99, 112]
)

assert_eq(
    diff_of_squares(2**64 + 3),
    [%n[19965236662, 19497792975], %n[9223372036854775810, 9223372036854775809]]
)

assert_eq(
    10.by { .sum_of_squares.len == 3 },
    %n[325, 425, 625, 650, 725, 845, 850, 925, 1025, 1250]
)

assert_eq(
    sum_of_squares(2**64 + 9),
    [%n[3, 4294967296], %n[1202590840, 4123168605], %n[1511828491, 4020089388], %n[2576980380, 3435973835]],
)

assert_eq(primitive_part(423, 2), 1)
assert_eq(primitive_part(13, 5040), 13)

assert_eq(
    20.of { .primitive_part(_) },
    %n[1, 1, 2, 3, 2, 5, 1, 7, 2, 3, 1, 11, 1, 13, 1, 1, 2, 17, 1, 19]
)

assert_eq(
    15.of { primitive_part(_+1, { .fib }) },
    %n[1, 1, 2, 3, 5, 4, 13, 7, 17, 11, 89, 6, 233, 29, 61]
)

assert_eq(
    22.of { primitive_part(_+1, { .factorial }) },
    %n[1, 2, 6, 12, 120, 60, 5040, 1680, 60480, 15120, 39916800, 55440, 6227020800, 8648640, 1816214400, 518918400, 355687428096000, 147026880, 121645100408832000, 55870214400, 1689515283456000, 14079294028800]
)

assert_eq(
    %n[1,2,3,4,5,6,7].nth_permutation(10),
    %n[1, 2, 3, 5, 7, 4, 6],
)

assert_eq(
    %n[1,2,3,4,5,6,7].nth_permutation(100),
    %n[1, 2, 7, 3, 6, 4, 5],
)

assert_eq(18446744073709551606 - (-100),18446744073709551706)
assert_eq((-100) - 18446744073709551606, -18446744073709551706)
assert_eq(100 - 18446744073709551606, -18446744073709551506)
assert_eq(99761490592 * -99761490592, -9952355005137704510464)
assert_eq(-99761490592 * 99761490592, -9952355005137704510464)
assert_eq(-99761490592 * -99761490592, 9952355005137704510464)

for k in (2..100, 200..210, 1234..1254, 9000..9010, 30_000..30_010, 100_000..100_010) {
    assert_eq(20.of { .pyramidal(k) }.grep { !.is_pyramidal(k) }, [])
    assert_eq(20.of { 1e6.irand.pyramidal(k) }.grep { !.is_pyramidal(k) }, [])

    assert_eq(20.of { .polygonal(k) }.grep { !.is_polygonal(k) }, [])
    assert_eq(20.of { 1e6.irand.polygonal(k) }.grep { !.is_polygonal(k) }, [])

    assert_eq(20.of { .centered_polygonal(k) }.grep { !.is_centered_polygonal(k) }, [])
    assert_eq(20.of { 1e6.irand.centered_polygonal(k) }.grep { !.is_centered_polygonal(k) }, [])
}

with (%n[20, 140, 405, 2856, 25296, 111720, 25984, 5474000, 237600, 223826688, 3852800, 268565760, 1834725376, 175861400000, 335674368, 2863363937280, 4383831556096, 206015846400, 3400704000, 938209120583680, 2981338216980480, 21463949229465600, 45410367307776, 72056803765911552]) {|arr|
    assert_eq(
        arr.map_kv {|k,v|
            var t = v.pyramidal_root(k+3).round(-40)
            assert(t.is_int)
            pyramidal(t, k+3)
        },
        arr
    )
}

with (%n[1, 4, 10, 460, 9010, 772210, 20120860, 1553569960, 85507715710, 14932196985010, 1033664429333260, 197628216951078460, 21266854897681220860, 7423007155473283614010, 3108276166302017120182510, 851452464506763307285599610, 32749388246772812069108696710]) {|arr|
    assert_eq(
        arr.map {|v|
            var t = v.centered_polygonal_root(3)
            assert(t.is_int)
            centered_polygonal(t, 3)
        },
        arr
    )
}

with (%n[64, 925, 2976, 93457, 866272, 11025, 3036880, 18412718645101, 9283470627432, 201580440699781, 92839099743040, 5236660451226975, 66779973961058176]) {|arr|
    assert_eq(
        arr.map_kv {|k,v|
            var t = v.centered_polygonal_root(k+3)
            assert(t.is_int)
            centered_polygonal(t, k+3)
        },
        arr
    )
}

with (%n[28, 16, 176, 4950, 8910, 1408, 346500, 277992, 7542080, 326656, 544320, 120400000, 145213440, 48549888, 4733575168, 536813568, 2149576704, 3057500160, 938539560960, 1358951178240, 36324805836800, 99956555776, 49212503949312, 118747221196800, 59461613912064, 13749193801728, 7526849672380416, 98516240758210560, 4969489493917696, 78673429816934400, 4467570822566903808, 1013309912383488000]) {|arr|
    assert_eq(
        arr.map_kv {|k,v|
            var t = v.polygonal_root(k+3)
            assert(t.is_int)
            polygonal(t, k+3)
        },
        arr
    )
}

for name,f,g in ([
    [:pyramidal, Num.method(:pyramidal_root), Num.method(:pyramidal)],
    [:polygonal, Num.method(:polygonal_root), Num.method(:polygonal)]
    [:central, Num.method(:centered_polygonal_root), Num.method(:centered_polygonal)]
]) {
    30.times {
        for n in (1e8.irand, 1e30.irand) {
            var k = irand(3, 100)
            var r = f(n, k)
            assert((g(r.floor, k) <= n) && (g(r.ceil, k) >= n), "[#{name}] Testing: #{n} #{k}")
        }
    }
}

assert_eq(
    20.of { .phi },
    %n[0, 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18]
)

assert_eq(
    20.of { .lambda },
    %n[0, 1, 1, 2, 2, 4, 2, 6, 2, 6, 4, 10, 2, 12, 6, 4, 4, 16, 6, 18]
)

assert_eq(
    20.of { .liouville },
    %n[0, 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1]
)

assert_eq(
    20.of { .moebius },
    20.of { .is_squarefree ? ((-1)**.omega) : 0 }
)

assert_eq(
    30.of { .bell },
    %n[1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570, 4213597, 27644437, 190899322, 1382958545, 10480142147, 82864869804, 682076806159, 5832742205057, 51724158235372, 474869816156751, 4506715738447323, 44152005855084346, 445958869294805289, 4638590332229999353, 49631246523618756274, 545717047936059989389, 6160539404599934652455, 71339801938860275191172],
)

assert_eq(
    30.of { .bellmod(4171) },
    %n[1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570, 4213597, 27644437, 190899322, 1382958545, 10480142147, 82864869804, 682076806159, 5832742205057, 51724158235372, 474869816156751, 4506715738447323, 44152005855084346, 445958869294805289, 4638590332229999353, 49631246523618756274, 545717047936059989389, 6160539404599934652455, 71339801938860275191172].map { .mod(4171) },
)

assert_eq(
    20.by { .is_composite },
    20.by { !.is_prime &&  _>=4 }
)

assert_eq(
    20.by { .is_squarefree },
    20.by { .moebius != 0 }
)

assert_eq(lpf(2**64 + 1), 274177)
assert_eq(gpf(2**64 + 1), 67280421310721)

assert_eq(10.of { .catalan }, %n[1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862])
assert_eq(10.of { euler_number(_) }, %n[1, 0, -1, 0, 5, 0, -61, 0, 1385, 0])
assert_eq(10.of { euler_number(2*_) }, %n[1, -1, 5, -61, 1385, -50521, 2702765, -199360981, 19391512145, -2404879675441])

assert_eq(15.of { .lucas },        %n[2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521, 843])
assert_eq(15.of { .fibonacci(2) }, %n[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377])
assert_eq(15.of { .fibonacci(3) }, %n[0, 0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927])
assert_eq(15.of { .fibonacci(4) }, %n[0, 0, 0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208, 401, 773])

assert_eq(1e5.fib(2), Math.linear_rec(2.of(1), 1.of(0)+[1], 1e5))
assert_eq(1e5.fib(3), Math.linear_rec(3.of(1), 2.of(0)+[1], 1e5))
assert_eq(1e5.fib(4), Math.linear_rec(4.of(1), 3.of(0)+[1], 1e5))

assert_eq(15.of {|n| n! }, %n[1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200])
assert_eq(15.of {|n| n!! }, %n[1, 1, 2, 3, 8, 15, 48, 105, 384, 945, 3840, 10395, 46080, 135135, 645120])
assert_eq(15.of {|n| n.multi_factorial(3) }, %n[1, 1, 2, 3, 4, 10, 18, 28, 80, 162, 280, 880, 1944, 3640, 12320])
assert_eq(15.of {|n| n.subfactorial }, %n[1, 0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961, 14684570, 176214841, 2290792932, 32071101049])

assert_eq(prime_sum(50, 100, 2), primes(50, 100).sum { _**2 })
assert_eq(composite_sum(50, 100, 2), composites(50, 100).sum { _**2 })

assert_eq(composite_sum(100, 97), 0)
assert_eq(composite_sum(100, 97), 0)
assert_eq(composite_sum(93, 93), 93)

assert_eq(composite_sum(1e2, 1e3), composites(1e2, 1e3).sum)
assert_eq(composite_sum(1e2.prev_prime, 1e3), composites(1e2.prev_prime, 1e3).sum)
assert_eq(composite_sum(1e2, 1e3.next_prime), composites(1e2, 1e3.next_prime).sum)
assert_eq(composite_sum(1e2.next_prime, 1e3.next_prime), composites(1e2.next_prime, 1e3.next_prime).sum)

assert_eq(20.by { .is_squarefree }, squarefree(20.nth_squarefree))
assert_eq(20.by { .is_cubefree }, 1..20.nth_cubefree -> grep { .is_powerfree(3) })

assert_eq(20.by { .is_nonsquarefree }, nonsquarefree(20.nth_nonsquarefree))
assert_eq(20.by { .is_noncubefree }, 1..20.nth_noncubefree -> grep { .is_nonpowerfree(3) })

assert_eq(20.by { .is_squarefull }, 2.powerful(20.nth_powerful(2)))
assert_eq(20.by { .is_cubefull }, 3.powerful(20.nth_powerful(3)))

assert_eq(1000!.lpf, 2)
assert_eq(1000!.make_coprime(2*3*5*7).lpf, 11)

assert_eq(1000!.gpf, 1000.prev_prime)
assert_eq(gpf(1000! * (2**64 + 1)), 67280421310721)

assert_eq(20.of{ .nonsquarefree_count },      20.of{ .is_nonsquarefree ? 1 : 0 }.acc)
assert_eq(20.of{ .nonsquarefree_sum },        20.of{ .is_nonsquarefree ? _ : 0 }.acc)
assert_eq(20.of{ .noncubefree_count },        20.of{ .is_noncubefree ? 1 : 0 }.acc)
assert_eq(20.of{ .noncubefree_sum },          20.of{ .is_noncubefree ? _ : 0 }.acc)

assert_eq(20.of{ 4.nonpowerfree_count(_) },   20.of{ .is_nonpowerfree(4) ? 1 : 0 }.acc)
assert_eq(20.of{ 4.nonpowerfree_sum(_) },     20.of{ .is_nonpowerfree(4) ? _ : 0 }.acc)

assert_eq(20.of{ 4.powerfree_count(_) },       20.of{ .is_powerfree(4) ? 1 : 0 }.acc)
assert_eq(20.of{ 4.powerfree_sum(_) },         20.of{ .is_powerfree(4) ? _ : 0 }.acc)

assert_eq(20.of { .squarefree_count },  20.of { .is_squarefree ? 1 : 0 }.acc)
assert_eq(20.of { .squarefree_sum },    20.of { .is_squarefree ? _ : 0 }.acc)
assert_eq(20.of { .cubefree_count },    20.of { .is_cubefree ? 1 : 0 }.acc)
assert_eq(20.of { .cubefree_sum },      20.of { .is_cubefree ? _ : 0 }.acc)
assert_eq(20.of { .squarefull_count },  20.of { .is_squarefull ? 1 : 0 }.acc)
assert_eq(20.of { .squarefull_sum },    20.of { .is_squarefull ? _ : 0 }.acc)
assert_eq(20.of { .cubefull_count },    20.of { .is_cubefull ? 1 : 0 }.acc)
assert_eq(20.of { .cubefull_sum },      20.of { .is_cubefull ? _ : 0 }.acc)
assert_eq(20.of { .composite_count },   20.of { .is_composite ? 1 : 0 }.acc)
assert_eq(20.of { .composite_sum },     20.of { .is_composite ? _ : 0 }.acc)
assert_eq(20.of { .mertens },           20.of { .moebius }.acc)
assert_eq(20.of { .liouville },         20.of { .is_zero ? 0 : ((-1)**bigomega(_)) })
assert_eq(20.of { .liouville_sum },     20.of { .liouville }.acc)
assert_eq(20.of { 2.power_count(_) },   20.of { .is_power(2)&&.is_pos ? 1 : 0 }.acc)
assert_eq(20.of { 3.power_count(_) },   20.of { .is_power(3)&&.is_pos ? 1 : 0 }.acc)
assert_eq(20.of { 2.power_sum(_) },     20.of { .is_power(2) ? _ : 0 }.acc)
assert_eq(20.of { 3.power_sum(_) },     20.of { .is_power(3) ? _ : 0 }.acc)
assert_eq(20.of { .pi },                20.of { .is_prime ? 1 : 0 }.acc)
assert_eq(20.of { .prime_sum },         20.of { .is_prime ? _ : 0 }.acc)
assert_eq(20.of { prime_sum(2,_,0) },   20.of { .is_prime ? 1 : 0 }.acc)
assert_eq(20.of { prime_sum(0,_) },     20.of { .is_prime ? _ : 0 }.acc)
assert_eq(20.of { .prime_sum(nil,2) },  20.of { .is_prime ? .sqr : 0 }.acc)
assert_eq(20.of { prime_sum(2,_,3) },   20.of { .is_prime ? .cube : 0 }.acc)
assert_eq(20.of { prime_sum(3,_,3) },   20.of { (.is_prime && _>=3) ? .cube : 0 }.acc)
assert_eq(20.of { .primorial },         20.of { .is_prime ? _ : 1 }.map_reduce{|a,b| a*b })
assert_eq(20.of { .pn_primorial },      20.of { .nth_prime }.map_reduce{|a,b| a*b })
assert_eq(20.of { .lpf_sum },           20.of { _ >= 2 ? .lpf : 0 }.acc)
assert_eq(20.of { .gpf_sum },           20.of { _ >= 2 ? .gpf : 0 }.acc)
assert_eq(20.of { .sigma0 },            20.of { .divisors.len })
assert_eq(20.of { .sigma(0) },          20.of { .divisors.len })
assert_eq(20.of { .sigma },             20.of { .divisors.sum })
assert_eq(20.of { .sigma(2) },          20.of { .divisors.sum { .sqr } })
assert_eq(20.of { .omega },             20.of { .rad.factor.len })
assert_eq(20.of { .omega },             20.of { .factor_exp.len })
assert_eq(20.of { .bigomega },          20.of { .factor.len })
assert_eq(20.of { .gpf },               20.of { .factor.tail \\ _ })
assert_eq(20.of { .lpf },               20.of { .factor.head \\ _ })
assert_eq(20.of { .sopf },              20.of { .rad.factor.sum })
assert_eq(20.of { .sopfr },             20.of { .factor.sum })
assert_eq(20.of { .factorial },         20.of { .is_zero ? 1 : _ }.map_reduce{|a,b| a*b })
assert_eq(20.of { .superfactorial },    20.of { .factorial }.map_reduce{|a,b| a*b })
assert_eq(20.of { .hyperfactorial },    20.of { _**_ }.map_reduce{|a,b| a*b })
assert_eq(20.of { .subfactorial },      20.of {|n| (-1)**n / n! }.acc.map_kv {|k,v| k! * v })

assert_eq(20.of { .perfect_power_count },     20.of { .is_perfect_power && .is_pos ? 1 : 0 }.acc)
assert_eq(20.of { .perfect_power_sum },       20.of { .is_perfect_power ? _ : 0 }.acc)

for k in (0..3) {
    assert_eq(20.of { k.almost_prime_sum(_) },            20.of { .is_almost_prime(k) ? _ : 0 }.acc)
    assert_eq(20.of { k.squarefree_almost_prime_sum(_) }, 20.of { .is_squarefree_almost_prime(k) ? _ : 0 }.acc)
    assert_eq(20.of { k.omega_prime_sum(_) },             20.of { .is_omega_prime(k) ? _ : 0 }.acc)

    assert_eq(20.of { k.almost_prime_count(_) },            20.of { .is_almost_prime(k) ? 1 : 0 }.acc)
    assert_eq(20.of { k.squarefree_almost_prime_count(_) }, 20.of { .is_squarefree_almost_prime(k) ? 1 : 0 }.acc)
    assert_eq(20.of { k.omega_prime_count(_) },             20.of { .is_omega_prime(k) ? 1 : 0 }.acc)

    assert_eq(20.of { .perfect_power_count(k) },  20.of { .is_power(k) && .is_pos ? 1 : 0 }.acc)
    assert_eq(20.of { .perfect_power_sum(k) },    20.of { .is_power(k) ? _ : 0 }.acc)

    assert_eq(20.of { k.powerful_count(_) },  20.of { .is_powerful(k) && .is_pos ? 1 : 0 }.acc)
    assert_eq(20.of { k.powerful_sum(_) },    20.of { .is_powerful(k) ? _ : 0 }.acc)
}

assert_eq(20.of { .digits(2).sum },     20.of { .hammingweight })
assert_eq(20.of { .digits(2).len },     20.of { .len(2) })
assert_eq(20.of { .digits.len },        20.of { .len })
assert_eq(20.of { .digits(10).len },    20.of { .len(10) })

assert_eq(20.of { .sumdigits },         20.of { .digits.sum })
assert_eq(20.of { .sumdigits(2) },      20.of { .digits(2).sum })
assert_eq(20.of { .sumdigits(3) },      20.of { .digits(3).sum })
assert_eq(20.of { .sumdigits(8) },      20.of { .digits(8).sum })
assert_eq(20.of { .sumdigits(9) },      20.of { .digits(9).sum })
assert_eq(20.of { .sumdigits(10) },     20.of { .digits(10).sum })
assert_eq(20.of { .sumdigits(11) },     20.of { .digits(11).sum })
assert_eq(20.of { .sumdigits(16) },     20.of { .digits(16).sum })

assert_eq(20.of { .digits.digits2num },            20.of { _ })
assert_eq(20.of { .digits(2).digits2num(2) },      20.of { _ })
assert_eq(20.of { .digits(3).digits2num(3) },      20.of { _ })
assert_eq(20.of { .digits(8).digits2num(8) },      20.of { _ })
assert_eq(20.of { .digits(9).digits2num(9) },      20.of { _ })
assert_eq(20.of { .digits(10).digits2num(10) },    20.of { _ })
assert_eq(20.of { .digits(11).digits2num(11) },    20.of { _ })
assert_eq(20.of { .digits(16).digits2num(16) },    20.of { _ })

assert_eq(20.of { composite_sum(0, _, 0) },  20.of { .composite_count })
assert_eq(20.of { composite_sum(2, _, 0) },  20.of { .is_composite ? 1 : 0 }.acc)
assert_eq(20.of { composite_sum(3, _, 1) },  20.of { .is_composite ? _ : 0 }.acc)
assert_eq(20.of { .composite_sum(nil, 2) },  20.of { .is_composite ? _**2 : 0 }.acc)
assert_eq(20.of { composite_sum(0, _, 2) },  20.of { .is_composite ? _**2 : 0 }.acc)
assert_eq(20.of { composite_sum(4, _, 3) },  20.of { .is_composite ? _**3 : 0 }.acc)
assert_eq(20.of { composite_sum(6, _, 3) },  20.of { (.is_composite && _>=6) ? _**3 : 0 }.acc)

for n in (2147483646, 2147483647, 2147483648, 4294967293, 4294967294, 4294967295, 4294967297, 9223372036854775808, 18446744073709551613, 18446744073709551614, 18446744073709551617) {

    assert_eq(n.digits.sum, n.sumdigits)
    assert_eq(n.digits.digits2num, n)

    assert_eq(n.digits(2).join.flip, n.as_bin)
    assert_eq(n.digits(8).join.flip, n.as_oct)
    assert_eq(n.digits(16).len, n.as_hex.len)

    var bases = [2, 8, 9, 10, 11, 16, 43, n, irand(2, n), irand(2, n**2), irand(2, n**3), irand(2, n.isqrt), irand(2, n.iroot(3)), 4294967294, 4294967295, 4294967296, 4294967297, 2147483646, 2147483647, 2147483648, 2147483649, 1073741823, 1073741829].uniq
    var nums  = [n, n.irand, irand(n**2), irand(n**3)].uniq

    for b in (bases), k in (nums) {
        var D = k.digits(b)
        assert_eq(D.len, k.len(b))
        assert_eq(D.digits2num(b), k)
        assert_eq(k.sumdigits(b), D.sum)
    }
}

assert_eq(lpf_sum(50, 100), 50..100 -> sum{.lpf})
assert_eq(gpf_sum(50, 100), 50..100 -> sum{.gpf})

assert_eq(convergents(3/2), [1, 3/2])
assert_eq(convergents(2/3), [0, 1, 2/3])
assert_eq(convergents(42), [42])
assert_eq(convergents(-42), [-42])
assert_eq(convergents(2**128 + 1), [2**128 + 1])

assert_eq(42.as_cfrac, [42])
assert_eq(as_cfrac(2**128 + 1), [2**128 + 1])
assert_eq(42.neg.as_cfrac, [-42])
assert_eq((3/4).as_cfrac, [0, 1, 3])
assert_eq((2/3).as_cfrac, [0, 1, 2])

assert_eq(18446744073709551606.next_prime, 18446744073709551629)
assert_eq(18446744073709551629.prev_prime, 18446744073709551557)

assert_eq(euler_numbers(0), %n[1])
assert_eq(euler_numbers(1), %n[1, 0])
assert_eq(euler_numbers(2), %n[1, 0, -1])
assert_eq(euler_numbers(3), %n[1, 0, -1, 0])
assert_eq(euler_numbers(4), %n[1, 0, -1, 0, 5])

assert_eq(bernoulli_numbers(0), %n[1])
assert_eq(bernoulli_numbers(1), %n[1, -1/2])
assert_eq(bernoulli_numbers(2), %n[1, -1/2, 1/6])
assert_eq(bernoulli_numbers(3), %n[1, -1/2, 1/6, 0])
assert_eq(bernoulli_numbers(4), %n[1, -1/2, 1/6, 0, -1/30])

assert_eq( euler_numbers(10),  10.inc.of { .euler_number })
assert_eq( euler_numbers(11),  11.inc.of { .euler_number })
assert_eq(euler_numbers(100), 100.inc.of { .euler_number })
assert_eq(euler_numbers(200), 200.inc.of { .euler_number })

assert_eq( bernoulli_numbers(10),  10.inc.of { .bernoulli(0) })
assert_eq( bernoulli_numbers(11),  11.inc.of { .bernoulli(0) })
assert_eq(bernoulli_numbers(100), 100.inc.of { .bernoulli(0) })
assert_eq(bernoulli_numbers(200), 200.inc.of { .bernoulli(0) })

assert_eq(fubini_numbers(10), %n[1, 1, 3, 13, 75, 541, 4683, 47293, 545835, 7087261, 102247563])
assert_eq(fubini_numbers(5), %n[1, 1, 3, 13, 75, 541])
assert_eq(21.of { .fubini }, %n[1, 1, 3, 13, 75, 541, 4683, 47293, 545835, 7087261, 102247563, 1622632573, 28091567595, 526858348381, 10641342970443, 230283190977853, 5315654681981355, 130370767029135901, 3385534663256845323, 92801587319328411133, 2677687796244384203115])

assert_eq(20.of { .liouville_sum }, %n[0, 1, 0, -1, 0, -1, 0, -1, -2, -1, 0, -1, -2, -3, -2, -1, 0, -1, -2, -3])

assert_eq(liouville_sum(123456), -92)
assert_eq(liouville_sum(12345), -49)

assert_eq(
    gather { primes_each(2**65, 2**65 + 1000, {|k| take(k) }) },
    prime_count(2**65, 2**65 + 1000).next_primes(2**65),
)

assert_eq(
    primes(2**63 - 10, 2**63 + 1000),
    prime_count(2**63 - 10, 2**63 + 1000).next_primes(2**63 - 10),
)

assert_eq(
    primes(2**64 - 10, 2**64 + 1000),
    prime_count(2**64 - 10, 2**64 + 1000).next_primes(2**64 - 10),
)

assert_eq(
    primes(2**65, 2**65 + 1000),
    prime_count(2**65, 2**65 + 1000).next_primes(2**65),
)

assert_eq(
    10.by {|n|
        n.is_composite && (divmod(powmod(10, n-1, n)-1, 9, n) == 0)
    },
    %n[91 259 451 481 703 1729 2821 2981 3367 4141],
)

assert_eq(
    (-8 .. 8 -> map { .is_power }),
    [true, false, false, false, false, false, false, true, true, true, false, false, true, false, false, false, true]
)

assert_eq(
    ((-8 .. 8) ~X (-8 .. 8) -> map_2d {|n,k| is_power(n,k) }),
    [false, false, false, false, false, false, false, false, false, true, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, true, false, true, false, true, false, true, false, true, false, true, false, true, false, true, false, false, false, false, false, false, false, false, false, false, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, false, true, false, false, false, false, false]
)

for n in ([@|(0..20), 120, 480, 5040, (2**64 + 1)**3, 2**64 + 1, 2**127 - 1, (2**127 - 1)**2]) {
    var D = n.divisors
    assert_eq(D.sum{|d| phi(d) }, n)
    assert_eq(D.sum{|d| phi(d,2) }, n**2)
    assert_eq(D.sum{|d| phi(d,3) }, n**3)

    var D = n.idivisors
    assert_eq(D.sum{|d| iphi(d) }, n)
    assert_eq(D.sum{|d| iphi(d,2) }, n**2)
    assert_eq(D.sum{|d| iphi(d,3) }, n**3)

    var D = n.udivisors
    assert_eq(D.sum{|d| uphi(d) }, n)
    assert_eq(D.sum{|d| uphi(d,2) }, n**2)
    assert_eq(D.sum{|d| uphi(d,3) }, n**3)
}

assert_eq(map(1..70, {.nuphi}), %n[1,2,3,2,5,6,7,4,6,10,11,6,13,14,15,8,17,12,19,10,21,22,23,12,20,26,18,14,29,30,31,16,33,34,35,12,37,38,39,20,41,42,43,22,30,46,47,24,42,40,51,26,53,36,55,28,57,58,59,30,61,62,42,32,65,66,67,34,69,70])
assert_eq(map(1..74, {.bphi}), %n[1,1,2,3,4,3,6,7,8,6,10,8,12,9,9,15,16,12,18,14,14,15,22,17,24,18,26,21,28,15,30,31,23,24,25,29,36,27,28,31,40,21,42,35,34,33,46,36,48,36,37,42,52,39,42,46,42,42,58,34,60,45,51,63,50,35,66,56,51,38,70,62,72,54])

assert_eq(znlog(232752345212475230211680, 23847293847923847239847098123812075234, 804842536444911030681947), 13)
assert_eq(znlog(5678, 5, 10007), 8620)

[
 [[5675,5,10000019], 2003974];  # 5675 = 5^2003974 mod 10000019
 [[18478760,5,314138927], 34034873];
 [[553521,459996,557057], 15471];
 [[7443282,4,13524947], 6762454];
 (Bool(Num.HAS_PRIME_UTIL) && (Num(Num.ULONG_MAX) > 2**32) ? [[32712908945642193,5,71245073933756341], 5945146967010377] : ());
].each_2d {|t, v|
    say "Testing: znlog(#{t.join(', ')}) == #{v}"
    assert_eq(znlog(t...), v)
}

# OEIS: https://oeis.org/A368350
assert_eq(
    36.of{|n| znlog(2**n - 5, 3, 2**(n+1)) }.join(' '),
    %n[NaN, 0, NaN, 1, 7, 3, 27, 43, 75, 139, 11, 779, 267, 1291, 3339, 7435, 32011, 48395, 81163, 146699, 277771, 15627, 1588491, 2637067, 539915, 4734219, 13122827, 63454475, 29900043, 231226635, 97008907, 902315275, 365444363, 1439186187, 3586669835, 7881637131].join(' ')
)

# GCD/LCM tests on negative inputs

assert_eq(lcm(-43), 43)
assert_eq(gcd(-43), 43)
assert_eq(gcud(-43), 43)

assert_eq(lcm(-42.7), 42)
assert_eq(gcd(-42.7), 42)
assert_eq(gcud(-42.7), 42)

assert_eq(lcm(-42.7, 1), 42)
assert_eq(gcd(-42.7, 0), 42)
assert_eq(gcud(-42.7, 0), 42)

assert_eq(lcm(-42.7, 1, 1), 42)
assert_eq(gcd(-42.7, 0, 0), 42)
assert_eq(gcud(-42.7, 0, 0), 42)

assert_eq(lcm(1, 1, -42.7), 42)
assert_eq(gcd(0, 0, -42.7), 42)
assert_eq(gcud(0, 0, -42.7), 42)

assert_eq(farey(5), [0, 1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5, 1]);
assert_eq(10.of{.farey.map{.nu}.sum}, %n[0, 1, 2, 5, 9, 19, 25, 46, 62, 89])
assert_eq(10.of{.farey.map{.de}.sum}, %n[1, 2, 4, 10, 18, 38, 50, 92, 124, 178])

assert_eq([farey_neighbors(5, 1/5)], [0, 1/4])
assert_eq([farey_neighbors(5, 1/4)], [1/5, 1/3])
assert_eq([farey_neighbors(5, 1/3)], [1/4, 2/5])
assert_eq([farey_neighbors(5, 2/5)], [1/3, 1/2])
assert_eq([farey_neighbors(5, 1/2)], [2/5, 3/5])
assert_eq([farey_neighbors(5, 3/5)], [1/2, 2/3])
assert_eq([farey_neighbors(5, 2/3)], [3/5, 3/4])
assert_eq([farey_neighbors(5, 3/4)], [2/3, 4/5])
assert_eq([farey_neighbors(5, 4/5)], [3/4, 1])

assert_eq([farey_neighbors(100, 89/99)], [80/89, 9/10])
assert_eq([farey_neighbors(100, 50/79)], [31/49, 19/30])
assert_eq([farey_neighbors(100, 19/30)], [50/79, 45/71])

assert_eq(egypt_greedy(9/10), %n[2, 3, 15])
assert_eq(egypt_greedy(5/121), %n[25, 757, 763309, 873960180913, 1527612795642093418846225])

assert_eq(egypt_greedy(43/97).sum{.inv}, 43/97)
assert_eq(egypt_greedy(97/43).sum{.inv}, 97/43)

say "** Tests passed!"