123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147 |
- #!/usr/bin/perl
- use 5.020;
- use warnings;
- use experimental qw(signatures);
- use Math::GMPz;
- use ntheory qw(:all);
- sub check_valuation ($n, $p) {
- #~ if ($p == 2) {
- #~ return valuation($n, $p) < 5;
- #~ }
- #~ if ($p == 3) {
- #~ return valuation($n, $p) < 3;
- #~ }
- #~ if ($p == 7) {
- #~ return valuation($n, $p) < 3;
- #~ }
- ($n % $p) != 0;
- }
- sub smooth_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 arithmetic_derivative {
- my ($n) = @_;
- my @factors = factor($n);
- my $sum = 0;
- foreach my $p (@factors) {
- $sum += divint($n, $p);
- }
- return $sum;
- }
- 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) {
- dynamicPreimage($N, cook_sigma($N, $k));
- }
- sub isok ($n) {
- my $d = arithmetic_derivative($n);
- foreach my $k (inverse_sigma($n)) {
- if ($d == $k) {
- return $k;
- }
- }
- return undef;
- }
- my @smooth_primes;
- foreach my $p (@{primes(1000)}) {
- if ($p == 2) {
- push @smooth_primes, $p;
- next;
- }
- if (is_smooth($p+1, 3)) {
- push @smooth_primes, $p;
- }
- }
- my $h = smooth_numbers(1e15, \@smooth_primes);
- say "\nFound: ", scalar(@$h), " terms\n";
- my %table;
- foreach my $n (@$h) {
- if (defined(my $k = isok($n))) {
- say "Found: $k -> ", join(' * ', map { "$_->[0]^$_->[1]" } factor_exp($k));
- }
- }
|