123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216 |
- #!/usr/bin/perl
- # Smallest number m such that the GCD of the x's that satisfy sigma(x)=m is n.
- # https://oeis.org/A241625
- # Known terms:
- # 1, 3, 4, 7, 6, 6187272, 8, 15, 13, 196602, 8105688, 28, 14
- # Upper-bounds for a(144):
- # a(144) <= {5334186182, 6653667020, 15812165226, 26999284026, 28752564542, 31717503246, 51526204158, 68240139110, 71641633206, 73884376566, 320692161764, 404350986572, 724592550812, 921030858722, 1138698021812, 7340065485194, 10562405788532, 405101345685452, 1058995989058904, 1293672253540082, 1991013491138282, 2643933035123354, 15645437418062582}
- # Other upper-bounds:
- # a(127) <= {4096, 1000153088, 1006583808}
- # Lower-bound for a(13):
- # a(13) > 1005756000
- # Upper-bounds for larger terms:
- # a(1444) <= 2667
- # a(4096) <= 8191
- # a(36100) <= 82677
- # a(94633984) <= 199753347
- use 5.014;
- use Math::GMPz;
- use Math::AnyNum qw(is_smooth);
- use ntheory qw(:all);
- use Memoize qw(memoize);
- use List::Util qw(uniq);
- use experimental qw(signatures);
- memoize('max_power');
- my @smooth_primes;
- sub is_smooth_for_e {
- my ($p, $e) = @_;
- is_smooth(Math::GMPz->new($p)**($e - 1), 7)
- and is_smooth(Math::GMPz->new($p)**($e + 1) - 1, 7);
- }
- sub p_is_smooth {
- my ($p) = @_;
- vecany {
- is_smooth_for_e($p, $_);
- }
- 1 .. 20;
- }
- sub max_power {
- my ($p) = @_;
- for (my $e = 20 ; $e >= 1 ; --$e) {
- if (is_smooth_for_e($p, $e)) {
- return $e;
- }
- }
- }
- forprimes {
- if ($_ == 2) {
- push @smooth_primes, $_;
- }
- else {
- if (p_is_smooth($_)) {
- push @smooth_primes, $_;
- }
- }
- } 4801;
- say "@smooth_primes";
- push @smooth_primes, (2, 13, 23, 29, 31, 41, 83, 127, 269, 1153);
- push @smooth_primes, (2, 3, 7, 11, 13, 31, 71, 151, 1009, 2833);
- @smooth_primes = uniq(@smooth_primes);
- @smooth_primes = sort { $a <=> $b } @smooth_primes;
- foreach my $p (@smooth_primes) {
- say "a($p) = ", max_power($p);
- }
- sub check_valuation ($n, $p) {
- $n % ($p * $p) != 0;
- }
- sub hamming_numbers ($limit, $primes) {
- my @h = (1);
- foreach my $p (@$primes) {
- say "Prime: $p";
- foreach my $n (@h) {
- if ($n * $p <= $limit and check_valuation($n, $p)) {
- push @h, $n * $p;
- }
- }
- }
- return \@h;
- }
- sub isok {
- my ($n) = @_;
- my $t = Math::GMPz->new(divisor_sum($n)) * Math::GMPz->new(euler_phi($n));
- is_power($t);
- }
- #use 5.016;
- #use Math::AnyNum qw(:overload);
- #die join ' ', inverse_sigma(19594645850);
- my $h = hamming_numbers(1e12, \@smooth_primes);
- say "Found: ", scalar(@$h), " terms";
- my %table;
- use utf8;
- use 5.020;
- use strict;
- use warnings;
- use integer;
- use ntheory qw(:all);
- use experimental qw(signatures);
- use List::Util qw(uniq);
- binmode(STDOUT, ':utf8');
- sub dynamicPreimage ($N, $L) {
- my %r = (1 => [1]);
- foreach my $l (@$L) {
- my %t;
- foreach my $pair (@$l) {
- my ($x, $y) = @$pair;
- foreach my $d (divisors(divint($N, $x))) {
- if (exists $r{$d}) {
- push @{$t{mulint($x, $d)}}, map { mulint($_, $y) } @{$r{$d}};
- }
- }
- }
- while (my ($k, $v) = each %t) {
- push @{$r{$k}}, @$v;
- }
- }
- return if not exists $r{$N};
- sort { $a <=> $b } @{$r{$N}};
- }
- sub cook_sigma ($N, $k) {
- my %L;
- foreach my $d (divisors($N)) {
- next if ($d == 1);
- foreach my $p (map { $_->[0] } factor_exp(subint($d, 1))) {
- my $q = addint(mulint($d, subint(powint($p, $k), 1)), 1);
- my $t = valuation($q, $p);
- next if ($t <= $k or ($t % $k) or $q != powint($p, $t));
- push @{$L{$p}}, [$d, powint($p, subint(divint($t, $k), 1))];
- }
- }
- [values %L];
- }
- sub inverse_sigma ($N, $k = 1) {
- ($N == 1) ? (1) : dynamicPreimage($N, cook_sigma($N, $k));
- }
- my %easy;
- @easy{
- 1, 2, 3, 4, 5, 7, 8, 9, 12, 13, 16, 25, 31, 80, 97, 18, 19, 22, 27, 29, 32, 36,
- 37, 43, 45, 49, 50, 61, 64, 67, 72, 73, 81, 91, 98, 100, 101, 106, 109, 121, 128, 129, 133, 134,
- 137, 146, 148, 149, 152, 157, 162, 163, 169, 171, 173, 192, 193, 197, 199, 200, 202, 211, 217, 218, 219
- } = ();
- my $count = 0;
- @$h = sort { $a <=> $b } @$h;
- #~ foreach my $n(@$h) {
- #~ if ($n == 5334186182) {
- #~ say "OK";
- #~ last;
- #~ }
- #~ }
- #~ say gcd(inverse_sigma(5334186182));
- #~ exit;
- foreach my $k (@$h) {
- say "Testing: $k" if (++$count % 1000 == 0);
- $k > 1e9 or next;
- my $t = gcd(inverse_sigma($k));
- if ($t >= 14 and $t <= 200) {
- say "\na($t) = $k\n" if not exists $easy{$t};
- if ($t == 14 or $t == 15) {
- say "Found: a($t) = $k";
- }
- }
- }
|