123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384 |
- #!/usr/bin/perl
- # Daniel "Trizen" Șuteu
- # Date: 27 August 2022
- # https://github.com/trizen
- # Generate all the Lucas-Carmichael numbers with n prime factors in a given range [a,b]. (not in sorted order)
- # See also:
- # https://en.wikipedia.org/wiki/Almost_prime
- use 5.020;
- 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 lucas_carmichael_numbers_in_range ($A, $B, $k, $callback) {
- $A = vecmax($A, pn_primorial($k));
- sub ($m, $lambda, $p, $k, $u = undef, $v = undef) {
- if ($k == 1) {
- printf("# $m -> ($u, $v) -> 10^%.3f\n", log($v - $u) / log(10)) if ($v - $u > 2e6);
- if ($v - $u > 1e10) {
- die "Range too large!\n";
- }
- foreach my $p (@{primes($u, $v)}) {
- my $t = $m * $p;
- if (($t + 1) % $lambda == 0 and ($t + 1) % ($p + 1) == 0) {
- $callback->($t);
- }
- }
- return;
- }
- my $s = rootint(divint($B, $m), $k);
- for (my $r ; $p <= $s ; $p = $r) {
- $r = next_prime($p);
- my $L = lcm($lambda, $p + 1);
- gcd($L, $m) == 1 or next;
- my $t = $m * $p;
- my $u = divceil($A, $t);
- my $v = divint($B, $t);
- if ($u <= $v) {
- __SUB__->($t, $L, $r, $k - 1, (($k == 2 && $r > $u) ? $r : $u), $v);
- }
- }
- }
- ->(1, 1, 3, $k);
- }
- my $min_k = 10; # mininum number of prime factors
- my $max_k = 1e4; # maxinum number of prime factors
- #my $from = 1; # generate terms >= this
- #my $upto = 196467189590024639; # generate terms <= this
- # Generate terms in this range
- my $from = 3614740529402519;
- my $upto = 20576473996736735;
- foreach my $k ($min_k .. $max_k) {
- last if pn_primorial($k) > $upto;
- say "# Generating with k = $k";
- lucas_carmichael_numbers_in_range($from, $upto, $k, sub ($n) { say $n });
- }
|