I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Achilles numbers

From Rosetta Code
Task
Achilles numbers
You are encouraged to solve this task according to the task description, using any language you may know.
This page uses content from Wikipedia. The original article was at Achilles number. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)


An Achilles number is a number that is powerful but imperfect. Named after Achilles, a hero of the Trojan war, who was also powerful but imperfect.


A positive integer n is a powerful number if, for every prime factor p of n, p2 is also a divisor.

In other words, every prime factor appears at least squared in the factorization.

All Achilles numbers are powerful. However, not all powerful numbers are Achilles numbers: only those that cannot be represented as mk, where m and k are positive integers greater than 1.


A strong Achilles number is an Achilles number whose Euler totient (𝜑) is also an Achilles number.


E.G.

108 is a powerful number. Its prime factorization is 22 × 33, and thus its prime factors are 2 and 3. Both 22 = 4 and 32 = 9 are divisors of 108. However, 108 cannot be represented as mk, where m and k are positive integers greater than 1, so 108 is an Achilles number.

360 is not an Achilles number because it is not powerful. One of its prime factors is 5 but 360 is not divisible by 52 = 25.

Finally, 784 is not an Achilles number. It is a powerful number, because not only are 2 and 7 its only prime factors, but also 22 = 4 and 72 = 49 are divisors of it. Nonetheless, it is a perfect power; its square root is an even integer, so it is not an Achilles number.


500 = 22 × 53 is a strong Achilles number as its Euler totient, 𝜑(500), is 200 = 23 × 52 which is also an Achilles number.


Task
  • Find and show the first 50 Achilles numbers.
  • Find and show at least the first 20 strong Achilles numbers.
  • For at least 2 through 5, show the count of Achilles numbers with that many digits.


See also



ALGOL 68[edit]

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32
BEGIN # find Achiles Numbers: numbers whose prime factors p appear at least   #
# twice (i.e. if p is a prime factor, so is p^2) and cannot be #
# expressed as m^k for any integer m, k > 1 #
# also find strong Achiles Numbers: Achiles Numbers where the Euler's #
# totient of the number is also Achiles #
# returns the number of integers k where 1 <= k <= n that are mutually #
# prime to n #
PROC totient = ( INT n )INT: # algorithm from the second Go sample #
IF n < 3 THEN 1 # in the Totient Function task #
ELIF n = 3 THEN 2
ELSE
INT result := n;
INT v := n;
INT i := 2;
WHILE i * i <= v DO
IF v MOD i = 0 THEN
WHILE v MOD i = 0 DO v OVERAB i OD;
result -:= result OVER i
FI;
IF i = 2 THEN
i := 1
FI;
i +:= 2
OD;
IF v > 1 THEN result -:= result OVER v FI;
result
FI # totient # ;
# find the numbers #
INT max number = 1 000 000; # max number we will consider #
PR read "primes.incl.a68" PR # include prime utilities #
[]BOOL prime = PRIMESIEVE max number; # construct a sieve of primes #
# table of numbers, will be set to TRUE for the Achiles Numbers #
[ 1 : max number ]BOOL achiles;
FOR a TO UPB achiles DO
achiles[ a ] := TRUE
OD;
# remove the numbers that don't have squared primes as factors #
achiles[ 1 ] := FALSE;
FOR a TO UPB achiles DO
IF prime[ a ] THEN
# have a prime, remove it and every multiple of it that isn't a #
# multiple of a squared #
INT a count := 0;
FOR j FROM a BY a TO UPB achiles DO
a count +:= 1;
IF a count = a THEN # have a multiple of i^2, keep the number #
a count := 0
ELSE # not a multiple of i^2, remove the number #
achiles[ j ] := FALSE
FI
OD
FI
OD;
# achiles now has TRUE for the powerful numbers, remove all m^k (m,k > 1) #
FOR m FROM 2 TO UPB achiles DO
INT mk := m;
INT max mk = UPB achiles OVER m; # avoid overflow if INT is 32 bit #
WHILE mk <= max mk DO
mk *:= m;
achiles[ mk ] := FALSE
OD
OD;
# achiles now has TRUE for imperfect powerful numbers #
# show the first 50 Achiles Numbers #
BEGIN
print( ( "First 50 Achiles Numbers:", newline ) );
INT a count := 0;
FOR a WHILE a count < 50 DO
IF achiles[ a ] THEN
a count +:= 1;
print( ( " ", whole( a, -6 ) ) );
IF a count MOD 10 = 0 THEN
print( ( newline ) )
FI
FI
OD
END;
# show the first 50 Strong Achiles numbers #
BEGIN
print( ( "First 20 Strong Achiles Numbers:", newline ) );
INT s count := 0;
FOR s WHILE s count < 20 DO
IF achiles[ s ] THEN
IF achiles[ totient( s ) ] THEN
s count +:= 1;
print( ( " ", whole( s, -6 ) ) );
IF s count MOD 10 = 0 THEN
print( ( newline ) )
FI
FI
FI
OD
END;
# count the number of Achiles Numbers by their digit counts #
BEGIN
INT a count := 0;
INT power of 10 := 100;
INT digit count := 2;
FOR a TO UPB achiles DO
IF achiles[ a ] THEN
# have an Achiles Number #
a count +:= 1
FI;
IF a = power of 10 THEN
# have reached a power of 10 #
print( ( "Achiles Numbers with ", whole( digit count, 0 )
, " digits: ", whole( a count, -6 )
, newline
)
);
digit count +:= 1;
power of 10 *:= 10;
a count := 0
FI
OD
END
END
Output:
First 50 Achiles Numbers:
     72    108    200    288    392    432    500    648    675    800
    864    968    972   1125   1152   1323   1352   1372   1568   1800
   1944   2000   2312   2592   2700   2888   3087   3200   3267   3456
   3528   3872   3888   4000   4232   4500   4563   4608   5000   5292
   5324   5400   5408   5488   6075   6125   6272   6728   6912   7200
