123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302 |
- #!/usr/bin/perl
- # Erdos construction method for Carmichael numbers:
- # 1. Choose an even integer L with many prime factors.
- # 2. Let P be the set of primes d+1, where d|L and d+1 does not divide L.
- # 3. Find a subset S of P such that prod(S) == 1 (mod L). Then prod(S) is a Carmichael number.
- # Alternatively:
- # 3. Find a subset S of P such that prod(S) == prod(P) (mod L). Then prod(P) / prod(S) is a Carmichael number.
- use 5.020;
- use warnings;
- use ntheory qw(:all);
- use List::Util qw(shuffle);
- use experimental qw(signatures);
- use Math::GMPz;
- # Modular product of a list of integers
- sub vecprodmod ($arr, $mod) {
- #~ if ($mod > ~0) {
- #~ my $prod = Math::GMPz->new(1);
- #~ foreach my $k(@$arr) {
- #~ $prod = ($prod * $k) % $mod;
- #~ }
- #~ return $prod;
- #~ }
- if ($mod < ~0) {
- my $prod = 1;
- foreach my $k(@$arr) {
- $prod = mulmod($prod, $k, $mod);
- }
- return $prod;
- }
- my $prod = 1;
- foreach my $k (@$arr) {
- $prod = Math::Prime::Util::GMP::mulmod($prod, $k, $mod);
- }
- Math::GMPz->new($prod);
- }
- # Primes p such that p-1 divides L and p does not divide L
- sub lambda_primes ($L) {
- #grep { ($L % $_) != 0 } grep { $_ > 2 and is_prime($_) } map { $_ + 1 }
- #map { Math::GMPz->new($_)}
- # divisors($L);
- grep { ($_ > 2) and (($L % $_) != 0) and is_prime($_) } map { ($_ >= ~0) ? (Math::GMPz->new($_)+1) : ($_ + 1) } divisors($L);
- }
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89);
- #my @prefix = (3,5,17,23, 29);
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89);
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 617, 2003, 2549);
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 617);
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 419, 449, 617);
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003);
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003);
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 617, 1409, 2003, 2549, 3137, 9857, 10193, 16073, 68993, 202049, 275969, 1500929, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2549, 3137, 4019, 4289, 10193, 16073, 21977, 38459, 50513, 52529, 76649, 93809, 97553]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2549, 8009, 9857, 23297, 50513, 68993, 275969, 375233, 1500929, 3232769, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2297, 2549, 3329, 8009, 10193, 16073, 23297, 50177, 93809, 202049, 275969, 656657, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297, 68993, 88397, 93809, 202049, 896897, 1500929, 2475089, 8044037, 18386369]
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019);
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353);
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019);
- #my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019);
- # primes(2, 15500).grep{ .dec.is_smooth(43) }.grep{ gcd(.phi, [3, 5, 17, 23, 29].prod) == 1 }
- # => [2, 3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 173, 197, 257, 353, 419, 449, 593, 617, 677, 683, 947, 1217, 1373, 1409, 1559, 1613, 2003, 2129, 2237, 2297, 2357, 2543, 2549, 2663, 2729, 2753, 2927, 3137, 3257, 3329, 3389, 3527, 3719, 4019, 4733, 4817, 5333, 5477, 6449, 6689, 6917, 7253, 7547, 7937, 8009, 8429, 8513, 8867, 9473, 9857, 9923, 10169, 10193, 12959, 13313, 13469, 13553, 13859, 14897, 15137, 15377, 15467]
- my @prefix = (
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297, 68993, 88397, 93809, 202049, 896897
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 9923, 10193, 13553, 308309, 486179, 656657
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 9923, 10193, 13553, 308309, 486179, 656657, 896897, 902903, 1500929, 1610753, 1944713, 2063777, 2540033, 2818817, 3232769, 6071297, 18386369
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 9923, 10193, 13553, 308309, 486179, 656657, 896897, 902903, 1500929, 1610753, 1944713, 2063777, 2540033, 2818817, 3232769, 6071297, 18386369, 22629377, 25281257, 31115393, 33020417
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369, 2059273217
- # 3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297, 68993, 88397, 93809, 202049, 896897, 1500929, 2475089, 8044037, 18386369, 2059273217
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 10193, 3232769, 6071297, 18386369
- # 3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 449, 617, 2003, 2297, 2549, 3137, 10193, 18386369
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297,
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 617, 1409, 2003, 2549, 3137, 9857, 10193, 16073, 68993, 18386369
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 4019, 4289, 7547, 9857, 10193, 16073, 68993, 88397, 93809
- #3, 5, 17, 23, 29, 53, 89, 197
- #3, 5, 17, 23, 29, 53, 89, 197, 617
- # 3, 5, 17, 23, 29, 197, 617, 1217, 6689, 14897, 16073, 20483, 46817, 50513, 68993, 76343, 81929,
- #3, 5, 17, 23, 29, 197, 419, 449, 617, 1217, 2129, 3137, 9857, 10193, 16073, 68993, 18386369
- #3, 5, 17, 23, 29, 43, 53, 89, 113, 127, 157, 257, 353, 397, 449, 617, 1093, 1499, 2731
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 353,
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 173, 197, 257, 263, 269, 293, 353, 359, 383
- # 3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 173, 197, 257, 353
- # 3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 173, 197, 257, 353, 419, 449, 593, 617, 677, 683, 947, 1217, 1373,
- # 3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 197, 257, 353, 419, 449, 593, 617, 677, 683, 1217, 1373, 1409, 1559, 1613
- # 3, 5, 17, 23, 29, 53, 83, 89, 113
- #~ 3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 173, 197, 257, 263, 269, 293, 353, 359, 383
- #~ A328691
- #~ A329460
- #~ #Constructing Carmichael numbers throughimproved subset-product algorithms
- #~ 1.893
- #~ The smallest Abundant number which is also pseudoprime (base-2).
- #~ Shyam Sunder Gupta
- #~ 171800042106877185
- #~ 10-11-2019
- #~ 3470207934739664512679701940114447720865
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019
- #3, 5, 17, 23, 29, 53, 83, 89, 2237, 2297, 3527, 8009, 15137, 24683, 26489,
- #3, 5, 17, 23, 29, 53, 83, 89, 2237, 2297, 3527, 8009, 15137, 24683, 26489, 49193, 50513, 56417, 77573, 93809, 98729, 105953, 202049, 250433
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 173, 197, 257, 353, 419, 449, 593, 617, 677, 683, 947, 1217, 1373, 1409, 1559
- #3, 5, 17, 23, 29, 53, 89, 113, 197, 257, 353, 419, 449, 617, 677, 1217, 1373, 1409, 2003, 2129, 2549, 2663, 2927, 3137, 3329, 3389, 3719, 4733, 6689, 6917, 7547, 8009, 8513, 9857, 10193, 13313, 13553, 14897
- #3, 5, 17, 23, 29, 53, 89, 113, 197, 257, 353, 419, 449, 617, 677, 1217, 1373, 1409, 2003, 2129, 2549, 2663, 2927, 3137, 3329, 3389, 3719, 4733, 6689, 6917, 7547, 8009, 8513, 9857, 10193, 13313, 13553, 14897, 15809, 17837, 18773, 19457, 20483, 21737, 23297, 25169, 27437, 30977, 31769, 34607, 35153, 38039, 40433, 46817,
- #3, 5, 17, 23, 29, 89, 113, 197, 257, 353, 449, 617, 1373, 1409, 2663, 3137, 3389, 7547, 9857, 13553, 30977, 50177, 65537, 68993, 114689, 120737, 130439, 166013, 170369, 268913, 275969, 470597, 495617, 521753, 614657, 913067, 1075649, 1103873
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 197, 257, 353, 419, 449, 593, 617, 677, 683, 1217, 1373, 1409, 1559, 1613, 2003, 2129, 2297, 2357, 2543, 2549, 2663, 2729, 2927, 3137, 3257, 3329, 3389, 3719, 4019, 4733, 5477, 6449, 6689
- #256392945019638502140697184374123648403371443140710905
- 3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 9923
- #3, 5, 17, 23, 29, 53, 83, 89, 113, 149, 173, 197, 257, 263, 269, 293, 353, 359, 383
- );
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 617, 1409, 2003, 2549, 3137, 9857, 10193, 16073, 68993, 202049, 275969, 1500929, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2549, 3137, 4019, 4289, 10193, 16073, 21977, 38459, 50513, 52529, 76649, 93809, 97553]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2549, 8009, 9857, 23297, 50513, 68993, 275969, 375233, 1500929, 3232769, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369, 2059273217]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2297, 2549, 3329, 8009, 10193, 16073, 23297, 50177, 93809, 202049, 275969, 656657, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2549, 8009, 9857, 23297, 50513, 68993, 275969, 375233, 1500929, 3232769, 18386369, 1176727553]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297, 68993, 88397, 93809, 202049, 896897, 1500929, 2475089, 8044037, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 9923, 10193, 13553, 308309, 486179, 656657, 896897, 902903, 1500929, 1610753, 1944713, 2063777, 2540033, 2818817, 3232769, 6071297, 18386369, 22629377, 25281257, 31115393, 33020417]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2549, 3137, 4019, 4289, 10193, 16073, 21977, 38459, 50513, 52529, 76649, 93809, 97553]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2549, 8009, 9857, 23297, 50513, 68993, 275969, 375233, 1500929, 3232769, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2297, 2549, 3329, 8009, 10193, 16073, 23297, 50177, 93809, 202049, 275969, 656657, 18386369]
- #~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297, 68993, 88397, 93809, 202049, 896897, 1500929, 2475089, 8044037, 18386369]
- #my @prefix = (3, 5, 17, 23, 29);
- my $prefix_prod = Math::GMPz->new(vecprod(@prefix));
- my $multiplier = lcm(map {$_-1} @prefix);
- sub method_1 ($L) { # smallest numbers first
- (vecall { ($L % ($_-1)) == 0 } @prefix) or return;
- my @P = lambda_primes($L);
- @P = grep {
- (($prefix_prod % $_) == 0)
- ? 1
- : Math::Prime::Util::GMP::gcd(Math::Prime::Util::GMP::totient(Math::Prime::Util::GMP::mulint($prefix_prod, $_)), Math::Prime::Util::GMP::mulint($prefix_prod, $_)) eq '1'
- } @P;
- #vecprodmod(@P, 3*5*17*23) == 0 or return;
- #vecprodmod(\@P, 3*5*17*23*29) == 0 or return;
- if (@prefix) {
- vecprodmod(\@P, $prefix_prod) == 0 or return;
- }
- #return if (vecprod(@P) < ~0);
- # @P = grep { $_ > 257 } @P;
- #@P = grep { $_ > $prefix[-1] } @P;
- @P = grep { gcd($prefix_prod, $_) == 1 } @P;
- say "# Testing: $L -- ", scalar(@P);
- my $n = scalar(@P);
- my @orig = @P;
- my $max = 1e3;
- my $max_k = 30;
- my $L_rem = invmod($prefix_prod, $L);
- foreach my $k (1 .. @P) {
- #next if (binomial($n, $k) > 1e6);
- next if ($k > $max_k);
- @P = @orig;
- my $count = 0;
- forcomb {
- if (vecprodmod([@P[@_]], $L) == $L_rem) {
- say vecprod(@P[@_], $prefix_prod);
- }
- lastfor if (++$count > $max);
- } $n, $k;
- next if (binomial($n, $k) < $max);
- @P = reverse(@P);
- $count = 0;
- forcomb {
- if (vecprodmod([@P[@_]], $L) == $L_rem) {
- say vecprod(@P[@_], $prefix_prod);
- }
- lastfor if (++$count > $max);
- } $n, $k;
- for (1..2) {
- @P = shuffle(@P);
- $count = 0;
- forcomb {
- if (vecprodmod([@P[@_]], $L) == $L_rem) {
- say vecprod(@P[@_], $prefix_prod);
- }
- lastfor if (++$count > $max);
- } $n, $k;
- }
- }
- my $B = vecprodmod([@P, $prefix_prod], $L);
- my $T = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P));
- foreach my $k (1 .. @P) {
- #next if (binomial($n, $k) > 1e6);
- last if ($k > $max_k);
- @P = @orig;
- my $count = 0;
- forcomb {
- if (vecprodmod([@P[@_]], $L) == $B) {
- my $S = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P[@_]));
- say vecprod($prefix_prod, $T / $S) if ($T != $S);
- }
- lastfor if (++$count > $max);
- } $n, $k;
- next if (binomial($n, $k) < $max);
- @P = reverse(@P);
- $count = 0;
- forcomb {
- if (vecprodmod([@P[@_]], $L) == $B) {
- my $S = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P[@_]));
- say vecprod($prefix_prod, $T / $S) if ($T != $S);
- }
- lastfor if (++$count > $max);
- } $n, $k;
- for (1..2) {
- @P = shuffle(@P);
- $count = 0;
- forcomb {
- if (vecprodmod([@P[@_]], $L) == $B) {
- my $S = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P[@_]));
- say vecprod($prefix_prod, $T / $S) if ($T != $S);
- }
- lastfor if (++$count > $max);
- } $n, $k;
- }
- }
- }
- use Math::GMPz;
- my %seen;
- my $count = 0;
- foreach my $n(1..1e4) {
- #method_1($n*3236000768);
- #method_1($n* 3236000768);
- #method_1($n * 18386368);
- # ($n == 1) or gcd($n, $multiplier) > 1 or next;
- method_1(mulint($n, $multiplier));
- }
|