Achilles numbers

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


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


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

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

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


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


E.G.

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

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

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


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


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


See also


AArch64 Assembly

Works with: as version Raspberry Pi 3B version Buster 64 bits
or android 64 bits with application Termux
/* 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"
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

ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32
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
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

Works with: as version Raspberry Pi
or android 32 bits with application Termux
/* 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"
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

BASIC256

Translation of: FreeBASIC
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


[You may run it online!]

C++

Translation of: Wren
Library: Boost
#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

using boost::multiprecision::uint128_t;

template <typename T> void unique_sort(std::vector<T>& vector) {
    std::sort(vector.begin(), vector.end());
    vector.erase(std::unique(vector.begin(), vector.end()), vector.end());
}

auto perfect_powers(uint128_t n) {
    std::vector<uint128_t> result;
    for (uint128_t i = 2, s = sqrt(n); i <= s; ++i)
        for (uint128_t p = i * i; p < n; p *= i)
            result.push_back(p);
    unique_sort(result);
    return result;
}

auto achilles(uint128_t from, uint128_t to, const std::vector<uint128_t>& pps) {
    std::vector<uint128_t> result;
    auto c = static_cast<uint128_t>(std::cbrt(static_cast<double>(to / 4)));
    auto s = sqrt(to / 8);
    for (uint128_t b = 2; b <= c; ++b) {
        uint128_t b3 = b * b * b;
        for (uint128_t a = 2; a <= s; ++a) {
            uint128_t p = b3 * a * a;
            if (p >= to)
                break;
            if (p >= from && !binary_search(pps.begin(), pps.end(), p))
                result.push_back(p);
        }
    }
    unique_sort(result);
    return result;
}

uint128_t totient(uint128_t n) {
    uint128_t tot = n;
    if ((n & 1) == 0) {
        while ((n & 1) == 0)
            n >>= 1;
        tot -= tot >> 1;
    }
    for (uint128_t p = 3; p * p <= n; p += 2) {
        if (n % p == 0) {
            while (n % p == 0)
                n /= p;
            tot -= tot / p;
        }
    }
    if (n > 1)
        tot -= tot / n;
    return tot;
}

int main() {
    auto start = std::chrono::high_resolution_clock::now();

    const uint128_t limit = 1000000000000000;

    auto pps = perfect_powers(limit);
    auto ach = achilles(1, 1000000, pps);

    std::cout << "First 50 Achilles numbers:\n";
    for (size_t i = 0; i < 50 && i < ach.size(); ++i)
        std::cout << std::setw(4) << ach[i] << ((i + 1) % 10 == 0 ? '\n' : ' ');

    std::cout << "\nFirst 50 strong Achilles numbers:\n";
    for (size_t i = 0, count = 0; count < 50 && i < ach.size(); ++i)
        if (binary_search(ach.begin(), ach.end(), totient(ach[i])))
            std::cout << std::setw(6) << ach[i]
                      << (++count % 10 == 0 ? '\n' : ' ');

    int digits = 2;
    std::cout << "\nNumber of Achilles numbers with:\n";
    for (uint128_t from = 1, to = 100; to <= limit; to *= 10, ++digits) {
        size_t count = achilles(from, to, pps).size();
        std::cout << std::setw(2) << digits << " digits: " << count << '\n';
        from = to;
    }

    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> duration(end - start);
    std::cout << "\nElapsed time: " << duration.count() << " seconds\n";
}
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

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

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

Elapsed time: 13.2644 seconds

Delphi

Works with: Delphi version 6.0


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;
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

Digits Counts:
Last Count: 100000
2 digits: 1
3 digits: 12
4 digits: 47
5 digits: 192


EasyLang

Translation of: BASIC256
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
.
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 to 5 digits:
1 12 47 192 

Factor

Works with: Factor version 0.99 2022-04-03
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
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

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
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

Translation of: Wren

Based on Version 2, takes around 19 seconds.

package main

import (
    "fmt"
    "math"
    "sort"
)

func totient(n int) int {
    tot := n
    i := 2
    for i*i <= n {
        if n%i == 0 {
            for n%i == 0 {
                n /= i
            }
            tot -= tot / i
        }
        if i == 2 {
            i = 1
        }
        i += 2
    }
    if n > 1 {
        tot -= tot / n
    }
    return tot
}

var pps = make(map[int]bool)

func getPerfectPowers(maxExp int) {
    upper := math.Pow(10, float64(maxExp))
    for i := 2; i <= int(math.Sqrt(upper)); i++ {
        fi := float64(i)
        p := fi
        for {
            p *= fi
            if p >= upper {
                break
            }
            pps[int(p)] = true
        }
    }
}

func getAchilles(minExp, maxExp int) map[int]bool {
    lower := math.Pow(10, float64(minExp))
    upper := math.Pow(10, float64(maxExp))
    achilles := make(map[int]bool)
    for b := 1; b <= int(math.Cbrt(upper)); b++ {
        b3 := b * b * b
        for a := 1; a <= int(math.Sqrt(upper)); a++ {
            p := b3 * a * a
            if p >= int(upper) {
                break
            }
            if p >= int(lower) {
                if _, ok := pps[p]; !ok {
                    achilles[p] = true
                }
            }
        }
    }
    return achilles
}

func main() {
    maxDigits := 15
    getPerfectPowers(maxDigits)
    achillesSet := getAchilles(1, 5) // enough for first 2 parts
    achilles := make([]int, len(achillesSet))
    i := 0
    for k := range achillesSet {
        achilles[i] = k
        i++
    }
    sort.Ints(achilles)

    fmt.Println("First 50 Achilles numbers:")
    for i = 0; i < 50; i++ {
        fmt.Printf("%4d ", achilles[i])
        if (i+1)%10 == 0 {
            fmt.Println()
        }
    }

    fmt.Println("\nFirst 30 strong Achilles numbers:")
    var strongAchilles []int
    count := 0
    for n := 0; count < 30; n++ {
        tot := totient(achilles[n])
        if _, ok := achillesSet[tot]; ok {
            strongAchilles = append(strongAchilles, achilles[n])
            count++
        }
    }
    for i = 0; i < 30; i++ {
        fmt.Printf("%5d ", strongAchilles[i])
        if (i+1)%10 == 0 {
            fmt.Println()
        }
    }

    fmt.Println("\nNumber of Achilles numbers with:")
    for d := 2; d <= maxDigits; d++ {
        ac := len(getAchilles(d-1, d))
        fmt.Printf("%2d digits: %d\n", d, ac)
    }
}
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800 
 864  968  972 1125 1152 1323 1352 1372 1568 1800 
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 

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

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

J

Implementation:

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

Task examples:

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

Explanation of the code:

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

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

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

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

5&p: is euler's totient function.

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

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);
        }
    }
}
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

Julia

using Primes

isAchilles(n) = (p = [x[2] for x in factor(n).pe]; all(>(1), p) && gcd(p) == 1)

isstrongAchilles(n) = isAchilles(n) && isAchilles(totient(n))

function teststrongachilles(nachilles = 50, nstrongachilles = 100)
    # task 1
    println("First $nachilles Achilles numbers:")
    n, found = 0, 0
    while found < nachilles
        if isAchilles(n)
            found += 1
            print(rpad(n, 5), found % 10 == 0 ? "\n" : "")
        end
        n += 1
    end
    # task 2
    println("\nFirst $nstrongachilles strong Achilles numbers:")
    n, found = 0, 0
    while found < nstrongachilles
        if isstrongAchilles(n)
            found += 1
            print(rpad(n, 7), found % 10 == 0 ? "\n" : "")
        end
        n += 1
    end
    # task 3
    println("\nCount of Achilles numbers for various intervals:")
    intervals = [10:99, 100:999, 1000:9999, 10000:99999, 100000:999999]
    for interval in intervals
        println(lpad(interval, 15), " ", count(isAchilles, interval))
    end
end

teststrongachilles()
Output:
First 50 Achilles numbers:
72   108  200  288  392  432  500  648  675  800
864  968  972  1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

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

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

Mathematica/Wolfram Language

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
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

Nim

Translation of: Wren
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}"
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

Perl

Borrowed, and lightly modified, code from Powerful_numbers

Library: ntheory
use strict;
use warnings;
use feature <say current_sub>;
use experimental 'signatures';
use List::AllUtils <max head>;
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";
}
Output:
First 50 Achilles numbers:
   72  108  200  288  392  432  500  648  675  800
  864  968  972 1125 1152 1323 1352 1372 1568 1800
 1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

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

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

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

use strict;
use warnings;

my %pps;
my $maxDigits = 9;

sub totient { 
   my $tot = my $n = shift; 
   my $i   = 2;
   while ($i*$i <= $n) {
      unless ($n % $i) {
         until($n % $i) { $n = int($n/$i) }
         $tot -= int($tot/$i)
      }
      if ($i == 2) { $i = 1 }
      $i += 2; 
   }
   if ($n > 1) { $tot -= int($tot/$n) }
   return $tot
}

sub getPerfectPowers {
   for my $i (2..int(sqrt(my $upper = 10**( shift )))) {
      my $p = $i;
      while (($p *= $i) < $upper) { $pps{$p}++ }
   }
}

sub getAchilles { 
   my ($lower, $upper) = map { 10** $_ } @_ ;
   my %achilles = (); 
   my $count = 0;
   for my $b (1..int($upper**(1/3))) {
      my ($b3,$p) = $b * $b * $b;
      for my $a (1..int(sqrt($upper))) {
         last if (($p = $b3 * $a * $a) >= $upper);
         $achilles{$p}++ if ($p >= $lower and !$pps{$p})  
      }
   }
   return keys %achilles
}

getPerfectPowers $maxDigits;

my @achilles = sort { $a <=> $b } getAchilles(1,5);
my %achillesSet;
@achillesSet{ @achilles } = undef; 

print "First 50 Achilles numbers:\n";
for (0..49) { printf "%5d".($_%10 == 9 ? "\n" : " "), $achilles[$_] }

my %strongAchilles;
my $count = my $n = 0;
for (my $count = my $n = 0; $count < 30; $n++) {
   if ( exists($achillesSet{ totient( $achilles[$n] ) })) {
      $strongAchilles{ $achilles[$n] }++;
      $count++
   }
}

my @strongAchilles30 = (sort { $a <=> $b } keys %strongAchilles)[0..29];

print "\nFirst 30 strong Achilles numbers:\n";
for (0..29) { printf "%5d".($_%10 == 9 ? "\n" : " "), $strongAchilles30[$_] }
 
print "\nNumber of Achilles numbers with:\n";
for my $d (2..$maxDigits) {
   printf "%2d digits: %d\n", $d, scalar getAchilles($d-1, $d)
}

Output is the same.

Phix

Library: Phix/online

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

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

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

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

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

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

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

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)
Output:
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

Quackery

gcd is defined at Greatest common divisor#Quackery.

primefactors is defined at Prime decomposition#Quackery.

totient is defined at Totient function#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
Output:
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

Raku

Timing is going to be system / OS dependent.

use Prime::Factor;
use Math::Root;

sub is-square-free (Int \n) {
    constant @p = ^100 .map: { next unless .is-prime; .ÂČ };
    for @p -> \p { return False if n %% p }
    True
}

sub powerful (\n, \k = 2) {
    my @p;
    p(1, 2*k - 1);
    sub p (\m, \r) {
        @p.push(m) and return if r < k;
        for 1 .. (n / m).&root(r) -> \v {
            if r > k {
                next unless is-square-free(v);
                next unless m gcd v == 1;
            }
            p(m * v ** r, r - 1)
        }
    }
    @p
}

my $achilles = powerful(10**9).hyper(:500batch).grep( {
    my $f = .&prime-factors.Bag;
    (+$f.keys > 1) && (1 == [gcd] $f.values) && (.sqrt.IntÂČ !== $_)
} ).classify: { .chars }

my \𝜑 = 0, |(1..*).hyper.map: -> \t { t × [×] t.&prime-factors.squish.map: { 1 - 1/$_ } }

my %as = Set.new: flat $achilles.values».list;

my $strong = lazy (flat $achilles.sort».value».list».sort).grep: { ?%as{𝜑[$_]} };

put "First 50 Achilles numbers:";
put (flat $achilles.sort».value».list».sort)[^50].batch(10)».fmt("%4d").join("\n");

put "\nFirst 30 strong Achilles numbers:";
put   $strong[^30].batch(10)».fmt("%5d").join("\n");

put "\nNumber of Achilles numbers with:";
say "$_ digits: " ~ +$achilles{$_} // 0 for 2..9;

printf "\n%.1f total elapsed seconds\n", now - INIT now;
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

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

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

6.1 total elapsed seconds

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

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

410.4 total elapsed seconds

RPL

Based on Wikipedia definition: n = p1a1p2a2
pkak is an Achilles number if min(a1, a2, 
, ak) ≄ 2 and gcd(a1, a2, 
, ak) = 1.

Works with: HP version 49g
â‰Ș FACTORS
   IF DUP SIZE 4 < THEN DROP 0 
   ELSE
      { } 2 3 PICK SIZE FOR j          @ keeps exponents only 
          OVER j GET +
      2 STEP NIP
      → a                                         
      â‰Ș a â‰Ș MIN ≫ STREAM 2 ≄ a â‰Ș GCD ≫ STREAM 1 == AND        
      ≫
   END
≫ 'ACH?' STO 

â‰Ș { } 1 → n 
  â‰Ș WHILE DUP SIZE 50 < REPEAT
        IF 'n' INCR ACH? THEN n + END
     END
≫ ≫ EVAL

â‰Ș { } 1 → n 
  â‰Ș WHILE DUP SIZE 20 < REPEAT
        IF 'n' INCR ACH? 
           THEN IF n EULER ACH? THEN n + END END
     END
≫ ≫ EVAL
Output:
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}

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.

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}" }
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

Rust

Translation of: Wren
fn perfect_powers(n: u128) -> Vec<u128> {
    let mut powers = Vec::<u128>::new();
    let sqrt = (n as f64).sqrt() as u128;
    for i in 2..=sqrt {
        let mut p = i * i;
        while p < n {
            powers.push(p);
            p *= i;
        }
    }
    powers.sort();
    powers.dedup();
    powers
}

fn bsearch<T: Ord>(vector: &Vec<T>, value: &T) -> bool {
    match vector.binary_search(value) {
        Ok(_) => true,
        _ => false,
    }
}

fn achilles(from: u128, to: u128, pps: &Vec<u128>) -> Vec<u128> {
    let mut result = Vec::<u128>::new();
    let cbrt = ((to / 4) as f64).cbrt() as u128;
    let sqrt = ((to / 8) as f64).sqrt() as u128;
    for b in 2..=cbrt {
        let b3 = b * b * b;
        for a in 2..=sqrt {
            let p = b3 * a * a;
            if p >= to {
                break;
            }
            if p >= from && !bsearch(&pps, &p) {
                result.push(p);
            }
        }
    }
    result.sort();
    result.dedup();
    result
}

fn totient(mut n: u128) -> u128 {
    let mut tot = n;
    if (n & 1) == 0 {
        while (n & 1) == 0 {
            n >>= 1;
        }
        tot -= tot >> 1;
    }
    let mut p = 3;
    while p * p <= n {
        if n % p == 0 {
            while n % p == 0 {
                n /= p;
            }
            tot -= tot / p;
        }
        p += 2;
    }
    if n > 1 {
        tot -= tot / n;
    }
    tot
}

fn main() {
    use std::time::Instant;
    let t0 = Instant::now();
    let limit = 1000000000000000u128;

    let pps = perfect_powers(limit);
    let ach = achilles(1, 1000000, &pps);

    println!("First 50 Achilles numbers:");
    for i in 0..50 {
        print!("{:4}{}", ach[i], if (i + 1) % 10 == 0 { "\n" } else { " " });
    }

    println!("\nFirst 50 strong Achilles numbers:");
    for (i, n) in ach
        .iter()
        .filter(|&x| bsearch(&ach, &totient(*x)))
        .take(50)
        .enumerate()
    {
        print!("{:6}{}", n, if (i + 1) % 10 == 0 { "\n" } else { " " });
    }
    println!();

    let mut from = 1u128;
    let mut to = 100u128;
    let mut digits = 2;
    while to <= limit {
        let count = achilles(from, to, &pps).len();
        println!("{:2} digits: {}", digits, count);
        from = to;
        to *= 10;
        digits += 1;
    }

    let duration = t0.elapsed();
    println!("\nElapsed time: {} milliseconds", duration.as_millis());
}
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

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

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

Elapsed time: 12608 milliseconds

Sidef

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}"
}
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

Wren

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

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 "./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 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 = Int.totient(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
}
Fmt.tprint("$4d", achilles, 10)

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
}
Fmt.tprint("$5d", strongAchilles, 10)

System.print("\nNumber of Achilles numbers with:")
var pow = 10
for (i in 2..maxDigits) {
    var count = 0
    for (j in pow..pow*10-1) {
        if (isAchilles.call(j)) count = count + 1
    }
    System.print("%(i) digits: %(count)")
    pow = pow * 10
}
Output:
First 50 Achilles numbers:
  72  108  200  288  392  432  500  648  675  800
 864  968  972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200

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

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

Version 2 (Much faster)

Library: Wren-set

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

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

import "./set" for Set
import "./seq" for Lst
import "./math" for Int
import "./fmt" for Fmt

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:")
Fmt.tprint("$4d", achilles.take(50), 10)

System.print("\nFirst 30 strong Achilles numbers:")
var strongAchilles = []
var count = 0
var n = 0
while (count < 30) {
    var tot = Int.totient(achilles[n])
    if (achillesSet.contains(tot)) {
        strongAchilles.add(achilles[n])
        count = count + 1
    }
    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)
}
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

func GCD(N, D);         \Return the greatest common divisor of N and D
int  N, D;              \numerator and denominator
int  R;
[if D > N then
    [R:= D;  D:= N;  N:= R];    \swap D and N
while D > 0 do
    [R:= rem(N/D);
    N:= D;
    D:= R;
    ];
return N;
];      \GCD

func Totient(N);        \Return the totient of N
int  N, Phi, M;
[Phi:= 0;
for M:= 1 to N do
    if GCD(M, N) = 1 then Phi:= Phi+1;
return Phi;
];

func Powerful(N0);      \Return 'true' if N0 is a powerful number
int  N0, N, F, Q, L;
[if N0 <= 1 then return false;
N:= N0;  F:= 2;
L:= sqrt(N0);
loop    [Q:= N/F;
        if rem(0) = 0 then      \found a factor
                [if rem(N0/(F*F)) then return false;
                N:= Q;
                if F>N then quit;
                ]
        else    [F:= F+1;
                if F > L then
                    [if rem(N0/(N*N)) then return false;
                    quit;
                    ];
                ];
        ];
return true;
];

func Achilles(N);       \Return 'true' if N is an Achilles number
int  N, M, A;
[if not Powerful(N) then return false;
M:= 2;
A:= M*M;
repeat  loop    [if A = N then return false;
                if A > N then quit;
                A:= A*M;
                ];
        M:= M+1;
        A:= M*M;
until   A > N;
return true;
];

int Cnt, N, Pwr, Start;
[Cnt:= 0;
N:= 1;
loop    [if Achilles(N) then
            [IntOut(0, N);
            Cnt:= Cnt+1;
            if Cnt >= 50 then quit;
            if rem(Cnt/10) then ChOut(0, 9) else CrLf(0);
            ];
        N:= N+1;
        ];
CrLf(0);  CrLf(0);
Cnt:= 0;
N:= 1;
loop    [if Achilles(N) then
            if Achilles(Totient(N)) then
                [IntOut(0, N);
                Cnt:= Cnt+1;
                if Cnt >= 20 then quit;
                if rem(Cnt/10) then ChOut(0, 9) else CrLf(0);
                ];
        N:= N+1;
        ];
CrLf(0);  CrLf(0);
for Pwr:= 1 to 6 do
    [IntOut(0, Pwr);  Text(0, ": ");
    Start:= fix(Pow(10.0, float(Pwr-1)));
    Cnt:= 0;
    for N:= Start to Start*10-1 do
        if Achilles(N) then Cnt:= Cnt+1;
    IntOut(0, Cnt);  CrLf(0);
    ];
]
Output:
72      108     200     288     392     432     500     648     675     800
864     968     972     1125    1152    1323    1352    1372    1568    1800
1944    2000    2312    2592    2700    2888    3087    3200    3267    3456
3528    3872    3888    4000    4232    4500    4563    4608    5000    5292
5324    5400    5408    5488    6075    6125    6272    6728    6912    7200

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

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