smooth_generate.pl 1.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596
  1. #!/usr/bin/perl
  2. # Least k such that phi(k) is an n-th power when k is the product of n distinct primes.
  3. # https://oeis.org/A281069
  4. # Upper-bounds:
  5. # a(11) <= 5703406198237930
  6. # a(12) <= 178435136773443810
  7. # a(13) <= 4961806417984478790
  8. # a(14) <= 1686255045377006404458210
  9. # a(15) <= 93685760130566580980309670
  10. # Other upper-bounds:
  11. # a(11) < 7829180118279030 < 7862926281269430
  12. # a(12) < 202848252169827810 < 771632451121296274230
  13. # a(13) < 4976768607520456110
  14. # a(14) < 1920340956699420394773270
  15. # a(15) < 35141760969505011678479285910
  16. use 5.020;
  17. use warnings;
  18. use experimental qw(signatures);
  19. use Math::GMPz;
  20. use ntheory qw(:all);
  21. sub check_valuation ($n, $p) {
  22. ($n % $p) != 0;
  23. }
  24. sub smooth_numbers ($limit, $primes) {
  25. my @h = (Math::GMPz->new(1));
  26. foreach my $p (@$primes) {
  27. say "Prime: $p";
  28. foreach my $n (@h) {
  29. if ($n * $p <= $limit and check_valuation($n, $p)) {
  30. push @h, $n * $p;
  31. }
  32. }
  33. }
  34. return \@h;
  35. }
  36. sub isok ($n) {
  37. is_power(euler_phi($n));
  38. }
  39. my @smooth_primes;
  40. foreach my $p (@{primes(1e8)}) {
  41. if ($p == 2) {
  42. push @smooth_primes, $p;
  43. next;
  44. }
  45. if (is_smooth($p-1, 3)) {
  46. push @smooth_primes, $p;
  47. }
  48. }
  49. my $multiple = 1254855;
  50. #my $multiple = 26535163830;
  51. #my $multiple = 7474975358110;
  52. #my $multiple = 1637019603426090;
  53. @smooth_primes = grep { gcd($multiple, $_) == 1 } @smooth_primes;
  54. my $h = smooth_numbers(~0, \@smooth_primes);
  55. say "\nFound: ", scalar(@$h), " terms";
  56. my %table;
  57. foreach my $k (@$h) {
  58. my $n = mulint($k, $multiple);
  59. my $p = isok($n);
  60. if ($p >= 11 and prime_omega($n) == $p) {
  61. say "a($p) = $n -> ", join(' * ', map { "$_->[0]^$_->[1]" } factor_exp($n));
  62. push @{$table{$p}}, $n;
  63. }
  64. }
  65. say '';
  66. foreach my $k (sort { $a <=> $b } keys %table) {
  67. say "a($k) <= ", vecmin(@{$table{$k}});
  68. }