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
91
92
93
94
95
96
97
98
99
|
#!/usr/bin/env perl
use warnings;
use strict;
use Math::Prime::Util qw/:all/;
use Math::GMPz;
# Sum of all n where is_power(euler_phi(n^2),3) = 1
# Simple but very slow way. The brute force method later in this file is
# basically the same thing, but using the more efficient ranged moebius and
# totient calls over intervals.
#
# Pari:
# s=0; for(n=2,limit,if(ispower(n*eulerphi(n),3),s=s+n)); print(s)
# Perl/MPU:
# my $s=0;
# for my $n (2..$limit) { $s += $n if is_power($n*euler_phi($n),3); }
# say $s;
#
# TIMING:
# 10^7 2*10^7 10^8 10^10
# Clever 0.06s 0.09s 0.24s 5s
# Brute 5.0s 10.2s 52.9s 5 hours
# Simple MPU 10.8s 24.6s 159s 1 day?
# Simple Pari 13.6s 33.4s 277s 5 days?
#
my $limit = shift || 10**10-1;
my $method = lc(shift || 'clever');
die "Method must be 'clever' or 'brute'\n" unless $method =~ /^(clever|brute)$/;
my $sum = 0;
if ($method eq 'clever') {
# About 5 seconds for 10^10-1
my $cblimit = int( ($limit*$limit) ** 0.3334 + 0.01 );
foreach my $k (2 .. $cblimit) {
next if $k & 1;
my($p, $e) = @{ (factor_exp($k))[-1] };
$e *= 3;
next unless $e & 1;
my $m = int($k / ($p ** int($e/3)));
$m **= 3;
next if $m % ($p-1);
$m = int($m / ($p-1));
my $n = $p ** (($e+1) >> 1);
next if $n >= $limit;
while ($m > 1) {
my ($p,$e) = @{ (factor_exp($m))[-1] };
last unless $e & 1;
last if $m % ($p-1);
$n *= $p ** (($e+1) >> 1);
last if $n >= $limit;
$m = int($m / ( ($p-1) * ($p**$e) ) );
}
if ($m == 1) {
#print "$n\n";
$sum += $n;
}
}
} else {
# About 5 hours for 10^10-1
my $interval = 10_000_000; # Window size for moebius/totient
#prime_precalc(10**9); # Slightly faster ranged phi
my($beg,$end) = (0,0);
while ($beg < $limit) {
$end = $beg + $interval - 1;
$end = $limit if $end > $limit;
my $start = ($beg<2)?2:$beg;
my $glim = int(~0 / $end);
my @m = moebius($beg, $end);
my @t = euler_phi($beg, $end);
if ($end <= $glim) { # Totient($n) * $n will always be < ~0
foreach my $n ($start .. $end) {
next unless $m[$n-$beg] == 0;
my $totn2 = $n * $t[$n-$beg];
if (is_power($totn2,3)) {
# print "$n\n";
$sum += $n
}
}
} else {
foreach my $n ($start .. $end) {
next unless $m[$n-$beg] == 0;
my $tot = $t[$n-$beg];
if ($tot <= $glim) {
print "$n\n" if is_power($n * $tot, 3);
} else {
$tot = Math::GMPz->new($n) * $tot;
print "$n\n" if Math::GMPz::Rmpz_perfect_power_p($tot) && is_power($tot,3);
}
}
}
$beg = $end+1;
}
}
print "$sum\n";
|
