dickson_linear_forms_prime_sieve_in_range_3.sf

#!/usr/bin/ruby

# Sieve for linear forms primes of the form `a_1*m + b_1`, `a_2*m + b_2`, ..., `a_k*m + b_k`.

# See also:
#   https://oeis.org/A088250
#   https://oeis.org/A318646
#   https://oeis.org/A372238/a372238.gp.txt
#   https://en.wikipedia.org/wiki/Dickson%27s_conjecture

func combine_crt(arr, M, p, S_p) {

    var Minv = invmod(M % p, p)

    gather {
        arr.each {|r|
            var r_mod_p = (r % p)
            S_p.each {|s|
                take(((((s - r_mod_p) % p) * Minv) % p)*M + r)
            }
        }
    }
}

func remainders_mod_p(p, terms) {
    var bad = []  # bad[m] = true means m is forbidden modulo p

    for n,k in (terms) {

        n %= p
        k %= p

        if (n == 0) {
            # Term is constant mod p
            if (k == 0) {
                # Always 0 mod p -> no admissible residue exists
                return []
            }
            next  # This term forbids no residue for this p
        }

        # Forbid the unique residue m ≡ -k * n^{-1} (mod p)
        var n_inv = invmod(n, p)
        var m_forbid = ((-k * n_inv) % p)
        bad[m_forbid] = true
    }

    ^p -> grep { !bad[_] }
}

func select_optimal_primes (A, B, terms) {

    var range = (B - A + 1)

    return [] if (range <= 0)

    var target_modulus = (1+range.iroot(5))**4

    var M = 1
    var primes = []

    for (var p = 2; M <= target_modulus; p.next_prime!) {
        var S_p = remainders_mod_p(p, terms)
        if (S_p.len == p) {
            next    # skip trivial primes
        }
        primes << [p, S_p]
        M *= p
    }

    return primes
}

func linear_form_primes_in_range(A, B, terms) {

    A > B && return []

    var primes = select_optimal_primes(A, B, terms)

    return [] if primes.is_empty

    var residues = [0]
    var M        = 1

    primes.each_2d {|p, S_p|
        residues = combine_crt(residues, M, p, S_p)
        M *= p
    }

    residues.isort!

    return [] if residues.is_empty

    var d = residues.diffs

    while (!d.is_empty && (d[0] == 0)) {
        d.shift
    }

    d << (residues[0] + M - residues[-1])

    var compute_small_values = false
    var small_values_limit   = min(500, B)
    var original_A           = nil

    if (A < small_values_limit) {
        original_A           = A
        A                    = small_values_limit+1
        compute_small_values = true
    }

    var m      = residues[0]
    var d_len  = d.len

    m += M*idiv(A, M)

    var j = 0

    while (m < A) {
        m += d[j++ % d_len]
    }

    var arr = []

    if (compute_small_values) {
        for v in (original_A .. small_values_limit) {
            if (terms.all {|k| is_prime(k[0]*v + k[1]) }) {
                arr << v
            }
        }
    }

    while (m <= B) {

        if (terms.all {|k| (k[0]*m + k[1]).is_prime }) {
            arr << m
        }

        m += d[j++ % d_len]
    }

    return arr
}

assert_eq(linear_form_primes_in_range(1, 41, [[1, 41]]), %n[2, 6, 12, 18, 20, 26, 30, 32, 38])
assert_eq(linear_form_primes_in_range(1, 50, [[1, 1]]), %n[1, 2, 4, 6, 10, 12, 16, 18, 22, 28, 30, 36, 40, 42, 46])
assert_eq(linear_form_primes_in_range(1, 100, [[1,1], [2,1]]), %n[1, 2, 6, 18, 30, 36, 78, 96])
assert_eq(linear_form_primes_in_range(1, 1000, [[1,1],[2,1],[3,1]]), %n[2, 6, 36, 210, 270, 306, 330, 336, 600, 726])
assert_eq(linear_form_primes_in_range(1, 10000, [[1,1],[2,1],[3,1],[4,1]]), %n[330, 1530, 3060, 4260, 4950, 6840])
assert_eq(linear_form_primes_in_range(1, 12000,  [[1,1],[2,1],[3,1],[4,1],[5,1]]), %n[10830])
assert_eq(linear_form_primes_in_range(9538620, 9993270,  [[1,1],[2,1],[3,1],[4,1],[5,1]]), %n[9538620, 9780870, 9783060, 9993270])
assert_eq(linear_form_primes_in_range(9538620+1, 9993270,  [[1,1],[2,1],[3,1],[4,1], [5,1]]), %n[9780870, 9783060, 9993270])

