-
Notifications
You must be signed in to change notification settings - Fork 0
/
smooth_generate.pl
90 lines (59 loc) · 1.62 KB
/
smooth_generate.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#!/usr/bin/perl
# Composite numbers for which the harmonic mean of proper divisors is an integer.
# https://oeis.org/A247077
# Known terms:
# 1645, 88473, 63626653506
# These are composite numbers n such that sigma(n)-1 divides n*(tau(n)-1).
# Conjecture: all terms are of the form n*(sigma(n)-1) where sigma(n)-1 is prime. - Chai Wah Wu, Dec 15 2020
# If the above conjecture is true, then a(4) > 10^14.
# This program assumes that the above conjecture is true.
use 5.020;
use warnings;
use experimental qw(signatures);
use Math::GMPz;
use ntheory qw(:all);
sub check_valuation ($n, $p) {
1
#~ 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 isok ($m) {
modint(mulint($m, divisor_sum($m, 0) - 1), divisor_sum($m) - 1) == 0;
}
my $h = smooth_numbers(1e13, primes(20));
say "\nTotal: ", scalar(@$h), " terms\n";
my %table;
foreach my $n (@$h) {
#$n > 1e7 || next;
my $p = divisor_sum($n) - 1;
is_prime($p) || next;
my $m = mulint($n, $p);
if (isok($m)) {
say "Found: $n -> $m";
}
}
__END__
Found: 35 -> 1645
Found: 231 -> 88473
Found: 171366 -> 63626653506
Found: 3662109375 -> 22351741783447265625