carmichael_with_n_prime_factors_from_prime_factors.pl 3.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123
  1. #!/usr/bin/perl
  2. # Daniel "Trizen" Șuteu
  3. # Date: 24 September 2022
  4. # https://github.com/trizen
  5. # See also:
  6. # https://en.wikipedia.org/wiki/Almost_prime
  7. # https://trizenx.blogspot.com/2020/08/pseudoprimes-construction-methods-and.html
  8. use 5.020;
  9. use warnings;
  10. use ntheory qw(:all);
  11. use experimental qw(signatures);
  12. use Math::GMPz;
  13. sub divceil ($x, $y) { # ceil(x/y)
  14. my $q = ($x / $y);
  15. ($q * $y == $x) ? $q : ($q + 1);
  16. }
  17. sub carmichael_in_range ($A, $B, $k, $primes, $callback) {
  18. $A = vecmax($A, pn_primorial($k));
  19. $A = Math::GMPz->new("$A");
  20. if ($A > $B) {
  21. return;
  22. }
  23. my $end = $#{$primes};
  24. my $k_exp = 1;
  25. my $congr = -1;
  26. sub ($m, $lambda, $j, $k) {
  27. my $y = rootint(($B / $m), $k);
  28. if ($k == 1) {
  29. my $x = divceil($A, $m);
  30. if ($primes->[-1] < $x) {
  31. return;
  32. }
  33. foreach my $i ($j .. $end) {
  34. my $p = $primes->[$i];
  35. last if ($p > $y);
  36. next if ($p < $x);
  37. #~ my $valuation = valuation($p - 1, 2);
  38. #~ ($valuation > $k_exp and powmod($base, (($p - 1) >> $valuation) << $k_exp, $p) == ($congr % $p)) || next;
  39. my $t = $m * $p;
  40. if (($t - 1) % $lambda == 0 and ($t - 1) % ($p-1) == 0) {
  41. $callback->($t);
  42. }
  43. }
  44. return;
  45. }
  46. foreach my $i ($j .. $end) {
  47. my $p = $primes->[$i];
  48. last if ($p > $y);
  49. #~ my $valuation = valuation($p - 1, 2);
  50. #~ $valuation > $k_exp or next;
  51. #~ powmod($base, (($p - 1) >> $valuation) << $k_exp, $p) == ($congr % $p) or next;
  52. my $L = lcm($lambda, $p-1);
  53. gcd($L, $m) == 1 or next;
  54. my $t = $m * $p;
  55. my $u = divceil($A, $t);
  56. my $v = ($B / $t);
  57. if ($u <= $v) {
  58. __SUB__->($t, $L, $i + 1, $k - 1);
  59. }
  60. }
  61. }
  62. ->(Math::GMPz->new(1), 1, 0, $k);
  63. }
  64. use IO::Handle;
  65. open my $fh, '>>', 'carmichael_many_factors.txt';
  66. $fh->autoflush(1);
  67. my %upper_bounds = (
  68. 36 => Math::GMPz->new("172830055680118494946407003033666507461304818401153193809383963715892256751681"),
  69. 37 => Math::GMPz->new("804470457257926449746758080269993968890016257754008080494336091899208072210478721"),
  70. 38 => Math::GMPz->new("244899124403114685817402147257255073631462537923865013235929258099059306044154477281"),
  71. 39 => Math::GMPz->new("2912560918714425750692738781370955872381272347556033831319694306259522835520469570883201"),
  72. 40 => Math::GMPz->new("27919230451074589715843311695264905349211077611052606444369590994069578293094749438742401"),
  73. );
  74. #foreach my $lambda (80000 .. 1e6) {
  75. foreach my $lambda (sort {$a<=>$b} 15120, 30240, 110880, 285120, 332640, 498960, 604800, 1441440, 1663200, 1738800, 1814400, 2217600, 5216400, 13305600, 43243200, 64864800, 648648000, 4034016000, 8951342400, 12070749600, 67541947200) {
  76. #812700, 139230, 3197250, 4709250, 4709250, 2174130, 8824410, 20396250, 10442250, 982800, 7068600, 116953200, 88, 360, 3024, 12852, 8400, 39984, 18900, 486864, 529200) {
  77. #while (<>) {
  78. # chomp(my $lambda = $_);
  79. #$lambda >= 96600 or next;
  80. say "# Generating: $lambda";
  81. my @primes = grep { $_ > 2 and $lambda % $_ != 0 and is_prime($_) } map { $_ + 1 } divisors($lambda);
  82. foreach my $k (36..40) {
  83. last if ($k > scalar(@primes));
  84. #if (binomial(scalar(@primes), $k) < 1e6) {
  85. carmichael_in_range(Math::GMPz->new(~0), $upper_bounds{$k}, $k, \@primes, sub ($n) { say $n; say $fh $n; });
  86. #}
  87. }
  88. }