assert_eq(linear_form_primes_in_range(1, 1000,  [[1,-1],[2,-1],[3,-1]]), %n[4, 6, 24, 30, 84, 90, 174, 234, 240, 294, 420, 660, 954])
assert_eq(linear_form_primes_in_range(1, 10000, [[1,-1],[2,-1],[3,-1],[4,-1]]), %n[6, 90, 1410, 1890])
assert_eq(linear_form_primes_in_range(1, 500, [[2,-1],[4,-1],[6,-1]]), %n[2, 3, 12, 15, 42, 45, 87, 117, 120, 147, 210, 330, 477])
assert_eq(linear_form_primes_in_range(1, 500, [[2,1],[4,3],[8,7]]), %n[2, 5, 20, 44, 89, 179, 254, 359])
assert_eq(linear_form_primes_in_range(1, 500, [[2,-1],[4,-1],[8,-1]]), %n[3, 6, 21, 45, 90, 180, 255, 360])
assert_eq(linear_form_primes_in_range(1, 500, [[2,-1],[4,-1],[8,-1],[16,-1]]), %n[3, 45, 90, 180, 255])
assert_eq(linear_form_primes_in_range(1, 500, [[17,1],[23,5]]), %n[18, 24, 66, 126, 186, 216, 378, 384, 426])

assert_eq(linear_form_primes_in_range(1, 500, [[17,4],[15,-8],[19, 2]]), %n[5, 9, 11, 65, 75, 105, 125, 159, 191, 221, 231, 291, 341, 369, 419, 461, 471, 479])
assert_eq(linear_form_primes_in_range(1, 500, [[17,4],[15,+8],[19, 2]]), %n[5, 11, 45, 65, 105, 159, 161, 189, 221, 275, 291, 299, 431, 479])

func f(n, multiple = 1, alpha = 1) {

    var terms = (1..n -> map {|k| [multiple * k, alpha] })

    var A = 1
    var B = 2*A

    loop {
        var arr = linear_form_primes_in_range(A, B, terms)

        if (arr) {
            return arr[0]
        }

        A = B+1
        B = 2*A
    }
}

assert_eq({ f(_, 1, +1) } `map` 1..8, [1, 1, 2, 330, 10830, 25410,  512820,  512820]);     # A088250
assert_eq({ f(_, 1, -1) } `map` 1..8, [3, 3, 4, 6,   6,     154770, 2894220, 2894220]);    # A088651
assert_eq({ f(_, 9, +1) } `map` 1..8, [2, 2, 4, 170, 9860,  23450,  56980,   56980]);      # A372238
assert_eq({ f(_, 2, -1) } `map` 1..8, [2, 2, 2, 3,   3,     77385,  1447110, 1447110]);    # A124492
assert_eq({ f(_, 2, +1) } `map` 1..8, [1, 1, 1, 165, 5415,  12705,  256410,  256410]);     # A071576

assert_eq({ f(_, _, +1) } `map` 1..8, [1, 1, 2, 765,  2166, 4235,  73260,  2780085]);
assert_eq({ f(_, _, -1) } `map` 1..8, [3, 2, 2, 3225, 18,   25795, 413460, 7505190]);

assert_eq({ f(_, _, -13) } `map` 1..6, [15, 8,  6,  15,  24, 2800]);
assert_eq({ f(_, _, +13) } `map` 1..6, [4,  12, 10, 90,  18, 40705]);
assert_eq({ f(_, _, -23) } `map` 1..6, [25, 13, 10, 255, 6,  5]);
assert_eq({ f(_, _, +23) } `map` 1..6, [6,  9,  10, 60,  48, 13300]);

assert_eq({ f(_, 1, +23) } `map` 1..6, [6, 18, 30, 210, 240, 79800]);
assert_eq({ f(_, 1, -23) } `map` 1..8, [25, 26, 30, 30, 30, 30, 142380, 1319010]);

assert_eq({ f(_, 1, +101) } `map` 1..6, [2,   6,   96,  180, 3990, 1683990]);
assert_eq({ f(_, 1, -101) } `map` 1..6, [103, 104, 104, 240, 3630, 78540]);

assert_eq(linear_form_primes_in_range(1, 1e3, [[2, 1], [4, 1], [6, 1]]), [1, 3, 18, 105, 135, 153, 165, 168, 300, 363, 585, 618, 648, 765, 828]);    # A124408
assert_eq(linear_form_primes_in_range(1, 1e4, [[2, 1], [4, 1], [6, 1], [8, 1]]),          [165, 765, 1530, 2130, 2475, 3420, 5415, 7695, 9060]);     # A124409
assert_eq(linear_form_primes_in_range(1, 1e5, [[2, 1], [4, 1], [6, 1], [8, 1], [10, 1]]), [5415, 12705, 13020, 44370, 82950, 98280]);                # A124410
assert_eq(linear_form_primes_in_range(1, 1e6, [[2, 1], [4, 1], [6, 1], [8, 1], [10, 1], [12, 1]]), [12705, 13020, 105525, 256410, 966840]);          # A124411

say "=> The least Chernick's \"universal form\" Carmichael number with n prime factors"

for n in (3..9) {

    var terms = [6, 12, 1 .. n-2 -> map{|k| 9 * 2**k }...].map { [_, 1] }

    var A = 1
    var B = 2*A

    loop {

        var arr = linear_form_primes_in_range(A, B, terms).grep {|k| k.valuation(2) >= n-4 }

        if (arr) {
            say "a(#{n}) = #{arr[0]}"
            break
        }

        A = B+1
        B = 2*A
    }
}

__END__
=> The least Chernick's "universal form" Carmichael number with n prime factors
a(3) = 1
a(4) = 1
a(5) = 380
a(6) = 380
a(7) = 780320
a(8) = 950560