123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112 |
- #!/usr/bin/perl
- # Given a positive integer `n`, this algorithm finds all the numbers k
- # such that sigma(k) = n, where `sigma(k)` is the sum of divisors of `k`.
- # Based on "invphi.gp" code by Max Alekseyev.
- # See also:
- # https://home.gwu.edu/~maxal/gpscripts/
- use utf8;
- use 5.020;
- use strict;
- use warnings;
- use Math::Prime::Util::GMP qw(:all);
- use List::Util qw(uniq);
- use experimental qw(signatures);
- binmode(STDOUT, ':utf8');
- sub inverse_sigma {
- my ($n) = @_;
- my %cache;
- my %factor_cache;
- my %divisor_cache;
- my $results = sub ($n, $m) {
- return [1] if ($n == 1);
- my $key = "$n $m";
- if (exists $cache{$key}) {
- return $cache{$key};
- }
- my (@R, @D);
- $divisor_cache{$n} //= [divisors($n)];
- foreach my $d (@{$divisor_cache{$n}}) {
- if ($d >= $m) {
- push @D, $d;
- $factor_cache{$d} //= do {
- my %factors;
- @factors{factor(subint($D[-1], 1))} = ();
- [keys %factors];
- };
- }
- }
- foreach my $d (@D) {
- foreach my $p (@{$factor_cache{$d}}) {
- my $r = addint(mulint($d, subint($p, 1)), 1);
- my $k = valuation($r, $p) - 1;
- next if ($k < 1);
- my $s = powint($p, $k + 1);
- next if ("$r" ne "$s");
- my $z = powint($p, $k);
- my $u = divint($n, $d);
- my $arr = __SUB__->($u, $d);
- foreach my $v (@$arr) {
- if (modint($v, $p) != 0) {
- push @R, mulint($v, $z);
- }
- }
- }
- }
- $cache{$key} = \@R;
- }->($n, 3);
- sort { $a <=> $b } uniq(@$results);
- }
- my %tests = (
- 6 => 6187272, 10 => 196602, 11 => 8105688, 16 => 2031554,
- 25 => 1355816, 31 => 8880128, 80 => 11532, 97 => 5488,
- );
- while (my ($n, $k) = each %tests) {
- my @inverse = inverse_sigma($k);
- say "σ−¹($k) = [@inverse]";
- if (gcd(@inverse) != $n) {
- die "Error for k = $k";
- }
- }
- use Test::More;
- plan tests => 4;
- is(join(' ', inverse_sigma(42)), join(' ', 20, 26, 41));
- is(join(' ', inverse_sigma(7688)), join(' ', 2800, 2928, 4575, 7687));
- is(join(' ', inverse_sigma("110680464442257309690")), "46116860184273879040");
- is(join(' ', inverse_sigma("9325257382230393314439814176")), "3535399776779654608221686964 4302950338161146561477374638 4637009852153025247015401018 4661529533007908774933879778 4884658628787348878992283572 5187814889839710566412258045 5311639278156872382400698772 5326520187917077557965023252 5328493035801953244119300732 5495240957385767488866317781 6208298641832871739558373002 6411114450283395403677372213 6417519160023938256228496989 6454455748546107757077838269 6799666841209661791779031135 6938435552254756930386764875 6992294299511863162400845113 7215972974344947207602237095 8501184728947212952861568533 8546137477181166087378779593 9130981186767260120388984667 9214242413394317203553625829 9323102747899933426890262757 9325091641128050246166715829 9325201015147968294835238387");
- __END__
- σ−¹(6187272) = [2855646 2651676]
- σ−¹(196602) = [105650 81920]
- σ−¹(8105688) = [4953454 4947723]
- σ−¹(2031554) = [845200 999424]
- σ−¹(8880128) = [6389751 7527079]
- σ−¹(5488) = [3783 2716]
- σ−¹(11532) = [4880 4400]
- σ−¹(1355816) = [457500 390000 811875 624700]
|