smooth_generate_quotient.pl 2.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129
  1. #!/usr/bin/perl
  2. # a(n) is the least k such that A348172(k) = n or 0 if no such k exists.
  3. # https://oeis.org/A348184
  4. # Known terms:
  5. # 4, 1, 9, 243, 3189375, 3176523, 10598252544
  6. # Conjectured upper-bounds:
  7. # a(8) <= 13011038208000
  8. # a(10) <= 1191916494900613125
  9. use 5.020;
  10. use warnings;
  11. use experimental qw(signatures);
  12. use Math::GMPz;
  13. use ntheory qw(:all);
  14. sub check_valuation ($n, $p) {
  15. if ($p == 2) {
  16. return valuation($n, $p) < 5;
  17. }
  18. if ($p == 3) {
  19. return valuation($n, $p) < 3;
  20. }
  21. if ($p == 7) {
  22. return valuation($n, $p) < 3;
  23. }
  24. ($n % $p) != 0;
  25. }
  26. sub smooth_numbers ($limit, $primes) {
  27. my @h = (1);
  28. foreach my $p (@$primes) {
  29. say "Prime: $p";
  30. foreach my $n (@h) {
  31. if ($n * $p < $limit) { #and check_valuation($n, $p)) {
  32. push @h, $n * $p;
  33. }
  34. }
  35. }
  36. return \@h;
  37. }
  38. #my $h = smooth_numbers(1e10, primes(7));
  39. my $h = smooth_numbers(6810951399432075+2, [2,3,5,7]);
  40. #my $h = smooth_numbers(13011038208000 * 31 + 1, primes(31));
  41. say "\nFound: ", scalar(@$h), " terms";
  42. my %best;
  43. foreach my $k (@$h) {
  44. #$k%2 == 0 or next;
  45. my %table;
  46. #for (my $i = 2; $i <= 10000; $i += 2) {
  47. for (my $i = 2; $i <= 1000; $i += 1) {
  48. my $n = mulint($i, $k);
  49. my $x = divisors($n);
  50. my $y = $n;
  51. my $g = gcd($x, $y);
  52. $x = divint($x, $g);
  53. $y = divint($y, $g);
  54. my $key = "$x/$y";
  55. $table{$key}{count}++;
  56. if (exists $table{$key}{max}) {
  57. if ($n < $table{$key}{max}) {
  58. $table{$key}{max} = $n;
  59. }
  60. }
  61. else {
  62. $table{$key}{max} = $n;
  63. }
  64. if ($table{$key}{count} >= 6) {
  65. if (exists $best{$table{$key}{count}}) {
  66. if ($table{$key}{max} < $best{$table{$key}{count}}) {
  67. $best{$table{$key}{count}} = $table{$key}{max};
  68. say "a($table{$key}{count}) = $best{$table{$key}{count}}";
  69. }
  70. }
  71. else {
  72. $best{$table{$key}{count}} = $table{$key}{max};
  73. say "a($table{$key}{count}) = $best{$table{$key}{count}}";
  74. }
  75. }
  76. }
  77. }
  78. __END__
  79. a(6) <= 7962624
  80. a(7) <= 362797056000
  81. a(6) <= 3176523
  82. a(8) <= 23544029528901
  83. a(10) <= 1191916494900613125
  84. a(6) <= 7962624
  85. a(7) <= 362797056000
  86. a(8) <= 13011038208000
  87. a(6) <= 384521484375
  88. a(7) <= 27751300048828125
  89. a(6) <= 3176523
  90. a(7) <= 689349609375
  91. a(8) <= 23544029528901
  92. a(10) = 1191916494900613125