prime_gaps_db.pl 2.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778
  1. #!/usr/bin/perl
  2. # Find prime gap records, from the factors of pseudoprimes.
  3. use 5.036;
  4. use Storable;
  5. use POSIX qw(ULONG_MAX);
  6. use Math::GMPz;
  7. use List::Util qw(uniq);
  8. use ntheory qw(:all);
  9. use Math::Prime::Util::GMP;
  10. use experimental qw(signatures);
  11. use constant {
  12. USE_MPFR => 0, # true to check primes greater than 300
  13. };
  14. eval { require GDBM_File };
  15. my $cache_db = "cache/factors.db";
  16. dbmopen(my %db, $cache_db, 0444)
  17. or die "Can't create/access database <<$cache_db>>: $!";
  18. my $z = Math::GMPz::Rmpz_init();
  19. my @results;
  20. my $max_merit = 0;
  21. use Math::MPFR;
  22. my $f = Math::MPFR::Rmpfr_init2(128);
  23. while (my ($key, $value) = each %db) {
  24. my @factors = split(' ', $value);
  25. $factors[-1] > 18361375334787046697 or next;
  26. foreach my $p (@factors) {
  27. if ($p > 18361375334787046697 and (USE_MPFR ? 1 : (length($p) <= 300))) {
  28. #$p = Math::Prime::Util::GMP::next_prime($p);
  29. #$p = Math::Prime::Util::GMP::prev_prime($p);
  30. my $gap = Math::Prime::Util::GMP::subint(Math::Prime::Util::GMP::next_prime($p), $p);
  31. my $merit = USE_MPFR
  32. ? do {
  33. Math::MPFR::Rmpfr_set_str($f, $p, 10, 0);
  34. Math::MPFR::Rmpfr_log($f, $f, 0);
  35. Math::MPFR::Rmpfr_ui_div($f, $gap, $f, 0);
  36. Math::MPFR::Rmpfr_get_d($f, 0);
  37. }
  38. : do { $gap / log($p) };
  39. if ($merit > $max_merit) {
  40. say "Merit: $merit with g = $gap for p = $p";
  41. $max_merit = $merit;
  42. }
  43. }
  44. }
  45. }
  46. __END__
  47. Merit: 1.25577454904363 with g = 88 for p = 2714804953442921912706470542501
  48. Merit: 3.05208217232875 with g = 148 for p = 1147056470270532373861
  49. Merit: 3.53787442283434 with g = 210 for p = 60077339080249258834663873
  50. Merit: 3.95233801271049 with g = 274 for p = 1282102862013350567377590730597
  51. Merit: 4.11829942456146 with g = 226 for p = 680435527254568598106721
  52. Merit: 4.38561379855734 with g = 220 for p = 6108790953688828846561
  53. Merit: 5.82605203757676 with g = 300 for p = 23070597828321430685567
  54. Merit: 6.14519150536919 with g = 678 for p = 823725213912413947530390077834174890593144661261
  55. Merit: 7.36779003891085 with g = 406 for p = 854427865432681725179473
  56. Merit: 8.57040810150965 with g = 462 for p = 257786279140373792788201
  57. Merit: 10.397727286833 with g = 700 for p = 172881565818578615860856969761
  58. Merit: 11.5899202566782 with g = 586 for p = 9087403522138538233297
  59. Merit: 14.5130453759809 with g = 798 for p = 758032929981873202950721