``````#!/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";
``````