[feature request] Sieve for primes in linear forms: `a_1*m + b_1`, `a_2*m + b_2`, ..., `a_k*m + b_k`

[trizen]

Would be nice to have an efficient sieve that takes a list of parameters [a_1, b_1], [a_2, b_2], ..., [a_k, b_k] and finds values of m in a given range [A, B], such that a_1*m + b_1, a_2*m + b_2, ..., a_k*m + b_k, are all prime numbers.

This would be a generalization of the sieve_prime_cluster function, such that:

    sieve_prime_cluster($A, $B, @C) == linear_form_primes_in_range($A, $B, [[1, 0], map { [1, $_] } @C]);

A Perl implementation of such a sieve is given below:

use utf8;
use 5.036;
use ntheory qw(:all);
use Test::More tests => 40;

sub remainders_mod_p($p, $terms) {
    my @bad;    # bad[m] = 1 means m is forbidden modulo p

    foreach my $pair (@$terms) {
        my ($n, $k) = @$pair;

        $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)
        my $n_inv    = invmod($n, $p);
        my $m_forbid = (-$k * $n_inv) % $p;
        $bad[$m_forbid] = 1;
    }

    return grep { !$bad[$_] } 0 .. $p - 1;
}

sub combine_crt($arr, $M, $p, $S_p) {

    my @res;
    my $Minv = invmod($M % $p, $p);

    foreach my $r (@$arr) {
        my $r_mod_p = $r % $p;
        foreach my $s (@$S_p) {
            push @res, (((($s - $r_mod_p) % $p) * $Minv) % $p) * $M + $r;
        }
    }

    return \@res;
}

sub remainders_for_primes($primes) {

    my $M        = 1;
    my $residues = [0];

    foreach my $pair (@$primes) {
        my ($p, $S_p) = @$pair;

        # Early return if no valid residues
        return ($M, []) unless @$S_p;

        $residues = combine_crt($residues, $M, $p, $S_p);
        $M *= $p;
    }

    return ($M, [sort { $a <=> $b } @$residues]);
}

sub deltas ($integers) {

    my @deltas;
    my $prev = 0;

    foreach my $n (@$integers) {
        push @deltas, $n - $prev;
        $prev = $n;
    }

    shift(@deltas);
    return \@deltas;
}

sub select_optimal_primes ($A, $B, $terms) {

    my $range = $B - $A + 1;
    return () if $range <= 0;

    my $target_modulus = (1 + rootint($range, 5))**4;

    my $M = 1;
    my @primes;

    for (my $p = 2 ; $M <= $target_modulus ; $p = next_prime($p)) {
        my @S_p = remainders_mod_p($p, $terms);

        if (scalar(@S_p) == $p) {
            next;    # skip trivial primes
        }

        push(@primes, [$p, \@S_p]);
        $M *= $p;
    }

    return @primes;
}

