123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272 |
- #!/usr/bin/perl
- use 5.020;
- use warnings;
- use strict;
- use ntheory qw(:all);
- use experimental qw(signatures);
- use List::Util qw(uniqnum);
- use File::Find qw(find);
- use Math::GMPz;
- use Math::AnyNum;
- use Math::Prime::Util::GMP qw(is_euler_pseudoprime);
- sub is_absolute_euler_pseudoprime ($n) {
- is_carmichael($n) and vecall { (($n-1)>>1) % ($_-1) == 0 } factor($n);
- }
- # Let b(n) be the smallest odd composite k such that q^((k-1)/2) == -1 (mod k) for every prime q <= prime(n).
- # b(1) = 3277
- # b(2) = 1530787
- # b(3) = 3697278427
- # b(4) = 118670087467
- # b(5) <= 2152302898747
- # b(6) <= 614796634515444067
- # b(7) <= 614796634515444067
- # 341, 1729, 1729, 46657
- my $N = 13;
- my $P = nth_prime($N);
- my $MAX = ~0;
- my @primes_bellow = @{primes($P)};
- #~ my $n = 1;
- #~ my $P = 2;
- my @primes = (2);
- sub non_residue {
- my ($n) = @_;
- foreach my $p (@primes) {
- #if ($p > $P) {
- # return -1;
- #}
- #sqrtmod($p, $n) // return $p;
- #sqrtmod($p, $n) // do { say "$p $n"; return $p};
- #kronecker($n, $p) == 1 and return $p;
- (vecall { kronecker($p, $_->[0]) == 1 } factor_exp($n)) || return $p;
- }
- return -1;
- }
- #~ foroddcomposites {
- #~ my $k = $_;
- #~ #if (powmod($q^((k-1)/2) == -1 and non_residue($k) == $P) {
- #~ if (powmod($P, ($k-1)>>1, $k) == $k-1) {
- #~ my $q = non_residue($k);
- #~ if ($q == $P) {
- #~ say $k;
- #~ $P = next_prime($P);
- #~ push @primes, $P;
- #~ exit if $P == 13;
- #~ }
- #~ }
- #~ } 1e9;
- #~ __END__
- sub isok {
- my ($k) = @_;
- if (powmod($P, ($k-1)>>1, $k) == $k-1) {
- my $q = non_residue($k);
- return ($P == $q);
- }
- return;
- #~ my $res;
- #~ my $p = nth_prime($n);
- #~ for(my $k = $A000229[$n-1]*$A000229[$n-1]; ; $k+=2) {
- #~ if (!is_prime($k) and powmod($p, ($k-1)>>1, $k) == $k-1) {
- #~ my $q = non_residue($k, $p);
- #~ if ($p == $q) {
- #~ return $k;
- #~ }
- #~ }
- #~ }
- #~ return $res;
- }
- #~ Found: 2095988771234401
- #~ Found: 35488761020833
- #~ Found: 4364724604902481
- #~ Found: 5087540155980241
- #~ Found: 619440406020833
- sub isok_12_may {
- my ($k) = @_;
- # Odd composite numbers n such that b^((n-1)/2) == +-1 (mod n) for every natural b < log(n).
- my $w = $k-1;
- my $z = $w>>1;
- my $upto = int(log("$k"));
- #~ if (lc($upto) eq 'inf') {
- #~ $upto = Math::AnyNum->new($k)->ilog;
- #~ }
- #for (my $p = 2; $p <= $upto; ++$p) {
- for(my $p = 2; $p <= $upto; $p = next_prime($p)) {
- my $t = powmod($p, $z, $k);
- ($t == 1) || ($t == $w) || return;
- }
- return 1;
- }
- #~ say isok_12_may(35488761020833);
- #~ say isok_12_may(2095988771234401);
- #~ say is_absolute_euler_pseudoprime(35488761020833);
- #~ __END__
- sub isok_12_may_stronger {
- my ($k) = @_;
- # Odd composite numbers n such that b^((n-1)/2) == +-1 (mod n) for every natural b < log(n).
- my $w = $k-1;
- my $z = $w>>1;
- my $upto = int(log("$k"));
- #~ if (lc($upto) eq 'inf') {
- #~ $upto = Math::AnyNum->new($k)->ilog;
- #~ }
- for(my $p = 2; $p <= $upto; $p = next_prime($p)) {
- my $t = powmod($p, $z, $k);
- ($t == 1) || return;
- }
- return 1;
- }
- #~ #foreach my $k(1..1e6) {
- #~ $| = 1;
- #~ #for (my $k = 3; $k <= 1e9; $k += 2) {
- #~ foroddcomposites {
- #~ my $k = $_;
- #~ #if (not is_prime($k) and isok($k)) {
- #~ if (isok($k)) {
- #~ print($k, ", ");
- #~ if (not isok_stronger($k)) {
- #~ die "Counter-example: $k";
- #~ }
- #~ }
- #~ } 1e9;
- #~ __END__
- sub isok_b {
- my ($k) = @_;
- #$k % 2 == 1 or return;
- vecall { powmod($_, ($k-1)>>1, $k) == $k-1 } @primes_bellow;
- }
- sub isok_b2 {
- my ($k) = @_;
- $k % 2 == 1 or return;
- vecall { powmod($_, ($k-1)>>1, $k) == 1 } 2..($P-1);
- }
- my %seen;
- sub process_file {
- my ($file) = @_;
- open my $fh, '<', $file;
- while (<$fh>) {
- next if /^\h*#/;
- /\S/ or next;
- my $n = (split(' ', $_))[-1];
- $n || next;
- #if ($n > $MAX or $n <= 2) {
- # next;
- #}
- #~ if (length($n) > 30) {
- #~ next;
- #~ }
- if ($n > 1.7 * 10**16) {
- next;
- }
- #if ($n > ~0) {
- #if (length($n) > 30) {
- # next;
- #}
- next if $seen{$n}++;
- if ($n > ~0) {
- $n = Math::GMPz->new("$n");
- }
- #next if is_prime($n);
- #if (isok_b($n)) {
- # if (not is_absolute_euler_pseudoprime($n) and isok_12_may($n) ) {
- if (isok_12_may($n) and not isok_12_may_stronger($n)) {
- say "Counter-example: $n";
- #$MAX = $n;
- # if (not isok_stronger($n)) {
- # die "Counter-example: $n";
- # }
- #~ if ($n < $MAX) {
- #~ $MAX = $n;
- #~ say "New record: $n";
- #~ }
- }
- }
- close $fh;
- }
- my $psp = "/home/swampyx/Other/Programare/experimental-projects/pseudoprimes/oeis-pseudoprimes";
- find({
- wanted => sub {
- if ( /\.txt\z/) {
- #say "Processing $_";
- process_file($_);
- }
- },
- no_chdir => 1,
- } => $psp);
|