First 20 Strong Achiles Numbers:
    500    864   1944   2000   2592   3456   5000  10125  10368  12348
  12500  16875  19652  19773  30375  31104  32000  33275  37044  40500
Achiles Numbers with 2 digits:      1
Achiles Numbers with 3 digits:     12
Achiles Numbers with 4 digits:     47
Achiles Numbers with 5 digits:    192
Achiles Numbers with 6 digits:    664

C++[edit]

Translation of: Wren
Library: Boost
#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>
 
#include <boost/multiprecision/cpp_int.hpp>
 
using boost::multiprecision::uint128_t;
 
template <typename T> void unique_sort(std::vector<T>& vector) {
std::sort(vector.begin(), vector.end());
vector.erase(std::unique(vector.begin(), vector.end()), vector.end());
}
 
auto perfect_powers(uint128_t n) {
std::vector<uint128_t> result;
for (uint128_t i = 2, s = sqrt(n); i <= s; ++i)
for (uint128_t p = i * i; p < n; p *= i)
result.push_back(p);
unique_sort(result);
return result;
}
 
auto achilles(uint128_t from, uint128_t to, const std::vector<uint128_t>& pps) {
std::vector<uint128_t> result;
auto c = static_cast<uint128_t>(std::cbrt(static_cast<double>(to / 4)));
auto s = sqrt(to / 8);
for (uint128_t b = 2; b <= c; ++b) {
uint128_t b3 = b * b * b;
for (uint128_t a = 2; a <= s; ++a) {
uint128_t p = b3 * a * a;
if (p >= to)
break;
if (p >= from && !binary_search(pps.begin(), pps.end(), p))
result.push_back(p);
}
}
unique_sort(result);
return result;
}
 
uint128_t totient(uint128_t n) {
uint128_t tot = n;
if ((n & 1) == 0) {
while ((n & 1) == 0)
n >>= 1;
tot -= tot >> 1;
}
for (uint128_t p = 3; p * p <= n; p += 2) {
if (n % p == 0) {
while (n % p == 0)
n /= p;
tot -= tot / p;
}
}
if (n > 1)
tot -= tot / n;
return tot;
}
 
int main() {
auto start = std::chrono::high_resolution_clock::now();
 
const uint128_t limit = 1000000000000000;
 
auto pps = perfect_powers(limit);
auto ach = achilles(1, 1000000, pps);
 
std::cout << "First 50 Achilles numbers:\n";
for (size_t i = 0; i < 50 && i < ach.size(); ++i)
std::cout << std::setw(4) << ach[i] << ((i + 1) % 10 == 0 ? '\n' : ' ');
 
std::cout << "\nFirst 50 strong Achilles numbers:\n";
for (size_t i = 0, count = 0; count < 50 && i < ach.size(); ++i)
if (binary_search(ach.begin(), ach.end(), totient(ach[i])))
std::cout << std::setw(6) << ach[i]
<< (++count % 10 == 0 ? '\n' : ' ');
 
int digits = 2;
std::cout << "\nNumber of Achilles numbers with:\n";
for (uint128_t from = 1, to = 100; to <= limit; to *= 10, ++digits) {
size_t count = achilles(from, to, pps).size();
std::cout << std::setw(2) << digits << " digits: " << count << '\n';
from = to;
}
 
auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> duration(end - start);
std::cout << "\nElapsed time: " << duration.count() << " seconds\n";
}
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

First 50 strong Achilles numbers:
   500    864   1944   2000   2592   3456   5000  10125  10368  12348
 12500  16875  19652  19773  30375  31104  32000  33275  37044  40500
 49392  50000  52488  55296  61731  64827  67500  69984  78608  80000
 81000  83349  84375  93312 108000 111132 124416 128000 135000 148176
151875 158184 162000 165888 172872 177957 197568 200000 202612 209952

Number of Achilles numbers with:
 2 digits: 1
 3 digits: 12
 4 digits: 47
 5 digits: 192
 6 digits: 664
 7 digits: 2242
 8 digits: 7395
 9 digits: 24008
10 digits: 77330
11 digits: 247449
12 digits: 788855
13 digits: 2508051
14 digits: 7960336
15 digits: 25235383

Elapsed time: 13.2644 seconds


FreeBASIC[edit]

Function GCD(n As Uinteger, d As Uinteger) As Uinteger
Return Iif(d = 0, n, GCD(d, n Mod d))
End Function
 
Function Totient(n As Integer) As Integer
Dim As Integer m, tot = 0
For m = 1 To n
If GCD(m, n) = 1 Then tot += 1
Next m
Return tot
End Function
 
Function isPowerful(m As Integer) As Boolean
Dim As Integer n = m, f = 2, q, l = Sqr(m)
 