sub linear_form_primes_in_range($A, $B, $terms) {

    return [] if ($A > $B);
    return [] if !@$terms;

    # Select an optimal subset of small primes for the modular sieve
    my @primes = select_optimal_primes($A, $B, $terms);
    return [] unless @primes;

    my ($M, $r) = remainders_for_primes(\@primes);

    # If there are no admissible residues modulo M, there can be no solutions
    return [] if !@$r;

    # Compute differences
    my @d = @{deltas($r)};

    # Remove leading zeros
    while (@d and $d[0] == 0) {
        shift @d;
    }

    # Add wraparound delta
    push @d, $r->[0] + $M - $r->[-1];

    my $compute_small_values = 0;
    my $small_values_limit   = vecmin(500, $B);
    my $original_A           = undef;

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

    my $m      = $r->[0];
    my $d_len  = scalar(@d);
    my $prev_m = $m;

    # Jump near to the start of the range
    $m += $M * divint($A, $M);

    my $j = 0;

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

    my @multiples = map { $_->[0] } @$terms;
    my @alphas    = map { $_->[1] } @$terms;
    my @range     = (0 .. $#multiples);

    my ($ok, @arr);

    # Compute small values if needed
    if ($compute_small_values) {
        foreach my $k ($original_A .. $small_values_limit) {

            $ok = 1;
            foreach my $i (@range) {
                if (!is_prime($multiples[$i] * $k + $alphas[$i])) {
                    $ok = 0;
                    last;
                }
            }

            $ok && push @arr, $k;
        }
    }

    while ($m <= $B) {

        $ok = 1;
        foreach my $i (@range) {
            if (!is_prime($multiples[$i] * $m + $alphas[$i])) {
                $ok = 0;
                last;
            }
        }

        $ok && push(@arr, $m);
        $m += $d[$j++ % $d_len];
    }

    return \@arr;
}

# Several tests and examples

is_deeply(linear_form_primes_in_range(1, 41, [[1, 41]]),                  [2, 6, 12, 18, 20, 26, 30, 32, 38]);
is_deeply(linear_form_primes_in_range(1, 50, [[1, 1]]),                   [1, 2, 4, 6, 10, 12, 16, 18, 22, 28, 30, 36, 40, 42, 46]);
is_deeply(linear_form_primes_in_range(1, 100, [[1, 1], [2, 1]]),          [1, 2, 6, 18, 30, 36, 78, 96]);
is_deeply(linear_form_primes_in_range(1, 1000, [[1, 1], [2, 1], [3, 1]]), [2, 6, 36, 210, 270, 306, 330, 336, 600, 726]);

is_deeply(linear_form_primes_in_range(1, 100, [[5, 5]]), []);
is_deeply(linear_form_primes_in_range(1, 100, []),       []);

is_deeply(linear_form_primes_in_range(270,         270,     [[1, 1], [2, 1], [3, 1]]),                 [270]);
is_deeply(linear_form_primes_in_range(1,           10000,   [[1, 1], [2, 1], [3, 1], [4, 1]]),         [330, 1530, 3060, 4260, 4950, 6840]);
is_deeply(linear_form_primes_in_range(1,           12000,   [[1, 1], [2, 1], [3, 1], [4, 1], [5, 1]]), [10830]);
is_deeply(linear_form_primes_in_range(9538620,     9993270, [[1, 1], [2, 1], [3, 1], [4, 1], [5, 1]]), [9538620, 9780870, 9783060, 9993270]);
is_deeply(linear_form_primes_in_range(9538620,     9538620, [[1, 1], [2, 1], [3, 1], [4, 1], [5, 1]]), [9538620]);
is_deeply(linear_form_primes_in_range(9538620 + 1, 9993270, [[1, 1], [2, 1], [3, 1], [4, 1], [5, 1]]), [9780870, 9783060, 9993270]);

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

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

sub f($n, $multiple = 1, $alpha = 1) {

    my @terms = map { [$multiple * $_, $alpha] } 1 .. $n;

    my $A = 1;
    my $B = 2 * $A;

    while (1) {
        my @arr = @{linear_form_primes_in_range($A, $B, \@terms)};

        if (@arr) {
            return $arr[0];
        }

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

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

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

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

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

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

is_deeply(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
is_deeply(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
is_deeply(linear_form_primes_in_range(1, 1e5, [[2, 1], [4, 1], [6, 1], [8, 1], [10, 1]]), [5415, 12705, 13020, 44370, 82950, 98280]);                # A124410
is_deeply(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 "\n=> The least Chernick's \"universal form\" Carmichael number with n prime factors";

foreach my $n (3 .. 9) {

    my $terms = [map { [$_, 1] } (6, 12, (map { 9 * (1 << $_) } 1 .. $n - 2))];

    my $A = 1;
    my $B = 2 * $A;

    while (1) {
        my @arr = @{linear_form_primes_in_range($A, $B, $terms)};

        @arr = grep { valuation($_, 2) >= $n - 4 } @arr;

        if (@arr) {
            say "a($n) = $arr[0]";
            last;
        }

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

say "\n=> Smallest number k such that r*k + 1 is prime for all r = 1 to n";

foreach my $n (1 .. 10) {
    say "a($n) = ", f($n, 1, 1);
}