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
<lang algol68>BEGIN # find Achilles 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 Achilles Numbers: Achilles Numbers where the Euler's #
# totient of the number is also Achilles #
# 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 Achilles 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 Achilles Numbers #
BEGIN
print( ( "First 50 Achilles 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 Achilles numbers #
BEGIN
print( ( "First 20 Strong Achilles 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 Achilles 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 Achilles Number #
a count +:= 1
FI;
IF a = power of 10 THEN
# have reached a power of 10 #
print( ( "Achilles 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</lang>
- 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
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
ARM Assembly
<lang ARM Assembly> /* ARM assembly Raspberry PI */ /* program achilleNumber.s */
/* REMARK 1 : this program use routines in a include file see task Include a file language arm assembly for the routine affichageMess conversion10 see at end of this program the instruction include */
/* for constantes see task include a file in arm assembly */ /************************************/ /* Constantes */ /************************************/ .include "../constantes.inc" .equ NBFACT, 33 .equ MAXI, 50 .equ MAXI1, 20 .equ MAXI2, 1000000
/*********************************/ /* Initialized data */ /*********************************/ .data szMessNumber: .asciz " @ " szCarriageReturn: .asciz "\n" szErrorGen: .asciz "Program error !!!\n" szMessPrime: .asciz "This number is prime.\n" szMessTitAchille: .asciz "First 50 Achilles Numbers:\n" szMessTitStrong: .asciz "First 20 Strong Achilles Numbers:\n" szMessDigitsCounter: .asciz "Numbers with @ digits : @ \n" /*********************************/ /* UnInitialized data */ /*********************************/ .bss sZoneConv: .skip 24 tbZoneDecom: .skip 8 * NBFACT // factor 4 bytes, number of each factor 4 bytes /*********************************/ /* code section */ /*********************************/ .text .global main main: @ entry of program
ldr r0,iAdrszMessTitAchille bl affichageMess mov r4,#1 @ start number mov r5,#0 @ total counter mov r6,#0 @ line display counter
1:
mov r0,r4 bl controlAchille cmp r0,#0 @ achille number ? beq 2f @ no mov r0,r4 ldr r1,iAdrsZoneConv bl conversion10 @ call décimal conversion ldr r0,iAdrszMessNumber ldr r1,iAdrsZoneConv @ insert conversion in message bl strInsertAtCharInc bl affichageMess @ display message add r5,r5,#1 @ increment counter add r6,r6,#1 @ increment indice line display cmp r6,#10 @ if = 10 new line bne 2f mov r6,#0 ldr r0,iAdrszCarriageReturn bl affichageMess
2:
add r4,r4,#1 @ increment number cmp r5,#MAXI blt 1b @ and loop ldr r0,iAdrszMessTitStrong bl affichageMess mov r4,#1 @ start number mov r5,#0 @ total counter mov r6,#0
3:
mov r0,r4 bl controlAchille cmp r0,#0 beq 4f mov r0,r4 bl computeTotient bl controlAchille cmp r0,#0 beq 4f mov r0,r4 ldr r1,iAdrsZoneConv bl conversion10 @ call décimal conversion ldr r0,iAdrszMessNumber ldr r1,iAdrsZoneConv @ insert conversion in message bl strInsertAtCharInc bl affichageMess @ display message add r5,r5,#1 add r6,r6,#1 cmp r6,#10 bne 4f mov r6,#0 ldr r0,iAdrszCarriageReturn bl affichageMess
4:
add r4,r4,#1 cmp r5,#MAXI1 blt 3b ldr r3,icstMaxi2 mov r4,#1 @ start number mov r6,#0 @ total counter 2 digits mov r7,#0 @ total counter 3 digits mov r8,#0 @ total counter 4 digits mov r9,#0 @ total counter 5 digits mov r10,#0 @ total counter 6 digits
5:
mov r0,r4 bl controlAchille cmp r0,#0 beq 6f mov r0,r4 ldr r1,iAdrsZoneConv bl conversion10 @ call décimal conversion r0 return digit number cmp r0,#6 addeq r10,r10,#1 beq 6f cmp r0,#5 addeq r9,r9,#1 beq 6f cmp r0,#4 addeq r8,r8,#1 beq 6f cmp r0,#3 addeq r7,r7,#1 beq 6f cmp r0,#2 addeq r6,r6,#1 beq 6f
6:
add r4,r4,#1 cmp r4,r3 blt 5b mov r0,#2 mov r1,r6 bl displayCounter mov r0,#3 mov r1,r7 bl displayCounter mov r0,#4 mov r1,r8 bl displayCounter mov r0,#5 mov r1,r9 bl displayCounter mov r0,#6 mov r1,r10 bl displayCounter b 100f
98:
ldr r0,iAdrszErrorGen bl affichageMess
100: @ standard end of the program
mov r0, #0 @ return code mov r7, #EXIT @ request to exit program svc #0 @ perform the system call
iAdrszCarriageReturn: .int szCarriageReturn iAdrszErrorGen: .int szErrorGen iAdrsZoneConv: .int sZoneConv iAdrtbZoneDecom: .int tbZoneDecom iAdrszMessNumber: .int szMessNumber iAdrszMessTitAchille: .int szMessTitAchille iAdrszMessTitStrong: .int szMessTitStrong icstMaxi2: .int MAXI2 /******************************************************************/ /* display digit counter */ /******************************************************************/ /* r0 contains limit */ /* r1 contains counter */ displayCounter:
push {r1-r3,lr} @ save registers
mov r2,r1
ldr r1,iAdrsZoneConv
bl conversion10 @ call décimal conversion
ldr r0,iAdrszMessDigitsCounter
ldr r1,iAdrsZoneConv @ insert conversion in message
bl strInsertAtCharInc
mov r3,r0
mov r0,r2
ldr r1,iAdrsZoneConv
bl conversion10 @ call décimal conversion
mov r0,r3
ldr r1,iAdrsZoneConv @ insert conversion in message
bl strInsertAtCharInc
bl affichageMess @ display message
100:
pop {r1-r3,pc} @ restaur registers
iAdrszMessDigitsCounter: .int szMessDigitsCounter /******************************************************************/ /* control if number is Achille number */ /******************************************************************/ /* r0 contains number */ /* r0 return 0 if not else return 1 */ controlAchille:
push {r1-r4,lr} @ save registers
mov r4,r0
ldr r1,iAdrtbZoneDecom
bl decompFact @ factor decomposition
cmp r0,#-1
beq 98f @ error ?
cmp r0,#1 @ one only factor ?
moveq r0,#0
beq 100f
mov r1,r0
ldr r0,iAdrtbZoneDecom
mov r2,r4
bl controlDivisor
b 100f
98:
ldr r0,iAdrszErrorGen bl affichageMess
100:
pop {r1-r4,pc} @ restaur registers
/******************************************************************/ /* control divisors function */ /******************************************************************/ /* r0 contains address of divisors area */ /* r1 contains the number of area items */ /* r2 contains number */ controlDivisor:
push {r1-r10,lr} @ save registers
cmp r1,#0
moveq r0,#0
beq 100f
mov r6,r1 @ factors number
mov r8,r2 @ save number
mov r9,#0 @ indice
mov r4,r0 @ save area address
add r5,r4,r9,lsl #3 @ compute address first factor
ldr r7,[r5,#4] @ load first exposant of factor
add r2,r9,#1
1:
add r5,r4,r2,lsl #3 @ compute address next factor ldr r3,[r5,#4] @ load exposant of factor cmp r3,r7 @ factor exposant <> ? bne 2f @ yes -> end verif add r2,r2,#1 @ increment indice cmp r2,r6 @ factor maxi ? blt 1b @ no -> loop mov r0,#0 b 100f @ all exposants are equals
2:
mov r10,r2 @ save indice
21:
movlt r2,r7 @ if r3 < r7 -> inversion movlt r7,r3 movlt r3,r2 @ r7 is the smaller exposant mov r0,r3 mov r1,r7 @ r7 < r3 bl computePgcd cmp r0,#1 beq 23f @ no commun multiple -> ne peux donc pas etre une puissance
22:
add r10,r10,#1 @ increment indice cmp r10,r6 @ factor maxi ? movge r0,#0 bge 100f @ yes -> all exposants are multiples to smaller add r5,r4,r10,lsl #3 ldr r3,[r5,#4] @ load exposant of next factor cmp r3,r7 beq 22b @ for next b 21b @ for compare the 2 exposants
23:
mov r9,#0 @ indice
3:
add r5,r4,r9,lsl #3 ldr r7,[r5] @ load factor mul r1,r7,r7 @ factor square mov r0,r8 @ number bl division cmp r3,#0 @ remainder null ? movne r0,#0 bne 100f add r9,#1 @ other factor cmp r9,r6 @ factors maxi ? blt 3b mov r0,#1 @ achille number ok
100:
pop {r1-r10,lr} @ restaur registers
bx lr @ return
/******************************************/ /* calcul du pgcd */ /*****************************************/ /* r0 number one */ /* r1 number two */ /* r0 result return */ computePgcd:
push {r2,lr} @ save registers
1:
cmp r0,#0 ble 2f cmp r1,r0 movgt r2,r0 movgt r0,r1 movgt r1,r2 sub r0,r1 b 1b
2:
mov r0,r1
pop {r2,pc} @ restaur registers
/******************************************************************/ /* compute totient of number */ /******************************************************************/ /* r0 contains number */ computeTotient:
push {r1-r5,lr} @ save registers
mov r4,r0 @ totient
mov r5,r0 @ save number
mov r1,#0 @ for first divisor
1: @ begin loop
mul r3,r1,r1 @ compute square cmp r3,r5 @ compare number bgt 4f @ end add r1,r1,#2 @ next divisor mov r0,r5 bl division cmp r3,#0 @ remainder null ? bne 3f
2: @ begin loop 2
mov r0,r5 bl division cmp r3,#0 moveq r5,r2 @ new value = quotient beq 2b mov r0,r4 @ totient bl division sub r4,r4,r2 @ compute new totient
3:
cmp r1,#2 @ first divisor ? moveq r1,#1 @ divisor = 1 b 1b @ and loop
4:
cmp r5,#1 @ final value > 1 ble 5f mov r0,r4 @ totient mov r1,r5 @ divide by value bl division sub r4,r4,r2 @ compute new totient
5:
mov r0,r4
100:
pop {r1-r5,pc} @ restaur registers
/******************************************************************/ /* factor decomposition */ /******************************************************************/ /* r0 contains number */ /* r1 contains address of divisors area */ /* r0 return divisors items in table */ decompFact:
push {r1-r8,lr} @ save registers
mov r5,r1
mov r8,r0 @ save number
bl isPrime @ prime ?
cmp r0,#1
beq 98f @ yes is prime
mov r4,#0 @ raz indice
mov r1,#2 @ first divisor
mov r6,#0 @ previous divisor
mov r7,#0 @ number of same divisors
2:
mov r0,r8 @ dividende bl division @ r1 divisor r2 quotient r3 remainder cmp r3,#0 bne 5f @ if remainder <> zero -> no divisor mov r8,r2 @ else quotient -> new dividende cmp r1,r6 @ same divisor ? beq 4f @ yes cmp r6,#0 @ no but is the first divisor ? beq 3f @ yes str r6,[r5,r4,lsl #2] @ else store in the table add r4,r4,#1 @ and increment counter str r7,[r5,r4,lsl #2] @ store counter add r4,r4,#1 @ next item mov r7,#0 @ and raz counter
3:
mov r6,r1 @ new divisor
4:
add r7,r7,#1 @ increment counter b 7f @ and loop /* not divisor -> increment next divisor */
5:
cmp r1,#2 @ if divisor = 2 -> add 1 addeq r1,#1 addne r1,#2 @ else add 2 b 2b /* divisor -> test if new dividende is prime */
7:
mov r3,r1 @ save divisor cmp r8,#1 @ dividende = 1 ? -> end beq 10f mov r0,r8 @ new dividende is prime ? mov r1,#0 bl isPrime @ the new dividende is prime ? cmp r0,#1 bne 10f @ the new dividende is not prime
cmp r8,r6 @ else dividende is same divisor ? beq 9f @ yes cmp r6,#0 @ no but is the first divisor ? beq 8f @ yes it is a first str r6,[r5,r4,lsl #2] @ else store in table add r4,r4,#1 @ and increment counter str r7,[r5,r4,lsl #2] @ and store counter add r4,r4,#1 @ next item
8:
mov r6,r8 @ new dividende -> divisor prec mov r7,#0 @ and raz counter
9:
add r7,r7,#1 @ increment counter b 11f
10:
mov r1,r3 @ current divisor = new divisor cmp r1,r8 @ current divisor > new dividende ? ble 2b @ no -> loop /* end decomposition */
11:
str r6,[r5,r4,lsl #2] @ store last divisor add r4,r4,#1 str r7,[r5,r4,lsl #2] @ and store last number of same divisors add r4,r4,#1 lsr r0,r4,#1 @ return number of table items mov r3,#0 str r3,[r5,r4,lsl #2] @ store zéro in last table item add r4,r4,#1 str r3,[r5,r4,lsl #2] @ and zero in counter same divisor b 100f
98:
//ldr r0,iAdrszMessPrime //bl affichageMess mov r0,#1 @ return code b 100f
99:
ldr r0,iAdrszErrorGen bl affichageMess mov r0,#-1 @ error code b 100f
100:
pop {r1-r8,lr} @ restaur registers
bx lr
iAdrszMessPrime: .int szMessPrime
/***************************************************/ /* check if a number is prime */ /***************************************************/ /* r0 contains the number */ /* r0 return 1 if prime 0 else */ @2147483647 @4294967297 @131071 isPrime:
push {r1-r6,lr} @ save registers
cmp r0,#0
beq 90f
cmp r0,#17
bhi 1f
cmp r0,#3
bls 80f @ for 1,2,3 return prime
cmp r0,#5
beq 80f @ for 5 return prime
cmp r0,#7
beq 80f @ for 7 return prime
cmp r0,#11
beq 80f @ for 11 return prime
cmp r0,#13
beq 80f @ for 13 return prime
cmp r0,#17
beq 80f @ for 17 return prime
1:
tst r0,#1 @ even ? beq 90f @ yes -> not prime mov r2,r0 @ save number sub r1,r0,#1 @ exposant n - 1 mov r0,#3 @ base bl moduloPuR32 @ compute base power n - 1 modulo n cmp r0,#1 bne 90f @ if <> 1 -> not prime mov r0,#5 bl moduloPuR32 cmp r0,#1 bne 90f mov r0,#7 bl moduloPuR32 cmp r0,#1 bne 90f mov r0,#11 bl moduloPuR32 cmp r0,#1 bne 90f mov r0,#13 bl moduloPuR32 cmp r0,#1 bne 90f mov r0,#17 bl moduloPuR32 cmp r0,#1 bne 90f
80:
mov r0,#1 @ is prime b 100f
90:
mov r0,#0 @ no prime
100: @ fin standard de la fonction
pop {r1-r6,lr} @ restaur des registres
bx lr @ retour de la fonction en utilisant lr
/********************************************************/ /* Calcul modulo de b puissance e modulo m */ /* Exemple 4 puissance 13 modulo 497 = 445 */ /* */ /********************************************************/ /* r0 nombre */ /* r1 exposant */ /* r2 modulo */ /* r0 return result */ moduloPuR32:
push {r1-r7,lr} @ save registers
cmp r0,#0 @ verif <> zero
beq 100f
cmp r2,#0 @ verif <> zero
beq 100f @ TODO: vérifier les cas erreur
1:
mov r4,r2 @ save modulo mov r5,r1 @ save exposant mov r6,r0 @ save base mov r3,#1 @ start result
mov r1,#0 @ division de r0,r1 par r2 bl division32R mov r6,r2 @ base <- remainder
2:
tst r5,#1 @ exposant even or odd beq 3f umull r0,r1,r6,r3 mov r2,r4 bl division32R mov r3,r2 @ result <- remainder
3:
umull r0,r1,r6,r6 mov r2,r4 bl division32R mov r6,r2 @ base <- remainder
lsr r5,#1 @ left shift 1 bit cmp r5,#0 @ end ? bne 2b mov r0,r3
100: @ fin standard de la fonction
pop {r1-r7,lr} @ restaur des registres
bx lr @ retour de la fonction en utilisant lr
/***************************************************/ /* division number 64 bits in 2 registers by number 32 bits */ /***************************************************/ /* r0 contains lower part dividende */ /* r1 contains upper part dividende */ /* r2 contains divisor */ /* r0 return lower part quotient */ /* r1 return upper part quotient */ /* r2 return remainder */ division32R:
push {r3-r9,lr} @ save registers
mov r6,#0 @ init upper upper part remainder !!
mov r7,r1 @ init upper part remainder with upper part dividende
mov r8,r0 @ init lower part remainder with lower part dividende
mov r9,#0 @ upper part quotient
mov r4,#0 @ lower part quotient
mov r5,#32 @ bits number
1: @ begin loop
lsl r6,#1 @ shift upper upper part remainder
lsls r7,#1 @ shift upper part remainder
orrcs r6,#1
lsls r8,#1 @ shift lower part remainder
orrcs r7,#1
lsls r4,#1 @ shift lower part quotient
lsl r9,#1 @ shift upper part quotient
orrcs r9,#1
@ divisor sustract upper part remainder
subs r7,r2
sbcs r6,#0 @ and substract carry
bmi 2f @ négative ?
@ positive or equal
orr r4,#1 @ 1 -> right bit quotient
b 3f
2: @ negative
orr r4,#0 @ 0 -> right bit quotient adds r7,r2 @ and restaur remainder adc r6,#0
3:
subs r5,#1 @ decrement bit size bgt 1b @ end ? mov r0,r4 @ lower part quotient mov r1,r9 @ upper part quotient mov r2,r7 @ remainder
100: @ function end
pop {r3-r9,lr} @ restaur registers
bx lr
/***************************************************/
/* ROUTINES INCLUDE */
/***************************************************/
.include "../affichage.inc"
</lang>
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 Numbers with 2 digits : 1 Numbers with 3 digits : 12 Numbers with 4 digits : 47 Numbers with 5 digits : 192 Numbers with 6 digits : 664
C++
<lang cpp>#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";
}</lang>
- 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
<lang freebasic>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</lang>
- 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
Based on Version 2, takes around 19 seconds. <lang go>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)
}
}</lang>
- 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
Implementation:
<lang J>achilles=: (*/ .>&1 * 1 = +./)@(1{__&q:)"0 strong=: achilles@(5&p:)</lang>
Task examples:
<lang J> 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</lang>
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
<lang>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()
</lang>
- 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
Mathematica/Wolfram Language
<lang Mathematica>ClearAll[PowerfulNumberQ, StrongAchillesNumberQ] PowerfulNumberQ[n_Integer] := AllTrue[FactorInteger[n]All, 2, GreaterEqualThan[2]] AchillesNumberQ[n_Integer] := Module[{divs},
If[PowerfulNumberQ[n], divs = Divisors[n]; If[Length[divs] > 2, divs = divs2 ;; -2; !AnyTrue[Log[#, n] & /@ divs, IntegerQ] , True ] , False ] ]
StrongAchillesNumberQ[n_] := AchillesNumberQ[n] \[And] AchillesNumberQ[EulerPhi[n]]
n = 0; i = 0; Reap[While[n < 50,
i++; If[AchillesNumberQ[i], n++; Sow[i]] ]]2, 1
n = 0; i = 0; Reap[While[n < 20,
i++; If[StrongAchillesNumberQ[i], n++; Sow[i]] ]]2, 1
Tally[IntegerLength /@ Select[Range[9999999], AchillesNumberQ]] // Grid</lang>
- 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}
2 1
3 12
4 47
5 192
6 664
7 2242
Perl
Borrowed, and lightly modified, code from Powerful_numbers
<lang perl>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";
}</lang>
- 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. <lang perl>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)
}</lang> Output is the same.
Phix
You can run this online here, though [slightly outdated and] you should expect a blank screen for about 9s.
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
Timing is going to be system / OS dependent. <lang perl6>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;</lang>
- 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
<lang rust>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());
}</lang>
- 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
Version 1 (Brute force)
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! <lang ecmascript>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
}</lang>
- 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)
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. <lang ecmascript>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)
}</lang>
- 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
<lang XPL0>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);
];
]</lang>
- 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