If m <= 1 Then Return false
Do
q = n/f
If (n Mod f) = 0 Then
If (m Mod(f*f)) Then Return false
n = q
If f > n Then Exit Do
Else
f += 1
If f > l Then
If (m Mod (n*n)) Then Return false
Exit Do
End If
End If
Loop
Return true
End Function
 
Function isAchilles(n As Integer) As Boolean
If Not isPowerful(n) Then Return false
Dim As Integer m = 2, a = m*m
Do
Do
If a = n Then Return false
If a > n Then Exit Do
a *= m
Loop
m += 1
a = m*m
Loop Until a > n
Return true
End Function
 
Dim As Integer num, n, i
Dim As Single inicio
Dim As Double t0 = Timer
 
Print "First 50 Achilles numbers:"
num = 0
n = 1
Do
If isAchilles(n) Then
Print Using "#####"; n;
num += 1
If num >= 50 Then Exit Do
If (num Mod 10) Then Print Space(3); Else Print
End If
n += 1
Loop
 
Print !"\n\nFirst 20 strong Achilles numbers:"
num = 0
n = 1
Do
If isAchilles(n) And isAchilles(Totient(n)) Then
Print Using "#####"; n;
num += 1
If num >= 20 Then Exit Do
If (num Mod 10) Then Print Space(3); Else Print
End If
n += 1
Loop
 
Print !"\n\nNumber of Achilles numbers with:"
For i = 2 To 6
inicio = Fix(10.0 ^ (i-1))
num = 0
For n = inicio To inicio*10-1
If isAchilles(n) Then num += 1
Next n
Print i; " digits:"; num
Next i
Sleep
Output:
First 50 Achilles numbers:
   72     108     200     288     392     432     500     648     675     800
  864     968     972    1125    1152    1323    1352    1372    1568    1800
 1944    2000    2312    2592    2700    2888    3087    3200    3267    3456
 3528    3872    3888    4000    4232    4500    4563    4608    5000    5292
 5324    5400    5408    5488    6075    6125    6272    6728    6912    7200
 
First 20 strong Achilles numbers:
  500     864    1944    2000    2592    3456    5000   10125   10368   12348
12500   16875   19652   19773   30375   31104   32000   33275   37044   40500

Number of Achilles numbers with:
 2 digits: 1
 3 digits: 12
 4 digits: 47
 5 digits: 192
 6 digits: 664


Go[edit]

Translation of: Wren

Based on Version 2, takes around 19 seconds.

package main
 
import (
"fmt"
"math"
"sort"
)
 
func totient(n int) int {
tot := n
i := 2
for i*i <= n {
if n%i == 0 {
for n%i == 0 {
n /= i
}
tot -= tot / i
}
if i == 2 {
i = 1
}
i += 2
}
if n > 1 {
tot -= tot / n
}
return tot
}
 
var pps = make(map[int]bool)
 
func getPerfectPowers(maxExp int) {
upper := math.Pow(10, float64(maxExp))
for i := 2; i <= int(math.Sqrt(upper)); i++ {
fi := float64(i)
p := fi
for {
p *= fi
if p >= upper {
break
}
pps[int(p)] = true
}
}
}
 
func getAchilles(minExp, maxExp int) map[int]bool {
lower := math.Pow(10, float64(minExp))
upper := math.Pow(10, float64(maxExp))
achilles := make(map[int]bool)
for b := 1; b <= int(math.Cbrt(upper)); b++ {
b3 := b * b * b
for a := 1; a <= int(math.Sqrt(upper)); a++ {
p := b3 * a * a
if p >= int(upper) {
break
}
if p >= int(lower) {
if _, ok := pps[p]; !ok {
achilles[p] = true
}
}
}
}
return achilles
}
 
func main() {
maxDigits := 15
getPerfectPowers(maxDigits)
achillesSet := getAchilles(1, 5) // enough for first 2 parts
achilles := make([]int, len(achillesSet))
i := 0
for k := range achillesSet {
achilles[i] = k
i++
}
sort.Ints(achilles)
 
fmt.Println("First 50 Achilles numbers:")
for i = 0; i < 50; i++ {
fmt.Printf("%4d ", achilles[i])
if (i+1)%10 == 0 {
fmt.Println()
}
}
 
fmt.Println("\nFirst 30 strong Achilles numbers:")
var strongAchilles []int
count := 0
for n := 0; count < 30; n++ {
tot := totient(achilles[n])
if _, ok := achillesSet[tot]; ok {
strongAchilles = append(strongAchilles, achilles[n])
count++
}
}
for i = 0; i < 30; i++ {
fmt.Printf("%5d ", strongAchilles[i])
if (i+1)%10 == 0 {
fmt.Println()
}
}
 
fmt.Println("\nNumber of Achilles numbers with:")
for d := 2; d <= maxDigits; d++ {
ac := len(getAchilles(d-1, d))
fmt.Printf("%2d digits: %d\n", d, ac)
}
}
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800 
 864  968  972 1125 1152 1323 1352 1372 1568 1800 
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 

First 30 strong Achilles numbers:
  500   864  1944  2000  2592  3456  5000 10125 10368 12348 
12500 16875 19652 19773 30375 31104 32000 33275 37044 40500 
49392 50000 52488 55296 61731 64827 67500 69984 78608 80000 

