search.pl 2.0 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889
  1. #!/usr/bin/perl
  2. # Smallest number m such that the GCD of the x's that satisfy sigma(x)=m is n.
  3. # https://oeis.org/A241625
  4. # a(127) = 4096
  5. # a(151) = 3465904
  6. # a(14) > 1017189995
  7. use utf8;
  8. use 5.020;
  9. use strict;
  10. use warnings;
  11. use ntheory qw(:all);
  12. use experimental qw(signatures);
  13. sub dynamicPreimage ($N, $L) {
  14. my %r = (1 => [1]);
  15. foreach my $l (@$L) {
  16. my %t;
  17. foreach my $pair (@$l) {
  18. my ($x, $y) = @$pair;
  19. foreach my $d (divisors(divint($N, $x))) {
  20. if (exists $r{$d}) {
  21. push @{$t{mulint($x, $d)}}, map { mulint($_, $y) } @{$r{$d}};
  22. }
  23. }
  24. }
  25. while (my ($k, $v) = each %t) {
  26. push @{$r{$k}}, @$v;
  27. }
  28. }
  29. return if not exists $r{$N};
  30. sort { $a <=> $b } @{$r{$N}};
  31. }
  32. sub cook_sigma ($N, $k) {
  33. my %L;
  34. foreach my $d (divisors($N)) {
  35. next if ($d == 1);
  36. foreach my $p (map { $_->[0] } factor_exp(subint($d, 1))) {
  37. my $q = addint(mulint($d, subint(powint($p, $k), 1)), 1);
  38. my $t = valuation($q, $p);
  39. next if ($t <= $k or ($t % $k) or $q != powint($p, $t));
  40. push @{$L{$p}}, [$d, powint($p, subint(divint($t, $k), 1))];
  41. }
  42. }
  43. [values %L];
  44. }
  45. sub inverse_sigma ($N, $k = 1) {
  46. ($N == 1) ? (1) : dynamicPreimage($N, cook_sigma($N, $k));
  47. }
  48. my %easy;
  49. @easy{
  50. 1, 2, 3, 4, 5, 7, 8, 9, 12, 13, 16, 25, 31, 80, 97, 18, 19, 22, 27, 29, 32, 36,
  51. 37, 43, 45, 49, 50, 61, 64, 67, 72, 73, 81, 91, 98, 100, 101, 106, 109, 121, 128, 129, 133, 134,
  52. 137, 146, 148, 149, 152, 157, 162, 163, 169, 171, 173, 192, 193, 197, 199, 200, 202, 211, 217, 218, 219
  53. } = ();
  54. my $count = 0;
  55. foreach my $k (1017189995 .. 1e10) {
  56. say "Testing: $k" if (++$count % 10000 == 0);
  57. my $t = gcd(inverse_sigma($k));
  58. if ($t >= 14 and $t <= 219) {
  59. say "\na($t) = $k\n" if not exists $easy{$t};
  60. die "Found: $k" if ($t == 14 or $t == 15);
  61. }
  62. }