test-pseudoprime_search_2.pl 4.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215
  1. #!/usr/bin/perl
  2. use 5.020;
  3. use warnings;
  4. use strict;
  5. use ntheory qw(:all);
  6. use experimental qw(signatures);
  7. use List::Util qw(uniqnum);
  8. use File::Find qw(find);
  9. use Math::GMPz;
  10. use Math::AnyNum qw(prod);
  11. use Math::Prime::Util::GMP qw(is_euler_pseudoprime);
  12. sub is_absolute_euler_pseudoprime ($n) {
  13. is_carmichael($n) and vecall { (($n-1)>>1) % ($_-1) == 0 } factor($n);
  14. }
  15. my $N = 21;
  16. my $P = nth_prime($N);
  17. my $MAX = 10**300;
  18. my @primes_bellow = @{primes($P)};
  19. #~ my $n = 1;
  20. #~ my $P = 2;
  21. my @primes = (2);
  22. sub non_residue {
  23. my ($n) = @_;
  24. foreach my $p (@primes) {
  25. sqrtmod($p, $n) // return $p;
  26. }
  27. return -1;
  28. }
  29. #~ foroddcomposites {
  30. #~ my $k = $_;
  31. #~ #if (powmod($q^((k-1)/2) == -1 and non_residue($k) == $P) {
  32. #~ if (powmod($P, ($k-1)>>1, $k) == $k-1) {
  33. #~ my $q = non_residue($k);
  34. #~ if ($q == $P) {
  35. #~ say $k;
  36. #~ $P = next_prime($P);
  37. #~ push @primes, $P;
  38. #~ exit if $P == 13;
  39. #~ }
  40. #~ }
  41. #~ } 1e9;
  42. #~ __END__
  43. sub isok {
  44. my ($k) = @_;
  45. my $w = ($k-1)>>1;
  46. vecall { powmod($_, $w, $k) == (kronecker($_, $k) % $k) } 1..$P;
  47. }
  48. # Carmichael numbers of the form (6*k+1)*(12*k+1)*(18*k+1), where 6*k+1, 12*k+1 and 18*k+1 are all primes.
  49. # Carmichael numbers of the form C = (30n-p)*(60n-(2p+1))*(90n-(3p+2)), where n is a natural number and p, 2p+1, 3p+2 are all three prime numbers.
  50. # Numbers of the form: (6*m + 1) * (12*m + 1) * Product_{i=1..k-2} (9 * 2^i * m + 1), where k >= 3, with the condition that each of the factors is prime and that m is divisible by 2^(k-4).
  51. #~ foreach my $k(1..1e7) {
  52. #~ my $x = 6*$k+1;
  53. #~ my $y = 12*$k + 1;
  54. #~ my $z = 18*$k+1;
  55. #~ #my $w = 9 * 2**2 * $k+1;
  56. #~ if (is_prime($x) and is_prime($y) and is_prime($z)
  57. #~ #and is_prime($w)
  58. #~ ) {
  59. #~ my $n = prod($x, $y, $z);
  60. #~ if (isok($n)) {
  61. #~ say "[$k] -> a($N) <= $n";
  62. #~ }
  63. #~ }
  64. #~ }
  65. #~ __END__
  66. #~ #foreach my $k(1..1e6) {
  67. #~ $| = 1;
  68. #~ #for (my $k = 3; $k <= 1e9; $k += 2) {
  69. #~ foroddcomposites {
  70. #~ my $k = $_;
  71. #~ #if (not is_prime($k) and isok($k)) {
  72. #~ if (isok($k)) {
  73. #~ print($k, ", ");
  74. #~ if (not isok_stronger($k)) {
  75. #~ die "Counter-example: $k";
  76. #~ }
  77. #~ }
  78. #~ } 1e9;
  79. #~ __END__
  80. my %seen;
  81. sub process_file {
  82. my ($file) = @_;
  83. open my $fh, '<', $file;
  84. while (<$fh>) {
  85. next if /^\h*#/;
  86. /\S/ or next;
  87. my $n = (split(' ', $_))[-1];
  88. $n || next;
  89. #if ($n > $MAX or $n <= 2) {
  90. # next;
  91. #}
  92. #~ if (length($n) > 30) {
  93. #~ next;
  94. #~ }
  95. #~ if ($n < 14469841 or $n > $MAX) {
  96. #~ next;
  97. #~ }
  98. #~ if ($n < ~0) {
  99. #~ next;
  100. #~ }
  101. #if ($n < 619440406020833) {
  102. # next;
  103. #}
  104. #if ($n < 1.7 * 10**16) {
  105. #~ next;
  106. #~ }
  107. #~ if ($n > ~0) {
  108. #~ next;
  109. #~ }
  110. #~ if ($n > 10**8) {
  111. #~ next;
  112. #~ }
  113. #if ($n > ~0) {
  114. #if (length($n) > 30) {
  115. # next;
  116. #}
  117. next if $seen{$n}++;
  118. if ($n > ~0) {
  119. $n = Math::GMPz->new("$n");
  120. }
  121. #next if is_prime($n);
  122. #if (isok_b($n)) {
  123. # if (not is_absolute_euler_pseudoprime($n) and isok_12_may($n) ) {
  124. # say "Testing: $n";
  125. #if (isok_12_may($n) and not is_carmichael($n)) {#not isok_12_may_stronger($n)) {
  126. # if (isok_12_may($n) and not isok_12_may_stronger($n)) {# not is_carmichael($n)) {
  127. if (isok($n)) {
  128. #say "Counter-example: $n";
  129. say "a($N) = $n";
  130. if ($n < $MAX) {
  131. $MAX = $n;
  132. }
  133. #last if ($n > 15851273401);
  134. #$MAX = $n;
  135. # if (not isok_stronger($n)) {
  136. # die "Counter-example: $n";
  137. # }
  138. #~ if ($n < $MAX) {
  139. #~ $MAX = $n;
  140. #~ say "New record: $n";
  141. #~ }
  142. }
  143. }
  144. close $fh;
  145. }
  146. my $psp = "/home/swampyx/Other/Programare/experimental-projects/pseudoprimes/oeis-pseudoprimes";
  147. find({
  148. wanted => sub {
  149. if (/\.txt\z/) {
  150. #say "Processing $_";
  151. process_file($_);
  152. }
  153. },
  154. no_chdir => 1,
  155. } => $psp);