Number of Achilles numbers with:
 2 digits: 1
 3 digits: 12
 4 digits: 47
 5 digits: 192
 6 digits: 664
 7 digits: 2242
 8 digits: 7395
 9 digits: 24008
10 digits: 77330
11 digits: 247449
12 digits: 788855
13 digits: 2508051
14 digits: 7960336
15 digits: 25235383

J[edit]

Implementation:

achilles=: (*/ .>&1 * 1 = +./)@(1{__&q:)"0
strong=: [email protected](5&p:)

Task examples:

   5 10$(#~ achilles) 1+i.10000  NB. first 50 achilles numbers
72 108 200 288 392 432 500 648 675 800
864 968 972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200
 
20{.(#~ strong * achilles) 1+i.100000 NB. first twenty strong achilles numbers
500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500
 
+/achilles (+i.)/1 9*10^<:2 NB. count of two digit achilles numbers
1
+/achilles (+i.)/1 9*10^<:3
12
+/achilles (+i.)/1 9*10^<:4
47
+/achilles (+i.)/1 9*10^<:5
192
+/achilles (+i.)/1 9*10^<:6
664

Explanation of the code:

(1{__&q:) is a function which returns the non-zero powers of the prime factors of a positive integer. (__&q: returns both the primes and their factors, but here we do not care about the primes themselves.)

+./ returns the greatest common divisor of a list, and 1=+./ is true if that gcd is 1 (0 if it's false).

*/ .>&1 is true if all the values in a list are greater than 1 (0 if not).

"0 maps a function onto the individual (rank 0) items of a list or array (we use this to avoid complexities: for example if we padded our lists of prime factor powers with zeros, we could still find the gcd, but our test that the powers are greater than 1 would fail). (Actually... we could change */ .>&1 to (0 = 1 e. ]) but padding would still be a bad idea here, for performance reasons. Perhaps we ought to have an option for q: to return a sparse array...)

5&p: is euler's totient function.

(#~ predicate) list selects the elements of list where predicate is true.

Julia[edit]

using Primes
 
isAchilles(n) = (p = [x[2] for x in factor(n).pe]; all(>(1), p) && gcd(p) == 1)
 
isstrongAchilles(n) = isAchilles(n) && isAchilles(totient(n))
 
function teststrongachilles(nachilles = 50, nstrongachilles = 100)
# task 1
println("First $nachilles Achilles numbers:")
n, found = 0, 0
while found < nachilles
if isAchilles(n)
found += 1
print(rpad(n, 5), found % 10 == 0 ? "\n" : "")
end
n += 1
end
# task 2
println("\nFirst $nstrongachilles strong Achilles numbers:")
n, found = 0, 0
while found < nstrongachilles
if isstrongAchilles(n)
found += 1
print(rpad(n, 7), found % 10 == 0 ? "\n" : "")
end
n += 1
end
# task 3
println("\nCount of Achilles numbers for various intervals:")
intervals = [10:99, 100:999, 1000:9999, 10000:99999, 100000:999999]
for interval in intervals
println(lpad(interval, 15), " ", count(isAchilles, interval))
end
end
 
teststrongachilles()
 
Output:
First 50 Achilles numbers:
72   108  200  288  392  432  500  648  675  800
864  968  972  1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

First 100 strong Achilles numbers:
500    864    1944   2000   2592   3456   5000   10125  10368  12348
12500  16875  19652  19773  30375  31104  32000  33275  37044  40500  
49392  50000  52488  55296  61731  64827  67500  69984  78608  80000  
81000  83349  84375  93312  108000 111132 124416 128000 135000 148176 
151875 158184 162000 165888 172872 177957 197568 200000 202612 209952 
219488 221184 237276 243000 246924 253125 266200 270000 273375 296352 
320000 324000 333396 364500 397953 405000 432000 444528 453789 455877 
493848 497664 500000 518616 533871 540000 555579 583443 605052 607500 
629856 632736 648000 663552 665500 666792 675000 691488 740772 750141 
790272 800000 810448 820125 831875 877952 949104 972000 987696 1000188

Count of Achilles numbers for various intervals:
          10:99 1
        100:999 12
      1000:9999 47
    10000:99999 192
  100000:999999 664

Perl[edit]

Borrowed, and lightly modified, code from Powerful_numbers

Library: ntheory
use strict;
use warnings;
use feature <say current_sub>;
use experimental 'signatures';
use List::AllUtils <max head uniqint>;
use ntheory <is_square_free is_power euler_phi>;
use Math::AnyNum <:overload idiv iroot ipow is_coprime>;
 
sub table { my $t = shift() * (my $c = 1 + length max @_); ( sprintf( ('%'.$c.'d')x@_, @_) ) =~ s/.{1,$t}\K/\n/gr }
 
sub powerful_numbers ($n, $k = 2) {
my @powerful;
sub ($m, $r) {
$r < $k and push @powerful, $m and return;
for my $v (1 .. iroot(idiv($n, $m), $r)) {
if ($r > $k) { next unless is_square_free($v) and is_coprime($m, $v) }
__SUB__->($m * ipow($v, $r), $r - 1);
}
}->(1, 2*$k - 1);
sort { $a <=> $b } @powerful;
}
 
my(@P, @achilles, %Ahash, @strong);
@P = uniqint @P, powerful_numbers(10**9, $_) for 2..9; shift @P;
!is_power($_) and push @achilles, $_ and $Ahash{$_}++ for @P;
$Ahash{euler_phi $_} and push @strong, $_ for @achilles;
 
say "First 50 Achilles numbers:\n" . table 10, head 50, @achilles;
say "First 30 strong Achilles numbers:\n" . table 10, head 30, @strong;
say "Number of Achilles numbers with:\n";
for my $l (2..9) {
my $c; $l == length and $c++ for @achilles;
say "$l digits: $c";
}
Output:
First 50 Achilles numbers:
   72  108  200  288  392  432  500  648  675  800
  864  968  972 1125 1152 1323 1352 1372 1568 1800
 1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

First 30 strong Achilles numbers:
   500   864  1944  2000  2592  3456  5000 10125 10368 12348
 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500
 49392 50000 52488 55296 61731 64827 67500 69984 78608 80000

Number of Achilles numbers with:
2 digits: 1
3 digits: 12
4 digits: 47
5 digits: 192
6 digits: 664
7 digits: 2242
8 digits: 7395
9 digits: 24008

Here is a translation from Wren version 2, as an alternative.

use strict;
use warnings;
 
my %pps;
my $maxDigits = 9;
 
sub totient {
my $tot = my $n = shift;
my $i = 2;
while ($i*$i <= $n) {
unless ($n % $i) {
until($n % $i) { $n = int($n/$i) }
$tot -= int($tot/$i)
}
if ($i == 2) { $i = 1 }
$i += 2;
}
if ($n > 1) { $tot -= int($tot/$n) }
return $tot
}
 
sub getPerfectPowers {
for my $i (2..int(sqrt(my $upper = 10**( shift )))) {
my $p = $i;
while (($p *= $i) < $upper) { $pps{$p}++ }
}
}
 
sub getAchilles {
my ($lower, $upper) = map { 10** $_ } @_ ;
my %achilles = ();
my $count = 0;
for my $b (1..int($upper**(1/3))) {
my ($b3,$p) = $b * $b * $b;
for my $a (1..int(sqrt($upper))) {
last if (($p = $b3 * $a * $a) >= $upper);
$achilles{$p}++ if ($p >= $lower and !$pps{$p})
}
}
return keys %achilles
}
 
getPerfectPowers $maxDigits;
 
my @achilles = sort { $a <=> $b } getAchilles(1,5);
my %achillesSet;
@achillesSet{ @achilles } = undef;
 
print "First 50 Achilles numbers:\n";
for (0..49) { printf "%5d".($_%10 == 9 ? "\n" : " "), $achilles[$_] }
 
my %strongAchilles;
my $count = my $n = 0;
for (my $count = my $n = 0; $count < 30; $n++) {
if ( exists($achillesSet{ totient( $achilles[$n] ) })) {
$strongAchilles{ $achilles[$n] }++;
$count++
}
}
 
my @strongAchilles30 = (sort { $a <=> $b } keys %strongAchilles)[0..29];
 
print "\nFirst 30 strong Achilles numbers:\n";
for (0..29) { printf "%5d".($_%10 == 9 ? "\n" : " "), $strongAchilles30[$_] }
 
print "\nNumber of Achilles numbers with:\n";
for my $d (2..$maxDigits) {
printf "%2d digits: %d\n", $d, scalar getAchilles($d-1, $d)
}

Output is the same.

Phix[edit]

Library: Phix/online

You can run this online here, though [slightly outdated and] you should expect a blank screen for about 9s.

Translation of: Wren
with javascript_semantics
requires("1.0.2") -- [join_by(fmt)]
atom t0 = time()
constant maxDigits = iff(platform()=JS?10:12)
integer pps = new_dict()
 
procedure getPerfectPowers(integer maxExp)
    atom hi = power(10, maxExp)
    integer imax = floor(sqrt(hi))
    for i=2 to imax do
        atom p = i
        while true do
            p *= i
            if p>=hi then exit end if
            setd(p,true,pps)
        end while
    end for
end procedure

function get_achilles(integer minExp, maxExp)
    atom lo10 = power(10,minExp),
         hi10 = power(10,maxExp)
    integer bmax = floor(power(hi10,1/3)),
            amax = floor(sqrt(hi10))
    sequence achilles = {}
    for b=2 to bmax do
        atom b3 = b * b * b
        for a=2 to amax do
            atom p = b3 * a * a
            if p>=hi10 then exit end if
            if p>=lo10 then
                integer node = getd_index(p,pps)
                if node=NULL then
                    achilles &= p
                end if
            end if
        end for
    end for
    achilles = unique(achilles)
    return achilles
end function
 
getPerfectPowers(maxDigits)
sequence achilles = get_achilles(1,5)

function strong_achilles(integer n)
    integer totient = sum(sq_eq(apply(true,gcd,{tagset(n),n}),1))
    return find(totient,achilles)
end function

sequence a = join_by(achilles[1..50],1,10," ",fmt:="%4d"),
         sa = filter(achilles,strong_achilles)[1..30],
         ssa = join_by(sa,1,10," ",fmt:="%5d")
 
printf(1,"First 50 Achilles numbers:\n%s\n",{a})
printf(1,"First 30 strong Achilles numbers:\n%s\n",{ssa})
for d=2 to maxDigits do
    printf(1,"Achilles numbers with %d digits:%d\n",{d,length(get_achilles(d-1,d))})
end for
?elapsed(time()-t0)
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

First 30 strong Achilles numbers:
  500   864  1944  2000  2592  3456  5000 10125 10368 12348
12500 16875 19652 19773 30375 31104 32000 33275 37044 40500
49392 50000 52488 55296 61731 64827 67500 69984 78608 80000

Achilles numbers with 2 digits:1
Achilles numbers with 3 digits:12
Achilles numbers with 4 digits:47
Achilles numbers with 5 digits:192
Achilles numbers with 6 digits:664
Achilles numbers with 7 digits:2242
Achilles numbers with 8 digits:7395
Achilles numbers with 9 digits:24008
Achilles numbers with 10 digits:77330
Achilles numbers with 11 digits:247449
Achilles numbers with 12 digits:788855
"30.7s"

Raku[edit]

Timing is going to be system / OS dependent.

use Prime::Factor;
use Math::Root;
 
sub is-square-free (Int \n) {
constant @p = ^100 .map: { next unless .is-prime; .² };
for @p -> \p { return False if n %% p }
True
}
 
sub powerful (\n, \k = 2) {
my @p;
p(1, 2*k - 1);
sub p (\m, \r) {
@p.push(m) and return if r < k;
for 1 .. (n / m).&root(r) -> \v {
if r > k {
next unless is-square-free(v);
next unless m gcd v == 1;
}
p(m * v ** r, r - 1)
}
}
@p
}
 
my $achilles = powerful(10**9).hyper(:500batch).grep( {
my $f = .&prime-factors.Bag;
(+$f.keys > 1) && (1 == [gcd] $f.values) && (.sqrt.Int² !== $_)
} ).classify: { .chars }
 
my \𝜑 = 0, |(1..*).hyper.map: -> \t { t × [×] t.&prime-factors.squish.map: { 1 - 1/$_ } }
 
my %as = Set.new: flat $achilles.values».list;
 
my $strong = lazy (flat $achilles.sort».value».list».sort).grep: { ?%as{𝜑[$_]} };
 
put "First 50 Achilles numbers:";
put (flat $achilles.sort».value».list».sort)[^50].batch(10)».fmt("%4d").join("\n");
 
put "\nFirst 30 strong Achilles numbers:";
put $strong[^30].batch(10)».fmt("%5d").join("\n");
 
put "\nNumber of Achilles numbers with:";
say "$_ digits: " ~ +$achilles{$_} // 0 for 2..9;
 
printf "\n%.1f total elapsed seconds\n", now - INIT now;
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

First 30 strong Achilles numbers:
  500   864  1944  2000  2592  3456  5000 10125 10368 12348
12500 16875 19652 19773 30375 31104 32000 33275 37044 40500
49392 50000 52488 55296 61731 64827 67500 69984 78608 80000

Number of Achilles numbers with:
2 digits: 1
3 digits: 12
4 digits: 47
5 digits: 192
6 digits: 664
7 digits: 2242
8 digits: 7395
9 digits: 24008

6.1 total elapsed seconds

Could go further but slows to a crawl and starts chewing up memory in short order.

10 digits: 77330
11 digits: 247449
12 digits: 788855

410.4 total elapsed seconds

Rust[edit]

Translation of: Wren
fn perfect_powers(n: u128) -> Vec<u128> {
let mut powers = Vec::<u128>::new();
let sqrt = (n as f64).sqrt() as u128;
for i in 2..=sqrt {
let mut p = i * i;
while p < n {
powers.push(p);
p *= i;
}
}
powers.sort();
powers.dedup();
powers
}
 
fn bsearch<T: Ord>(vector: &Vec<T>, value: &T) -> bool {
match vector.binary_search(value) {
Ok(_) => true,
_ => false,
}
}
 
fn achilles(from: u128, to: u128, pps: &Vec<u128>) -> Vec<u128> {
let mut result = Vec::<u128>::new();
let cbrt = ((to / 4) as f64).cbrt() as u128;
let sqrt = ((to / 8) as f64).sqrt() as u128;
for b in 2..=cbrt {
let b3 = b * b * b;
for a in 2..=sqrt {
let p = b3 * a * a;
if p >= to {
break;
}
if p >= from && !bsearch(&pps, &p) {
result.push(p);
}
}
}
result.sort();
result.dedup();
result
}
 
fn totient(mut n: u128) -> u128 {
let mut tot = n;
if (n & 1) == 0 {
while (n & 1) == 0 {
n >>= 1;
}
tot -= tot >> 1;
}
let mut p = 3;
while p * p <= n {
if n % p == 0 {
while n % p == 0 {
n /= p;
}
tot -= tot / p;
}
p += 2;
}
if n > 1 {
tot -= tot / n;
}
tot
}
 
fn main() {
use std::time::Instant;
let t0 = Instant::now();
let limit = 1000000000000000u128;
 
let pps = perfect_powers(limit);
let ach = achilles(1, 1000000, &pps);
 
println!("First 50 Achilles numbers:");
for i in 0..50 {
print!("{:4}{}", ach[i], if (i + 1) % 10 == 0 { "\n" } else { " " });
}
 
println!("\nFirst 50 strong Achilles numbers:");
for (i, n) in ach
.iter()
.filter(|&x| bsearch(&ach, &totient(*x)))
.take(50)
.enumerate()
{
print!("{:6}{}", n, if (i + 1) % 10 == 0 { "\n" } else { " " });
}
println!();
 
let mut from = 1u128;
let mut to = 100u128;
let mut digits = 2;
while to <= limit {
let count = achilles(from, to, &pps).len();
println!("{:2} digits: {}", digits, count);
from = to;
to *= 10;
digits += 1;
}
 
let duration = t0.elapsed();
println!("\nElapsed time: {} milliseconds", duration.as_millis());
}
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

First 50 strong Achilles numbers:
   500    864   1944   2000   2592   3456   5000  10125  10368  12348
 12500  16875  19652  19773  30375  31104  32000  33275  37044  40500
 49392  50000  52488  55296  61731  64827  67500  69984  78608  80000
 81000  83349  84375  93312 108000 111132 124416 128000 135000 148176
151875 158184 162000 165888 172872 177957 197568 200000 202612 209952

 2 digits: 1
 3 digits: 12
 4 digits: 47
 5 digits: 192
 6 digits: 664
 7 digits: 2242
 8 digits: 7395
 9 digits: 24008
10 digits: 77330
11 digits: 247449
12 digits: 788855
13 digits: 2508051
14 digits: 7960336
15 digits: 25235383

Elapsed time: 12608 milliseconds

Wren[edit]

Library: Wren-math
Library: Wren-seq
Library: Wren-fmt

Version 1 (Brute force)[edit]

This finds the number of 6 digit Achilles numbers in 2.5 seconds, 7 digits in 51 seconds but 8 digits needs a whopping 21 minutes!

import "./math" for Int
import "./seq" for Lst
import "./fmt" for Fmt
 
var maxDigits = 8
var limit = 10.pow(maxDigits)
var c = Int.primeSieve(limit-1, false)
 
var totient = Fn.new { |n|
var tot = n
var i = 2
while (i*i <= n) {
if (n%i == 0) {
while(n%i == 0) n = (n/i).floor
tot = tot - (tot/i).floor
}
if (i == 2) i = 1
i = i + 2
}
if (n > 1) tot = tot - (tot/n).floor
return tot
}
 
var isPerfectPower = Fn.new { |n|
if (n == 1) return true
var x = 2
while (x * x <= n) {
var y = 2
var p = x.pow(y)
while (p > 0 && p <= n) {
if (p == n) return true
y = y + 1
p = x.pow(y)
}
x = x + 1
}
return false
}
 
var isPowerful = Fn.new { |n|
while (n % 2 == 0) {
var p = 0
while (n % 2 == 0) {
n = (n/2).floor
p = p + 1
}
if (p == 1) return false
}
var f = 3
while (f * f <= n) {
var p = 0
while (n % f == 0) {
n = (n/f).floor
p = p + 1
}
if (p == 1) return false
f = f + 2
}
return n == 1
}
 
var isAchilles = Fn.new { |n| c[n] && isPowerful.call(n) && !isPerfectPower.call(n) }
 
var isStrongAchilles = Fn.new { |n|
if (!isAchilles.call(n)) return false
var tot = totient.call(n)
return isAchilles.call(tot)
}
 
System.print("First 50 Achilles numbers:")
var achilles = []
var count = 0
var n = 2
while (count < 50) {
if (isAchilles.call(n)) {
achilles.add(n)
count = count + 1
}
n = n + 1
}
for (chunk in Lst.chunks(achilles, 10)) Fmt.print("$4d", chunk)
 
System.print("\nFirst 30 strong Achilles numbers:")
var strongAchilles = []
count = 0
n = achilles[0]
while (count < 30) {
if (isStrongAchilles.call(n)) {
strongAchilles.add(n)
count = count + 1
}
n = n + 1
}
for (chunk in Lst.chunks(strongAchilles, 10)) Fmt.print("$5d", chunk)
 
System.print("\nNumber of Achilles numbers with:")
var pow = 10
for (i in 2..maxDigits) {
var count = 0
for (j in pow..pow*10-1) {
if (isAchilles.call(j)) count = count + 1
}
System.print("%(i) digits: %(count)")
pow = pow * 10
}
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

First 30 strong Achilles numbers:
  500   864  1944  2000  2592  3456  5000 10125 10368 12348
12500 16875 19652 19773 30375 31104 32000 33275 37044 40500
49392 50000 52488 55296 61731 64827 67500 69984 78608 80000

Number of Achilles numbers with:
2 digits: 1
3 digits: 12
4 digits: 47
5 digits: 192
6 digits: 664
7 digits: 2242
8 digits: 7395

Version 2 (Much faster)[edit]

Library: Wren-set

Here we make use of the fact that powerful numbers are always of the form a²b³, where a and b > 0, to generate such numbers up to a given limit. We also generate in advance all perfect powers up to the same limit.

Ridiculously fast compared to the previous method: 12 digits can now be reached in 1.03 seconds, 13 digits in 3.7 seconds, 14 digits in 12.2 seconds and 15 digits in 69 seconds.

import "./set" for Set
import "./seq" for Lst
import "./fmt" for Fmt
 
var totient = Fn.new { |n|
var tot = n
var i = 2
while (i*i <= n) {
if (n%i == 0) {
while(n%i == 0) n = (n/i).floor
tot = tot - (tot/i).floor
}
if (i == 2) i = 1
i = i + 2
}
if (n > 1) tot = tot - (tot/n).floor
return tot
}
 
var pps = Set.new()
 
var getPerfectPowers = Fn.new { |maxExp|
var upper = 10.pow(maxExp)
for (i in 2..upper.sqrt.floor) {
var p = i
while ((p = p * i) < upper) pps.add(p)
}
}
 
var getAchilles = Fn.new { |minExp, maxExp|
var lower = 10.pow(minExp)
var upper = 10.pow(maxExp)
var achilles = Set.new() // avoids duplicates
for (b in 1..upper.cbrt.floor) {
var b3 = b * b * b
for (a in 1..upper.sqrt.floor) {
var p = b3 * a * a
if (p >= upper) break
if (p >= lower) {
if (!pps.contains(p)) achilles.add(p)
}
}
}
return achilles
}
 
var maxDigits = 15
getPerfectPowers.call(maxDigits)
 
var achillesSet = getAchilles.call(1, 5) // enough for first 2 parts
var achilles = achillesSet.toList
achilles.sort()
 
System.print("First 50 Achilles numbers:")
for (chunk in Lst.chunks(achilles[0..49], 10)) Fmt.print("$4d", chunk)
 
System.print("\nFirst 30 strong Achilles numbers:")
var strongAchilles = []
var count = 0
var n = 0
while (count < 30) {
var tot = totient.call(achilles[n])
if (achillesSet.contains(tot)) {
strongAchilles.add(achilles[n])
count = count + 1
}
n = n + 1
}
for (chunk in Lst.chunks(strongAchilles, 10)) Fmt.print("$5d", chunk)
 
System.print("\nNumber of Achilles numbers with:")
for (d in 2..maxDigits) {
var ac = getAchilles.call(d-1, d).count
Fmt.print("$2d digits: $d", d, ac)
}
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

First 30 strong Achilles numbers:
  500   864  1944  2000  2592  3456  5000 10125 10368 12348
12500 16875 19652 19773 30375 31104 32000 33275 37044 40500
49392 50000 52488 55296 61731 64827 67500 69984 78608 80000

Number of Achilles numbers with:
 2 digits: 1
 3 digits: 12
 4 digits: 47
 5 digits: 192
 6 digits: 664
 7 digits: 2242
 8 digits: 7395
 9 digits: 24008
10 digits: 77330
11 digits: 247449
12 digits: 788855
13 digits: 2508051
14 digits: 7960336
15 digits: 25235383

XPL0[edit]

func GCD(N, D);         \Return the greatest common divisor of N and D
int N, D; \numerator and denominator
int R;
[if D > N then
[R:= D; D:= N; N:= R]; \swap D and N
while D > 0 do
[R:= rem(N/D);
N:= D;
D:= R;
];
return N;
]; \GCD
 
func Totient(N); \Return the totient of N
int N, Phi, M;
[Phi:= 0;
for M:= 1 to N do
if GCD(M, N) = 1 then Phi:= Phi+1;
return Phi;
];
 
func Powerful(N0); \Return 'true' if N0 is a powerful number
int N0, N, F, Q, L;
[if N0 <= 1 then return false;
N:= N0; F:= 2;
L:= sqrt(N0);
loop [Q:= N/F;
if rem(0) = 0 then \found a factor
[if rem(N0/(F*F)) then return false;
N:= Q;
if F>N then quit;
]
else [F:= F+1;
if F > L then
[if rem(N0/(N*N)) then return false;
quit;
];
];
];
return true;
];
 
func Achilles(N); \Return 'true' if N is an Achilles number
int N, M, A;
[if not Powerful(N) then return false;
M:= 2;
A:= M*M;
repeat loop [if A = N then return false;
if A > N then quit;
A:= A*M;
];
M:= M+1;
A:= M*M;
until A > N;
return true;
];
 
int Cnt, N, Pwr, Start;
[Cnt:= 0;
N:= 1;
loop [if Achilles(N) then
[IntOut(0, N);
Cnt:= Cnt+1;
if Cnt >= 50 then quit;
if rem(Cnt/10) then ChOut(0, 9) else CrLf(0);
];
N:= N+1;
];
CrLf(0); CrLf(0);
Cnt:= 0;
N:= 1;
loop [if Achilles(N) then
if Achilles(Totient(N)) then
[IntOut(0, N);
Cnt:= Cnt+1;
if Cnt >= 20 then quit;
if rem(Cnt/10) then ChOut(0, 9) else CrLf(0);
];
N:= N+1;
];
CrLf(0); CrLf(0);
for Pwr:= 1 to 6 do
[IntOut(0, Pwr); Text(0, ": ");
Start:= fix(Pow(10.0, float(Pwr-1)));
Cnt:= 0;
for N:= Start to Start*10-1 do
if Achilles(N) then Cnt:= Cnt+1;
IntOut(0, Cnt); CrLf(0);
];
]
Output:
72      108     200     288     392     432     500     648     675     800
864     968     972     1125    1152    1323    1352    1372    1568    1800
1944    2000    2312    2592    2700    2888    3087    3200    3267    3456
3528    3872    3888    4000    4232    4500    4563    4608    5000    5292
5324    5400    5408    5488    6075    6125    6272    6728    6912    7200

500     864     1944    2000    2592    3456    5000    10125   10368   12348
12500   16875   19652   19773   30375   31104   32000   33275   37044   40500

1: 0
2: 1
3: 12
4: 47
5: 192
6: 664