123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121 |
- #!/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
|