experimental-projects/oeis-research/Tim Johannes Ohrtmann/Least Carmichael number that is divisible by the n-th cyclic number A003277(n), or 0 if no such number exists/prog.pl

#!/usr/bin/perl

# Least Carmichael number that is divisible by the n-th cyclic number A003277(n), or 0 if no such number exists
# https://oeis.org/A253595

use 5.020;
use warnings;

use ntheory      qw(:all);
use experimental qw(signatures);

sub divceil ($x, $y) {    # ceil(x/y)
    my $q = divint($x, $y);
    ($q * $y == $x) ? $q : ($q + 1);
}

sub carmichael_from_multiple ($A, $B, $m, $L, $lo, $k, $callback) {

    my $hi = rootint(divint($B, $m), $k);

    if ($lo > $hi) {
        return;
    }

    if ($k == 1) {

        $lo = vecmax($lo, divceil($A, $m));
        $lo > $hi && return;

        my $t = invmod($m, $L) // return;
        $t > $hi && return;
        $t += $L while ($t < $lo);

        for (my $p = $t ; $p <= $hi ; $p += $L) {
            if ($m % $p != 0 and is_prime($p)) {
                my $n = $m * $p;
                if (($n - 1) % ($p - 1) == 0) {
                    $callback->($n);
                }
            }
        }

        return;
    }

    foreach my $p (@{primes($lo, $hi)}) {

        $m % $p == 0 and next;
        gcd($m, $p - 1) == 1 or next;

        __SUB__->($A, $B, $m * $p, lcm($L, $p - 1), $p + 1, $k - 1, $callback);
    }
}

sub carmichael_divisible_by ($m) {

    $m >= 1 or return;
    $m % 2 == 0 and return;
    is_square_free($m) || return;
    gcd($m, euler_phi($m)) == 1 or return;

    my $A = vecmax(561, $m);
    my $B = 2 * $A;

    my $L = vecmax(1, lcm(map { $_ - 1 } factor($m)));

    my @found;

    for (; ;) {
        for my $k ((is_prime($m) ? 2 : 1) .. 1000) {

            my @P;
            for (my $p = 3 ; scalar(@P) < $k ; $p = next_prime($p)) {
                if ($m % $p != 0 and $L % $p != 0) {
                    push @P, $p;
                }
            }

            last if (vecprod(@P, $m) > $B);

            my $callback = sub ($n) {
                push @found, $n;
                $B = vecmin($B, $n);
            };

            carmichael_from_multiple($A, $B, $m, $L, $P[0], $k, $callback);
        }

        last if @found;

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

    vecmin(@found);
}

carmichael_divisible_by(3) == 561             or die;
carmichael_divisible_by(3 * 5) == 62745       or die;
carmichael_divisible_by(7 * 19) == 1729       or die;
carmichael_divisible_by(47 * 89) == 62745     or die;
carmichael_divisible_by(5 * 47 * 89) == 62745 or die;
carmichael_divisible_by(3 * 47 * 89) == 62745 or die;
carmichael_divisible_by(3 * 89) == 62745      or die;

my $count = 1;
foreach my $n (1..1e6) {
    my $r = carmichael_divisible_by($n);
    if (defined($r)) {
        ++$count;
        say "$count $r";
    }
}

#say join(', ', map { carmichael_divisible_by($_) } @{primes(3, 50)});
#say join(', ', map { carmichael_divisible_by($_) } 1 .. 60);

__END__
561, 1105, 1729, 561, 1105, 561, 1729, 6601, 2465, 2821, 29341, 6601, 334153, 62745
561, 561, 1105, 1729, 561, 1105, 62745, 561, 1729, 6601, 2465, 2821, 561, 825265, 29341, 6601, 334153, 62745, 561, 2433601, 74165065