Achilles numbers: Difference between revisions
→{{header|jq}}
m (→Version 2 (Much faster): Timing correction and other minor changes.) |
|||
| (47 intermediate revisions by 23 users not shown) | |||
Line 1:
{{
{{Wikipedia|Achilles number}}
Line 46:
=={{header|AArch64 Assembly}}==
{{works with|as|Raspberry Pi 3B version Buster 64 bits <br> or android 64 bits with application Termux }}
<syntaxhighlight lang="aarch64 assembly">
/* ARM assembly AARCH64 Raspberry PI 3B */
/* program achilleNumber.s */
/************************************/
/* Constantes */
/************************************/
.include "../includeConstantesARM64.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"
szMessErrGen: .asciz "Error end program.\n"
szMessNbPrem: .asciz "This number is prime !!!.\n"
szMessOverflow: .asciz "Overflow function isPrime.\n"
szMessError: .asciz "\033[31mError !!!\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 16 * NBFACT // factor 8 bytes, number of each factor 8 bytes
/*********************************/
/* code section */
/*********************************/
.text
.global main
main: // entry of program
ldr x0,qAdrszMessTitAchille
bl affichageMess
mov x4,#1 // start number
mov x5,#0 // total counter
mov x6,#0 // line display counter
1:
mov x0,x4
bl controlAchille
cmp x0,#0 // achille number ?
beq 2f // no
mov x0,x4
ldr x1,qAdrsZoneConv
bl conversion10 // call décimal conversion
ldr x0,qAdrszMessNumber
ldr x1,qAdrsZoneConv // insert conversion in message
bl strInsertAtCharInc
bl affichageMess // display message
add x5,x5,#1 // increment counter
add x6,x6,#1 // increment indice line display
cmp x6,#10 // if = 10 new line
bne 2f
mov x6,#0
ldr x0,qAdrszCarriageReturn
bl affichageMess
2:
add x4,x4,#1 // increment number
cmp x5,#MAXI
blt 1b // and loop
ldr x0,qAdrszMessTitStrong
bl affichageMess
mov x4,#1 // start number
mov x5,#0 // total counter
mov x6,#0
3:
mov x0,x4
bl controlAchille
cmp x0,#0
beq 4f
mov x0,x4
bl computeTotient
bl controlAchille
cmp x0,#0
beq 4f
mov x0,x4
ldr x1,qAdrsZoneConv
bl conversion10 // call décimal conversion
ldr x0,qAdrszMessNumber
ldr x1,qAdrsZoneConv // insert conversion in message
bl strInsertAtCharInc
bl affichageMess // display message
add x5,x5,#1
add x6,x6,#1
cmp x6,#10
bne 4f
mov x6,#0
ldr x0,qAdrszCarriageReturn
bl affichageMess
4:
add x4,x4,#1
cmp x5,#MAXI1
blt 3b
ldr x3,icstMaxi2
mov x4,#1 // start number
mov x6,#0 // total counter 2 digits
mov x7,#0 // total counter 3 digits
mov x8,#0 // total counter 4 digits
mov x9,#0 // total counter 5 digits
mov x10,#0 // total counter 6 digits
5:
mov x0,x4
bl controlAchille
cmp x0,#0
beq 10f
mov x0,x4
ldr x1,qAdrsZoneConv
bl conversion10 // call décimal conversion x0 return digit number
cmp x0,#6
bne 6f
add x10,x10,#1
beq 10f
6:
cmp x0,#5
bne 7f
add x9,x9,#1
b 10f
7:
cmp x0,#4
bne 8f
add x8,x8,#1
b 10f
8:
cmp x0,#3
bne 9f
add x7,x7,#1
b 10f
9:
cmp x0,#2
bne 10f
add x6,x6,#1
10:
add x4,x4,#1
cmp x4,x3
blt 5b
mov x0,#2
mov x1,x6
bl displayCounter
mov x0,#3
mov x1,x7
bl displayCounter
mov x0,#4
mov x1,x8
bl displayCounter
mov x0,#5
mov x1,x9
bl displayCounter
mov x0,#6
mov x1,x10
bl displayCounter
b 100f
98:
ldr x0,qAdrszErrorGen
bl affichageMess
100: // standard end of the program
mov x0, #0 // return code
mov x8,EXIT
svc #0 // perform the system call
qAdrszCarriageReturn: .quad szCarriageReturn
qAdrszErrorGen: .quad szErrorGen
qAdrsZoneConv: .quad sZoneConv
qAdrtbZoneDecom: .quad tbZoneDecom
qAdrszMessNumber: .quad szMessNumber
qAdrszMessTitAchille: .quad szMessTitAchille
qAdrszMessTitStrong: .quad szMessTitStrong
icstMaxi2: .quad MAXI2
/******************************************************************/
/* display digit counter */
/******************************************************************/
/* x0 contains limit */
/* x1 contains counter */
displayCounter:
stp x1,lr,[sp,-16]! // save registers
stp x2,x3,[sp,-16]! // save registers
mov x2,x1
ldr x1,qAdrsZoneConv
bl conversion10 // call décimal conversion
ldr x0,qAdrszMessDigitsCounter
ldr x1,qAdrsZoneConv // insert conversion in message
bl strInsertAtCharInc
mov x3,x0
mov x0,x2
ldr x1,qAdrsZoneConv
bl conversion10 // call décimal conversion
mov x0,x3
ldr x1,qAdrsZoneConv // insert conversion in message
bl strInsertAtCharInc
bl affichageMess // display message
100:
ldp x2,x3,[sp],16 // restaur registers
ldp x1,lr,[sp],16 // restaur registers
ret
qAdrszMessDigitsCounter: .quad szMessDigitsCounter
/******************************************************************/
/* control if number is Achille number */
/******************************************************************/
/* x0 contains number */
/* x0 return 0 if not else return 1 */
controlAchille:
stp x1,lr,[sp,-16]! // save registers
stp x2,x3,[sp,-16]! // save registers
stp x4,x5,[sp,-16]! // save registers
mov x4,x0
ldr x1,qAdrtbZoneDecom
bl decompFact // factor decomposition
cmp x0,#-1
beq 99f // error ?
cmp x0,#1 // one only factor or prime ?
ble 98f
mov x1,x0
ldr x0,qAdrtbZoneDecom
mov x2,x4
bl controlDivisor
b 100f
98:
mov x0,#0
b 100f
99:
ldr x0,qAdrszErrorGen
bl affichageMess
100:
ldp x4,x5,[sp],16 // restaur registers
ldp x2,x3,[sp],16 // restaur registers
ldp x1,lr,[sp],16 // restaur registers
ret
/******************************************************************/
/* control divisors function */
/******************************************************************/
/* x0 contains address of divisors area */
/* x1 contains the number of area items */
/* x2 contains number */
controlDivisor:
stp x1,lr,[sp,-16]! // save registers
stp x2,x3,[sp,-16]! // save registers
stp x4,x5,[sp,-16]! // save registers
stp x6,x7,[sp,-16]! // save registers
stp x8,x9,[sp,-16]! // save registers
stp x10,x11,[sp,-16]! // save registers
mov x6,x1 // factors number
mov x8,x2 // save number
mov x9,#0 // indice
mov x4,x0 // save area address
add x5,x4,x9,lsl #4 // compute address first factor
ldr x7,[x5,#8] // load first exposant of factor
add x2,x9,#1
1:
add x5,x4,x2,lsl #4 // compute address next factor
ldr x3,[x5,#8] // load exposant of factor
cmp x3,x7 // factor exposant <> ?
bne 2f // yes -> end verif
add x2,x2,#1 // increment indice
cmp x2,x6 // factor maxi ?
blt 1b // no -> loop
mov x0,#0
b 100f // all exposants are equals
2:
mov x10,x2 // save indice
21:
bge 22f
mov x2,x7 // if x3 < x7 -> inversion
mov x7,x3
mov x3,x2 // x7 is the smaller exposant
22:
mov x0,x3
mov x1,x7 // x7 < x3
bl calPGCDmod
cmp x0,#1
beq 24f // no commun multiple -> ne peux donc pas etre une puissance
23:
add x10,x10,#1 // increment indice
cmp x10,x6 // factor maxi ?
bge 99f // yes -> all exposants are multiples to smaller
add x5,x4,x10,lsl #4
ldr x3,[x5,#8] // load exposant of next factor
cmp x3,x7
beq 23b // for next
b 21b // for compare the 2 exposants
24:
mov x9,#0 // indice
3:
add x5,x4,x9,lsl #4
ldr x7,[x5] // load factor
mul x1,x7,x7 // factor square
udiv x2,x8,x1
msub x3,x1,x2,x8 // compute remainder
cmp x3,#0 // remainder null ?
bne 99f
add x9,x9,#1 // other factor
cmp x9,x6 // factors maxi ?
blt 3b
mov x0,#1 // achille number ok
b 100f
99: // achille not ok
mov x0,0
100:
ldp x10,x11,[sp],16 // restaur registers
ldp x8,x9,[sp],16 // restaur registers
ldp x6,x7,[sp],16 // restaur registers
ldp x4,x5,[sp],16 // restaur registers
ldp x2,x3,[sp],16 // restaur registers
ldp x1,lr,[sp],16 // restaur registers
ret
/***************************************************/
/* Compute pgcd modulo use */
/***************************************************/
/* x0 contains first number */
/* x1 contains second number */
/* x0 return PGCD */
/* if error carry set to 1 */
calPGCDmod:
stp x1,lr,[sp,-16]! // save registres
stp x2,x3,[sp,-16]! // save registres
cbz x0,99f // if = 0 error
cbz x1,99f
cmp x0,0
bgt 1f
neg x0,x0 // if negative inversion number 1
1:
cmp x1,0
bgt 2f
neg x1,x1 // if negative inversion number 2
2:
cmp x0,x1 // compare two numbers
bgt 3f
mov x2,x0 // inversion
mov x0,x1
mov x1,x2
3:
udiv x2,x0,x1 // division
msub x0,x2,x1,x0 // compute remainder
cmp x0,0
bgt 2b // loop
mov x0,x1
cmn x0,0 // clear carry
b 100f
99: // error
mov x0,0
cmp x0,0 // set carry
100:
ldp x2,x3,[sp],16 // restaur des 2 registres
ldp x1,lr,[sp],16 // restaur des 2 registres
ret // retour adresse lr x30
/******************************************************************/
/* compute totient of number */
/******************************************************************/
/* x0 contains number */
computeTotient:
stp x1,lr,[sp,-16]! // save registers
stp x2,x3,[sp,-16]! // save registers
stp x4,x5,[sp,-16]! // save registers
mov x4,x0 // totient
mov x5,x0 // save number
mov x1,#0 // for first divisor
1: // begin loop
mul x3,x1,x1 // compute square
cmp x3,x5 // compare number
bgt 4f // end
add x1,x1,#2 // next divisor
udiv x2,x5,x1
msub x3,x1,x2,x5 // compute remainder
cmp x3,#0 // remainder null ?
bne 3f
2: // begin loop 2
udiv x2,x5,x1
msub x3,x1,x2,x5 // compute remainder
cmp x3,#0
csel x5,x2,x5,eq // new value = quotient
beq 2b
udiv x2,x4,x1 // divide totient
sub x4,x4,x2 // compute new totient
3:
cmp x1,#2 // first divisor ?
mov x0,1
csel x1,x0,x1,eq // divisor = 1
b 1b // and loop
4:
cmp x5,#1 // final value > 1
ble 5f
mov x0,x4 // totient
mov x1,x5 // divide by value
udiv x2,x4,x5 // totient divide by value
sub x4,x4,x2 // compute new totient
5:
mov x0,x4
100:
ldp x4,x5,[sp],16 // restaur registers
ldp x2,x3,[sp],16 // restaur registers
ldp x1,lr,[sp],16 // restaur registers
ret
/******************************************************************/
/* factor decomposition */
/******************************************************************/
/* x0 contains number */
/* x1 contains address of divisors area */
/* x0 return divisors items in table */
decompFact:
stp x1,lr,[sp,-16]! // save registers
stp x2,x3,[sp,-16]! // save registers
stp x4,x5,[sp,-16]! // save registers
stp x6,x7,[sp,-16]! // save registers
stp x8,x9,[sp,-16]! // save registers
mov x5,x1
mov x8,x0 // save number
bl isPrime // prime ?
cmp x0,#1
beq 98f // yes is prime
mov x4,#0 // raz indice
mov x1,#2 // first divisor
mov x6,#0 // previous divisor
mov x7,#0 // number of same divisors
2:
udiv x2,x8,x1 // divide number or other result
msub x3,x2,x1,x8 // compute remainder
cmp x3,#0
bne 5f // if remainder <> zero -> no divisor
mov x8,x2 // else quotient -> new dividende
cmp x1,x6 // same divisor ?
beq 4f // yes
cmp x6,#0 // no but is the first divisor ?
beq 3f // yes
str x6,[x5,x4,lsl #3] // else store in the table
add x4,x4,#1 // and increment counter
str x7,[x5,x4,lsl #3] // store counter
add x4,x4,#1 // next item
mov x7,#0 // and raz counter
3:
mov x6,x1 // new divisor
4:
add x7,x7,#1 // increment counter
b 7f // and loop
/* not divisor -> increment next divisor */
5:
cmp x1,#2 // if divisor = 2 -> add 1
mov x0,#1
mov x3,#2 // else add 2
csel x3,x0,x3,eq
add x1,x1,x3
b 2b
/* divisor -> test if new dividende is prime */
7:
mov x3,x1 // save divisor
cmp x8,#1 // dividende = 1 ? -> end
beq 10f
mov x0,x8 // new dividende is prime ?
mov x1,#0
bl isPrime // the new dividende is prime ?
cmp x0,#1
bne 10f // the new dividende is not prime
cmp x8,x6 // else dividende is same divisor ?
beq 9f // yes
cmp x6,#0 // no but is the first divisor ?
beq 8f // yes it is a first
str x6,[x5,x4,lsl #3] // else store in table
add x4,x4,#1 // and increment counter
str x7,[x5,x4,lsl #3] // and store counter
add x4,x4,#1 // next item
8:
mov x6,x8 // new dividende -> divisor prec
mov x7,#0 // and raz counter
9:
add x7,x7,#1 // increment counter
b 11f
10:
mov x1,x3 // current divisor = new divisor
cmp x1,x8 // current divisor > new dividende ?
ble 2b // no -> loop
/* end decomposition */
11:
str x6,[x5,x4,lsl #3] // store last divisor
add x4,x4,#1
str x7,[x5,x4,lsl #3] // and store last number of same divisors
add x4,x4,#1
lsr x0,x4,#1 // return number of table items
mov x3,#0
str x3,[x5,x4,lsl #3] // store zéro in last table item
add x4,x4,#1
str x3,[x5,x4,lsl #3] // and zero in counter same divisor
b 100f
98:
//ldr x0,qAdrszMessPrime
//bl affichageMess
mov x0,#0 // return code 0 = number is prime
b 100f
99:
ldr x0,qAdrszMessErrGen
bl affichageMess
mov x0,#-1 // error code
b 100f
100:
ldp x8,x9,[sp],16 // restaur registers
ldp x6,x7,[sp],16 // restaur registers
ldp x4,x5,[sp],16 // restaur registers
ldp x2,x3,[sp],16 // restaur registers
ldp x1,lr,[sp],16 // restaur registers
ret
qAdrszMessErrGen: .quad szMessErrGen
/***************************************************/
/* Verification si un nombre est premier */
/***************************************************/
/* x0 contient le nombre à verifier */
/* x0 retourne 1 si premier 0 sinon */
isPrime:
stp x1,lr,[sp,-16]! // save registres
stp x2,x3,[sp,-16]! // save registres
mov x2,x0
sub x1,x0,#1
cmp x2,0
beq 99f // retourne zéro
cmp x2,2 // pour 1 et 2 retourne 1
ble 2f
mov x0,#2
bl moduloPur64
bcs 100f // erreur overflow
cmp x0,#1
bne 99f // Pas premier
cmp x2,3
beq 2f
mov x0,#3
bl moduloPur64
blt 100f // erreur overflow
cmp x0,#1
bne 99f
cmp x2,5
beq 2f
mov x0,#5
bl moduloPur64
bcs 100f // erreur overflow
cmp x0,#1
bne 99f // Pas premier
cmp x2,7
beq 2f
mov x0,#7
bl moduloPur64
bcs 100f // erreur overflow
cmp x0,#1
bne 99f // Pas premier
cmp x2,11
beq 2f
mov x0,#11
bl moduloPur64
bcs 100f // erreur overflow
cmp x0,#1
bne 99f // Pas premier
cmp x2,13
beq 2f
mov x0,#13
bl moduloPur64
bcs 100f // erreur overflow
cmp x0,#1
bne 99f // Pas premier
cmp x2,17
beq 2f
mov x0,#17
bl moduloPur64
bcs 100f // erreur overflow
cmp x0,#1
bne 99f // Pas premier
2:
cmn x0,0 // carry à zero pas d'erreur
mov x0,1 // premier
b 100f
99:
cmn x0,0 // carry à zero pas d'erreur
mov x0,#0 // Pas premier
100:
ldp x2,x3,[sp],16 // restaur des 2 registres
ldp x1,lr,[sp],16 // restaur des 2 registres
ret // retour adresse lr x30
/**************************************************************/
/********************************************************/
/* Calcul modulo de b puissance e modulo m */
/* Exemple 4 puissance 13 modulo 497 = 445 */
/********************************************************/
/* x0 nombre */
/* x1 exposant */
/* x2 modulo */
moduloPur64:
stp x1,lr,[sp,-16]! // save registres
stp x3,x4,[sp,-16]! // save registres
stp x5,x6,[sp,-16]! // save registres
stp x7,x8,[sp,-16]! // save registres
stp x9,x10,[sp,-16]! // save registres
cbz x0,100f
cbz x1,100f
mov x8,x0
mov x7,x1
mov x6,1 // resultat
udiv x4,x8,x2
msub x9,x4,x2,x8 // contient le reste
1:
tst x7,1
beq 2f
mul x4,x9,x6
umulh x5,x9,x6
//cbnz x5,99f
mov x6,x4
mov x0,x6
mov x1,x5
bl divisionReg128U
cbnz x1,99f // overflow
mov x6,x3
2:
mul x8,x9,x9
umulh x5,x9,x9
mov x0,x8
mov x1,x5
bl divisionReg128U
cbnz x1,99f // overflow
mov x9,x3
lsr x7,x7,1
cbnz x7,1b
mov x0,x6 // result
cmn x0,0 // carry à zero pas d'erreur
b 100f
99:
ldr x0,qAdrszMessOverflow
bl affichageMess
cmp x0,0 // carry à un car erreur
mov x0,-1 // code erreur
100:
ldp x9,x10,[sp],16 // restaur des 2 registres
ldp x7,x8,[sp],16 // restaur des 2 registres
ldp x5,x6,[sp],16 // restaur des 2 registres
ldp x3,x4,[sp],16 // restaur des 2 registres
ldp x1,lr,[sp],16 // restaur des 2 registres
ret // retour adresse lr x30
qAdrszMessOverflow: .quad szMessOverflow
/***************************************************/
/* division d un nombre de 128 bits par un nombre de 64 bits */
/***************************************************/
/* x0 contient partie basse dividende */
/* x1 contient partie haute dividente */
/* x2 contient le diviseur */
/* x0 retourne partie basse quotient */
/* x1 retourne partie haute quotient */
/* x3 retourne le reste */
divisionReg128U:
stp x6,lr,[sp,-16]! // save registres
stp x4,x5,[sp,-16]! // save registres
mov x5,#0 // raz du reste R
mov x3,#128 // compteur de boucle
mov x4,#0 // dernier bit
1:
lsl x5,x5,#1 // on decale le reste de 1
tst x1,1<<63 // test du bit le plus à gauche
lsl x1,x1,#1 // on decale la partie haute du quotient de 1
beq 2f
orr x5,x5,#1 // et on le pousse dans le reste R
2:
tst x0,1<<63
lsl x0,x0,#1 // puis on decale la partie basse
beq 3f
orr x1,x1,#1 // et on pousse le bit de gauche dans la partie haute
3:
orr x0,x0,x4 // position du dernier bit du quotient
mov x4,#0 // raz du bit
cmp x5,x2
blt 4f
sub x5,x5,x2 // on enleve le diviseur du reste
mov x4,#1 // dernier bit à 1
4:
// et boucle
subs x3,x3,#1
bgt 1b
lsl x1,x1,#1 // on decale le quotient de 1
tst x0,1<<63
lsl x0,x0,#1 // puis on decale la partie basse
beq 5f
orr x1,x1,#1
5:
orr x0,x0,x4 // position du dernier bit du quotient
mov x3,x5
100:
ldp x4,x5,[sp],16 // restaur des 2 registres
ldp x6,lr,[sp],16 // restaur des 2 registres
ret // retour adresse lr x30
/***************************************************/
/* ROUTINES INCLUDE */
/***************************************************/
.include "../includeARM64.inc"
</syntaxhighlight>
<pre>
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
</pre>
=={{header|ALGOL 68}}==
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
{{libheader|ALGOL 68-primes}}
<syntaxhighlight 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</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|ARM Assembly}}==
{{works with|as|Raspberry Pi <br> or android 32 bits with application Termux}}
<syntaxhighlight 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"
</syntaxhighlight>
<pre>
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
</pre>
=={{header|BASIC256}}==
{{trans|FreeBASIC}}
<syntaxhighlight lang="vb">function GCD(n, d)
if(d = 0) then return n else return GCD(d, n % d)
end function
function Totient(n)
tot = 0
for m = 1 to n
if GCD(m, n) = 1 then tot += 1
next m
return tot
end function
function isPowerful(m)
n = m
f = 2
l = sqr(m)
if m <= 1 then return false
while true
q = n/f
if (n % f) = 0 then
if (m %(f*f)) then return false
n = q
if f > n then exit while
else
f += 1
if f > l then
if (m % (n*n)) then return false
exit while
end if
end if
end while
return true
end function
function isAchilles(n)
if not isPowerful(n) then return false
m = 2
a = m*m
do
do
if a = n then return false
a *= m
until a > n
m += 1
a = m*m
until a > n
return true
end function
print "First 50 Achilles numbers:"
num = 0
n = 1
do
if isAchilles(n) then
print rjust(n, 5);
num += 1
if (num % 10) <> 0 then print " "; else print
end if
n += 1
until num >= 50
print : print
print "First 20 strong Achilles numbers:"
num = 0
n = 1
do
if isAchilles(n) and isAchilles(Totient(n)) then
print rjust(n, 5);
num += 1
if (num % 10) <> 0 then print " "; else print
end if
n += 1
until num >= 20
print : print
print "Number 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</syntaxhighlight>
[[https://dotnetfiddle.net/TCD2WC You may run it online!]]
=={{header|C++}}==
{{trans|Wren}}
{{libheader|Boost}}
<syntaxhighlight 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";
}</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Delphi}}==
{{works with|Delphi|6.0}}
{{libheader|SysUtils,StdCtrls}}
<syntaxhighlight lang="Delphi">
function GetTotient(N: integer): integer;
{Calculate Euler's Totient}
var M: integer;
begin
Result:= 0;
for M:= 1 to N do
if GreatestCommonDivisor(M, N) = 1 then
Result:= Result+1;
end;
function IsPowerfulNum(N: integer): boolean;
{Is a powerful number i.e. all prime factors square are divisor}
var I: integer;
var IA: TIntegerDynArray;
begin
Result:=False;
GetPrimeFactors(N,IA);
for I:=0 to High(IA) do
if (N mod (IA[I]*IA[I]))<>0 then exit;
Result:=True;
end;
function CanBeMtoK(N: integer): boolean;
{Can N be represented as M^K?}
var M, A: integer;
begin
Result:=False;
M:= 2;
A:= M*M;
repeat
begin
while true do
begin
if A = N then exit;
if A > N then break;
A:= A*M;
end;
M:= M+1;
A:= M*M;
end
until A > N;
Result:=True;
end;
function IsAchilles(N: integer): boolean;
{Achilles = Is Powerful and can be M^K}
begin
Result:=IsPowerfulNum(N) and CanBeMtoK(N);
end;
procedure AchillesNumbers(Memo: TMemo);
var I,Cnt,Digits: integer;
var S: string;
var DigCnt: array [0..5] of integer;
begin
Memo.Lines.Add('First 50 Achilles numbers:');
Cnt:=0; S:='';
for I:=2 to high(Integer) do
if IsAchilles(I) then
begin
Inc(Cnt);
S:=S+Format('%6d',[I]);
if (Cnt mod 10)=0 then S:=S+CRLF;
if Cnt>=50 then break;
end;
Memo.Lines.Add(S);
Memo.Lines.Add('First 20 Strong Achilles Numbers:');
Cnt:=0; S:='';
for I:=2 to high(Integer) do
if IsAchilles(I) then
if IsAchilles(GetTotient(I)) then
begin
Inc(Cnt);
S:=S+Format('%6d',[I]);
if (Cnt mod 10)=0 then S:=S+CRLF;
if Cnt>=20 then break;
end;
Memo.Lines.Add(S);
Memo.Lines.Add('Digits Counts:');
for I:=0 to High(DigCnt) do DigCnt[I]:=0;
for I:=2 to high(Integer) do
begin
Digits:=NumberOfDigits(I);
if Digits>High(DigCnt) then break;
if IsAchilles(I) then Inc(DigCnt[Digits]);
end;
Memo.Lines.Add('Last Count: '+IntToStr(I));
for I:=0 to High(DigCnt) do
if DigCnt[I]>0 then
begin
Memo.Lines.Add(Format('%d digits: %d',[I,DigCnt[I]]));
end
end;
</syntaxhighlight>
{{out}}
<pre>
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
Digits Counts:
Last Count: 100000
2 digits: 1
3 digits: 12
4 digits: 47
5 digits: 192
</pre>
=={{header|EasyLang}}==
{{trans|BASIC256}}
<syntaxhighlight>
func gcd n d .
if d = 0
return n
.
return gcd d (n mod d)
.
func totient n .
for m = 1 to n
if gcd m n = 1
tot += 1
.
.
return tot
.
func isPowerful m .
n = m
f = 2
l = sqrt m
if m <= 1
return 0
.
while 1 = 1
q = n div f
if n mod f = 0
if m mod (f * f) <> 0
return 0
.
n = q
if f > n
return 1
.
else
f += 1
if f > l
if m mod (n * n) <> 0
return 0
.
return 1
.
.
.
.
func isAchilles n .
if isPowerful n = 0
return 0
.
m = 2
a = m * m
repeat
repeat
if a = n
return 0
.
a *= m
until a > n
.
m += 1
a = m * m
until a > n
.
return 1
.
print "First 50 Achilles numbers:"
n = 1
repeat
if isAchilles n = 1
write n & " "
num += 1
.
n += 1
until num >= 50
.
print ""
print ""
print "First 20 strong Achilles numbers:"
num = 0
n = 1
repeat
if isAchilles n = 1 and isAchilles totient n = 1
write n & " "
num += 1
.
n += 1
until num >= 20
.
print ""
print ""
print "Number of Achilles numbers with 2 to 5 digits:"
a = 10
b = 100
for i = 2 to 5
num = 0
for n = a to b - 1
if isAchilles n = 1
num += 1
.
.
write num & " "
a = b
b *= 10
.
</syntaxhighlight>
{{out}}
<pre>
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 to 5 digits:
1 12 47 192
</pre>
=={{header|Factor}}==
{{works with|Factor|0.99 2022-04-03}}
<syntaxhighlight lang="factor">USING: assocs combinators.short-circuit formatting grouping io
kernel lists lists.lazy math math.functions math.primes.factors
prettyprint ranges sequences ;
: achilles? ( n -- ? )
group-factors values {
[ [ 1 > ] all? ]
[ unclip-slice [ simple-gcd ] reduce 1 = ]
} 1&& ;
: achilles ( -- list )
2 lfrom [ achilles? ] lfilter ;
: strong-achilles ( -- list )
achilles [ totient achilles? ] lfilter ;
: show ( n list -- ) ltake list>array 10 group simple-table. ;
: <order-of-magnitude> ( n -- range )
1 - 10^ dup 10 * [a..b) ;
"First 50 Achilles numbers:" print
50 achilles show nl
"First 30 strong Achilles numbers:" print
30 strong-achilles show nl
"Number of Achilles numbers with" print
{ 2 3 4 5 } [
dup <order-of-magnitude> [ achilles? ] count
"%d digits: %d\n" printf
] each</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|FreeBASIC}}==
<syntaxhighlight 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</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|Go}}==
{{trans|Wren}}
Based on Version 2, takes around 19 seconds.
<syntaxhighlight 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)
}
}</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|J}}==
Implementation:
<
strong=: achilles@(5&p:)</
Task examples:
<
72 108 200 288 392 432 500 648 675 800
864 968 972 1125 1152 1323 1352 1372 1568 1800
Line 75 ⟶ 2,377:
192
+/achilles (+i.)/1 9*10^<:6
664</
Explanation of the code:
Line 85 ⟶ 2,387:
<tt>*/ .>&1</tt> is true if all the values in a list are greater than 1 (0 if not).
<tt>"0</tt> 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 <tt>*/ .>&1</tt> to <tt>(0 = 1 e. ])</tt> 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...)
<tt>5&p:</tt> is euler's totient function.
<tt>(#~ predicate) list</tt> selects the elements of <tt>list</tt> where <tt>predicate</tt> is true.
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
public class AchillesNumbers {
private Map<Integer, Boolean> pps = new HashMap<>();
public int totient(int n) {
int tot = n;
int i = 2;
while (i * i <= n) {
if (n % i == 0) {
while (n % i == 0) {
n /= i;
}
tot -= tot / i;
}
if (i == 2) {
i = 1;
}
i += 2;
}
if (n > 1) {
tot -= tot / n;
}
return tot;
}
public void getPerfectPowers(int maxExp) {
double upper = Math.pow(10, maxExp);
for (int i = 2; i <= Math.sqrt(upper); i++) {
double fi = i;
double p = fi;
while (true) {
p *= fi;
if (p >= upper) {
break;
}
pps.put((int) p, true);
}
}
}
public Map<Integer, Boolean> getAchilles(int minExp, int maxExp) {
double lower = Math.pow(10, minExp);
double upper = Math.pow(10, maxExp);
Map<Integer, Boolean> achilles = new HashMap<>();
for (int b = 1; b <= (int) Math.cbrt(upper); b++) {
int b3 = b * b * b;
for (int a = 1; a <= (int) Math.sqrt(upper); a++) {
int p = b3 * a * a;
if (p >= (int) upper) {
break;
}
if (p >= (int) lower) {
if (!pps.containsKey(p)) {
achilles.put(p, true);
}
}
}
}
return achilles;
}
public static void main(String[] args) {
AchillesNumbers an = new AchillesNumbers();
int maxDigits = 8;
an.getPerfectPowers(maxDigits);
Map<Integer, Boolean> achillesSet = an.getAchilles(1, 5);
List<Integer> achilles = new ArrayList<>(achillesSet.keySet());
Collections.sort(achilles);
System.out.println("First 50 Achilles numbers:");
for (int i = 0; i < 50; i++) {
System.out.printf("%4d ", achilles.get(i));
if ((i + 1) % 10 == 0) {
System.out.println();
}
}
System.out.println("\nFirst 30 strong Achilles numbers:");
List<Integer> strongAchilles = new ArrayList<>();
int count = 0;
for (int n = 0; count < 30; n++) {
int tot = an.totient(achilles.get(n));
if (achillesSet.containsKey(tot)) {
strongAchilles.add(achilles.get(n));
count++;
}
}
for (int i = 0; i < 30; i++) {
System.out.printf("%5d ", strongAchilles.get(i));
if ((i + 1) % 10 == 0) {
System.out.println();
}
}
System.out.println("\nNumber of Achilles numbers with:");
for (int d = 2; d <= maxDigits; d++) {
int ac = an.getAchilles(d - 1, d).size();
System.out.printf("%2d digits: %d\n", d, ac);
}
}
}
</syntaxhighlight>
{{ out }}
<pre>
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
</pre>
=={{header|jq}}==
'''Adapted from the "fast" [[#Wren|Wren]] entry'''
'''Works with jq, the C implementation of jq'''
'''Works with gojq, the Go implementation of jq'''
`nwise` is included below because, as of this writing, gojq does not include `_nwise`.
<syntaxhighlight lang="jq">
# Require $n > 0
def nwise($n):
def _n: if length <= $n then . else .[:$n] , (.[$n:] | _n) end;
if $n <= 0 then "nwise: argument should be non-negative" else _n end;
### Part 1 - generic functions
# Ensure $x is in the input sorted array
def ensure($x):
bsearch($x) as $i
| if $i >= 0 then .
else (-1-$i) as $i
| .[:$i] + [$x] + .[$i:]
end ;
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
def table($wide; $pad):
nwise($wide) | map(lpad($pad)) | join(" ");
def count(s): reduce s as $x (0; .+1);
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);
# jq optimizes the recursive call of _gcd in the following:
def gcd($a;$b):
def _gcd:
if .[1] != 0 then [.[1], .[0] % .[1]] | _gcd else .[0] end;
[$a,$b] | _gcd ;
def totient:
. as $n
| count( range(0; .) | select( gcd($n; .) == 1) );
# Emit a sorted array
def getPerfectPowers( $maxExp ):
(10 | power($maxExp)) as $upper
| reduce range( 2; 1 + ($upper|sqrt|floor)) as $i ({pps: []};
.p = $i * $i
| until (.p >= $upper;
.pps += [ .p ]
| .p *= $i) )
| .pps
| sort;
# Input: a sufficiently long array of perfect powers in order
def getAchilles($minExp; $maxExp):
def cbrt: pow(.; 1/3);
. as $pps
| (10 | power($minExp)) as $lower
| (10 | power($maxExp)) as $upper
| ($upper | sqrt | floor) as $sqrtupper
| reduce range(1; 1 + ($upper|cbrt|floor)) as $b ({achilles: []};
($b | .*.*.) as $b3
| .done = false
| .a = 1
| until(.done or (.a > $sqrtupper);
($b3 * .a * .a) as $p
| if $p >= $upper then .done = true
elif $p >= $lower and ($pps | bsearch($p) < 0)
then .achilles |= ensure($p)
end
| .a += 1 ) )
| .achilles;
def task($maxDigits):
getPerfectPowers($maxDigits)
| . as $perfectPowers
| getAchilles(1; 5)
| . as $achilles
| "First 50 Achilles numbers:",
(.[:50] | table(10;5)),
"\nFirst 30 strong Achilles numbers:",
({ strongAchilles:[], count:0, n:0 }
| until (.count >= 30;
$achilles[.n] as $a
| ($a | totient) as $tot
| if ($achilles | bsearch($tot)) >= 0
then .strongAchilles |= ensure($a)
| .count += 1
end
| .n += 1 )
| (.strongAchilles | table(10;5) ),
"\nNumber of Achilles numbers with:",
( range(2; 1+$maxDigits) as $d
| ($perfectPowers|getAchilles($d-1; $d)|length) as $ac
| "\($d) digits: \($ac)" ) )
;
task(10)
</syntaxhighlight>
{{output}}
<pre>
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
</pre>
=={{header|Julia}}==
<syntaxhighlight lang="text">using Primes
isAchilles(n) = (p = [x[2] for x in factor(n).pe]; all(>(1), p) && gcd(p) == 1)
Line 128 ⟶ 2,692:
teststrongachilles()
</
<pre>
First 50 Achilles numbers:
Line 156 ⟶ 2,720:
100000:999999 664
</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight 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 = divs[[2 ;; -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</syntaxhighlight>
{{out}}
<pre>{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</pre>
=={{header|Nim}}==
{{trans|Wren}}
<syntaxhighlight lang="Nim">import std/[algorithm, sets, math, sequtils, strformat, strutils]
const MaxDigits = 15
func getPerfectPowers(maxExp: int): HashSet[int] =
let upper = 10^maxExp
for i in 2..int(sqrt(upper.toFloat)):
var p = i
while p < upper div i:
p *= i
result.incl p
let pps = getPerfectPowers(MaxDigits)
proc getAchilles(minExp, maxExp: int): HashSet[int] =
let lower = 10^minExp
let upper = 10^maxExp
for b in 1..int(cbrt(upper.toFloat)):
let b3 = b * b * b
for a in 1..int(sqrt(upper.toFloat)):
let p = b3 * a * a
if p >= upper: break
if p >= lower:
if p notin pps: result.incl p
### Part 1 ###
let achillesSet = getAchilles(1, 6)
let achilles = sorted(achillesSet.toSeq)
echo "First 50 Achilles numbers:"
for i in 0..49:
let n = achilles[i]
stdout.write &"{n:>4}"
stdout.write if i mod 10 == 9: '\n' else: ' '
### Part 2 ###
func totient(n: int): int =
var n = n
result = n
var i = 2
while i * i <= n:
if n mod i == 0:
while n mod i == 0:
n = int(n / i)
result -= int(result / i)
if i == 2: i = 1
inc i, 2
if n > 1:
result -= int(result / n)
echo "\nFirst 50 strong Achilles numbers:"
var strongAchilles: seq[int]
var count = 0
for n in achilles:
let tot = totient(n)
if tot in achillesSet:
strongAchilles.add n
inc count
if count == 50: break
for i, n in strongAchilles:
stdout.write &"{n:>6}"
stdout.write if i mod 10 == 9: '\n' else: ' '
### Part 3 ###
echo "\nNumber of Achilles numbers with:"
for d in 2..MaxDigits:
let ac = getAchilles(d - 1, d).len
echo &"{d:>2} digits: {ac}"
</syntaxhighlight>
{{out}}
<pre>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
</pre>
=={{header|Perl}}==
Borrowed, and lightly modified, code from [[Powerful_numbers]]
{{libheader|ntheory}}
<syntaxhighlight lang="perl">use strict;
use warnings;
use feature <say current_sub>;
use experimental 'signatures';
use List::AllUtils <max head>;
use ntheory <is_square_free euler_phi>;
use Math::AnyNum <:overload idiv is_power 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 (@achilles, %Ahash, @strong);
my @P = powerful_numbers(10**9, 2);
!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";
}</syntaxhighlight>
{{out}}
<pre>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</pre>
Here is a translation from Wren version 2, as an alternative.
<syntaxhighlight 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)
}</syntaxhighlight>
Output is the same.
=={{header|Phix}}==
{{libheader|Phix/online}}
You can run this online [http://phix.x10.mx/p2js/achilles.htm here], though [slightly outdated and] you should expect a blank screen for about 9s.
{{trans|Wren}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #
<span style="color: #
<span style="color: #004080;">atom</span> <span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()</span>
<span style="color: #
<span style="color:
<span style="color: #008080;">procedure</span> <span style="color: #000000;">getPerfectPowers</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">maxExp</span><span style="color: #0000FF;">)</span>
<span style="color: #
<span style="color: #008080;">if</span> <span style="color: #000000;">
<span style="color: #7060A8;">setd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #004600;">true</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pps</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #
<span style="color: #004080;">integer</span> <span style="color: #000000;">bmax</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">hi10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">amax</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">hi10</span><span style="color: #0000FF;">))</span>
<span style="color: #
<span style="color: #
<span style="color: #008080;">
<span style="color: #004080;">integer</span> <span style="color: #000000;">
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">
<span style="color: #000000;">achilles</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">unique</span><span style="color: #0000FF;">(</span><span style="color: #000000;">achilles</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">achilles</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #000000;">getPerfectPowers</span><span style="color: #0000FF;">(</span><span style="color: #000000;">maxDigits</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">achilles</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">get_achilles</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">strong_achilles</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">totient</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_eq</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">apply</span><span style="color: #0000FF;">(</span><span style="color: #004600;">true</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">gcd</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span><span style="color: #000000;">n</span><span style="color: #0000FF;">}),</span><span style="color: #000000;">1</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">find</span><span style="color: #0000FF;">(</span><span style="color: #000000;">totient</span><span style="color: #0000FF;">,</span><span style="color: #000000;">achilles</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">join_by</span><span style="color: #0000FF;">(</span><span style="color: #000000;">achilles</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">50</span><span style="color: #0000FF;">],</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" "</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">:=</span><span style="color: #008000;">"%4d"</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">sa</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">filter</span><span style="color: #0000FF;">(</span><span style="color: #000000;">achilles</span><span style="color: #0000FF;">,</span><span style="color: #000000;">strong_achilles</span><span style="color: #0000FF;">)[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">30</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">ssa</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">join_by</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sa</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" "</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">:=</span><span style="color: #008000;">"%5d"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"First 50 Achilles numbers:\n%s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"First 30 strong Achilles numbers:\n%s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">ssa</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">maxDigits</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Achilles numbers with %d digits:%d\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">d</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">get_achilles</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">))})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #0000FF;">?</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)</span>
<!--</
{{out}}
<pre>
Line 268 ⟶ 3,096:
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"
</pre>
=={{header|Python}}==
<syntaxhighlight lang="python">from math import gcd
from sympy import factorint
def is_Achilles(n):
p = factorint(n).values()
return all(i > 1 for i in p) and gcd(*p) == 1
def is_strong_Achilles(n):
return is_Achilles(n) and is_Achilles(totient(n))
def test_strong_Achilles(nachilles, nstrongachilles):
# task 1
print('First', nachilles, 'Achilles numbers:')
n, found = 0, 0
while found < nachilles:
if is_Achilles(n):
found += 1
print(f'{n: 8,}', end='\n' if found % 10 == 0 else '')
n += 1
# task 2
print('\nFirst', nstrongachilles, 'strong Achilles numbers:')
n, found = 0, 0
while found < nstrongachilles:
if is_strong_Achilles(n):
found += 1
print(f'{n: 9,}', end='\n' if found % 10 == 0 else '')
n += 1
# task 3
print('\nCount of Achilles numbers for various intervals:')
intervals = [[10, 99], [100, 999], [1000, 9999], [10000, 99999], [100000, 999999]]
for interval in intervals:
print(f'{interval}:', sum(is_Achilles(i) for i in range(*interval)))
test_strong_Achilles(50, 100)
</syntaxhighlight>{{out}}
<pre>
First 50 Achilles numbers:
72 108 200 288 392 432 500 648 675 800
864 968 972 1,125 1,152 1,323 1,352 1,372 1,568 1,800
1,944 2,000 2,312 2,592 2,700 2,888 3,087 3,200 3,267 3,456
3,528 3,872 3,888 4,000 4,232 4,500 4,563 4,608 5,000 5,292
5,324 5,400 5,408 5,488 6,075 6,125 6,272 6,728 6,912 7,200
First 100 strong Achilles numbers:
500 864 1,944 2,000 2,592 3,456 5,000 10,125 10,368 12,348
12,500 16,875 19,652 19,773 30,375 31,104 32,000 33,275 37,044 40,500
49,392 50,000 52,488 55,296 61,731 64,827 67,500 69,984 78,608 80,000
81,000 83,349 84,375 93,312 108,000 111,132 124,416 128,000 135,000 148,176
151,875 158,184 162,000 165,888 172,872 177,957 197,568 200,000 202,612 209,952
219,488 221,184 237,276 243,000 246,924 253,125 266,200 270,000 273,375 296,352
320,000 324,000 333,396 364,500 397,953 405,000 432,000 444,528 453,789 455,877
493,848 497,664 500,000 518,616 533,871 540,000 555,579 583,443 605,052 607,500
629,856 632,736 648,000 663,552 665,500 666,792 675,000 691,488 740,772 750,141
790,272 800,000 810,448 820,125 831,875 877,952 949,104 972,000 987,696 1,000,188
Count of Achilles numbers for various intervals:
[10, 99]: 1
[100, 999]: 12
[1000, 9999]: 47
[10000, 99999]: 192
[100000, 999999]: 664
</pre>
=={{header|Quackery}}==
<code>gcd</code> is defined at [[Greatest common divisor#Quackery]].
<code>primefactors</code> is defined at [[Prime decomposition#Quackery]].
<code>totient</code> is defined at [[Totient function#Quackery]].
<syntaxhighlight lang="Quackery"> [ ' [ 1 ] swap
behead swap witheach
[ dup dip
[ = iff
[ -1 pluck
1+ join ]
else
[ 1 join ] ] ]
drop ] is runs ( [ --> [ )
[ 1 over find swap found not ] is powerful ( [ --> b )
[ behead swap witheach gcd
1 = ] is imperfect ( [ --> b )
[ dup 2 < iff
[ drop false ] done
primefactors runs
dup powerful iff
imperfect
else [ drop false ] ] is achilles ( [ --> b )
[ dup achilles iff
[ totient achilles ]
else [ drop false ] ] is strong ( [ --> b )
[] 0
[ 1+ dup achilles if
[ tuck join swap ]
over size 50 = until ]
drop
say "First fifty achilles numbers:" cr
echo
cr cr
[] 0
[ 1+ dup strong if
[ tuck join swap ]
over size 20 = until ]
drop
say "First twenty strong achilles numbers:" cr
echo
cr cr
0 100 times
[ i^ achilles if 1+ ]
say "Achilles numbers with 2 digits: " echo
cr
0 900 times
[ i^ 100 + achilles if 1+ ]
say "Achilles numbers with 3 digits: " echo
cr
0 9000 times
[ i^ 1000 + achilles if 1+ ]
say "Achilles numbers with 4 digits: " echo
cr
0 90000 times
[ i^ 10000 + achilles if 1+ ]
say "Achilles numbers with 5 digits: " echo
cr
</syntaxhighlight>
{{out}}
<pre>First fifty 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 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 numbers with 2 digits: 1
Achilles numbers with 3 digits: 12
Achilles numbers with 4 digits: 47
Achilles numbers with 5 digits: 192
</pre>
=={{header|Raku}}==
Timing is going to be system / OS dependent.
<syntaxhighlight lang="raku"
use Math::Root;
Line 319 ⟶ 3,299:
say "$_ digits: " ~ +$achilles{$_} // 0 for 2..9;
printf "\n%.1f total elapsed seconds\n", now - INIT now;</
{{out}}
<pre>First 50 Achilles numbers:
Line 351 ⟶ 3,331:
410.4 total elapsed seconds
</pre>
=={{header|RPL}}==
Based on Wikipedia definition: n = p<sub>1</sub><sup>a<sub>1</sub></sup>p<sub>2</sub><sup>a<sub>2</sub></sup>…p<sub>k</sub><sup>a<sub>k</sub></sup> is an Achilles number if min(a<sub>1</sub>, a<sub>2</sub>, …, a<sub>k</sub>) ≥ 2 and gcd(a<sub>1</sub>, a<sub>2</sub>, …, a<sub>k</sub>) = 1.
{{works with|HP|49g}}
≪ FACTORS
'''IF''' DUP SIZE 4 < '''THEN''' DROP 0
'''ELSE'''
{ } 2 3 PICK SIZE '''FOR''' j <span style="color:grey">@ keeps exponents only</span>
OVER j GET +
2 '''STEP''' NIP
→ a
≪ a ≪ MIN ≫ STREAM 2 ≥ a ≪ GCD ≫ STREAM 1 == AND
≫
'''END'''
≫ '<span style="color:blue">ACH?</span>' STO
≪ { } 1 → n
≪ '''WHILE''' DUP SIZE 50 < '''REPEAT'''
'''IF''' 'n' INCR <span style="color:blue">ACH?</span> '''THEN''' n + '''END'''
'''END'''
≫ ≫ EVAL
≪ { } 1 → n
≪ '''WHILE''' DUP SIZE 20 < '''REPEAT'''
'''IF''' 'n' INCR <span style="color:blue">ACH?</span>
'''THEN''' IF n EULER <span style="color:blue">ACH?</span> '''THEN''' n + '''END END'''
'''END'''
≫ ≫ EVAL
{{out}}
<pre>
2: {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}
1: {500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500}
</pre>
First 50 Achilles numbers found in 53 seconds on the ''iHP48'' emulator running on a iPhone Xr; first 20 strong Achilles numbers found in 5 minutes 13 seconds on the same platform.
=={{header|Ruby}}==
<syntaxhighlight lang="ruby">require 'prime'
def achilles?(n)
exponents = n.prime_division.map(&:last)
exponents.none?(1) && exponents.inject(&:gcd) == 1
end
def 𝜑(n)
n.prime_division.inject(1){|res, (pr, exp)| res * (pr-1) * pr**(exp-1) }
end
achilleses = (2..).lazy.select{|n| achilles?(n) }
n = 50
puts "First #{n} Achilles numbers:"
achilleses.first(n).each_slice(10){|s| puts "%9d"*s.size % s}
puts "\nFirst #{n} strong Achilles numbers:"
achilleses.select{|ach| achilles?(𝜑(ach)) }.first(n).each_slice(10){|s| puts "%9d"*s.size % s }
puts
counts = achilleses.take_while{|ach| ach < 1000000}.map{|a| a.digits.size }.tally
counts.each{|k, v| puts "#{k} digits: #{v}" }
</syntaxhighlight>
{{out}}
<pre>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
</pre>
=={{header|Rust}}==
{{trans|Wren}}
<syntaxhighlight 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());
}</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Sidef}}==
<syntaxhighlight lang="ruby">var P = 2.powerful(1e9)
var achilles = Set(P.grep{ !.is_power }...)
var strong_achilles = achilles.grep { achilles.has(.phi) }
say "First 50 Achilles numbers:"
achilles.sort.first(50).slices(10).each { .map{'%4s'%_}.join(' ').say }
say "\nFirst 30 strong Achilles numbers:"
strong_achilles.sort.first(30).slices(10).each { .map{'%5s'%_}.join(' ').say }
say "\nNumber of Achilles numbers with:"
achilles.to_a.group_by{.len}.sort_by{|k| k.to_i }.each_2d{|a,b|
say "#{a} digits: #{b.len}"
}</syntaxhighlight>
{{out}}
<pre>
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
</pre>
Line 359 ⟶ 3,603:
===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!
<
import "./seq" for Lst
import "./fmt" for Fmt
Line 366 ⟶ 3,610:
var limit = 10.pow(maxDigits)
var c = Int.primeSieve(limit-1, false)
var isPerfectPower = Fn.new { |n|
Line 424 ⟶ 3,653:
var isStrongAchilles = Fn.new { |n|
if (!isAchilles.call(n)) return false
var tot =
return isAchilles.call(tot)
}
Line 439 ⟶ 3,668:
n = n + 1
}
Fmt.tprint("$4d", achilles, 10)
System.print("\nFirst 30 strong Achilles numbers:")
Line 452 ⟶ 3,681:
n = n + 1
}
Fmt.tprint("$5d", strongAchilles, 10)
System.print("\nNumber of Achilles numbers with:")
Line 463 ⟶ 3,692:
System.print("%(i) digits: %(count)")
pow = pow * 10
}</
{{out}}
Line 490 ⟶ 3,719:
===Version 2 (Much faster)===
{{libheader|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
Ridiculously fast compared to the previous method:
<
import "./seq" for Lst
import "./math" for Int
import "./fmt" for Fmt
var
var getPerfectPowers = Fn.new { |maxExp|
var
while ((p = p
}
}
var
var
var achilles = Set.new() // avoids duplicates
for (b in 1..upper.cbrt.floor) {
var b3 = b * b * b
Line 533 ⟶ 3,747:
if (p >= upper) break
if (p >= lower) {
if (!
}
}
Line 539 ⟶ 3,753:
return achilles
}
var maxDigits = 15
getPerfectPowers.call(maxDigits)
var achillesSet = getAchilles.call(1, 5) // enough for first 2 parts
Line 545 ⟶ 3,762:
System.print("First 50 Achilles numbers:")
Fmt.tprint("$4d", achilles.take(50), 10)
System.print("\nFirst 30 strong Achilles numbers:")
Line 552 ⟶ 3,769:
var n = 0
while (count < 30) {
var tot =
if (achillesSet.contains(tot)) {
strongAchilles.add(achilles[n])
Line 559 ⟶ 3,776:
n = n + 1
}
Fmt.tprint("$5d", strongAchilles, 10)
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)
}</
{{out}}
Line 593 ⟶ 3,809:
10 digits: 77330
11 digits: 247449
12 digits: 788855
13 digits: 2508051
14 digits: 7960336
15 digits: 25235383
</pre>
=={{header|XPL0}}==
<syntaxhighlight 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);
];
]</syntaxhighlight>
{{out}}
<pre>
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
</pre>
| |||