Totient function
You are encouraged to solve this task according to the task description, using any language you may know.
The totient function is also known as:
- Euler's totient function
- Euler's phi totient function
- phi totient function
- Φ function (uppercase Greek phi)
- φ function (lowercase Greek phi)
- Definitions (as per number theory)
The totient function:
- counts the integers up to a given positive integer n that are relatively prime to n
- counts the integers k in the range 1 ≤ k ≤ n for which the greatest common divisor gcd(n,k) is equal to 1
- counts numbers ≤ n and prime to n
If the totient number (for N) is one less than N, then N is prime.
- Task
Create a totient function and:
- Find and display (1 per line) for the 1st 25 integers:
- the integer (the index)
- the totient number for that integer
- indicate if that integer is prime
- Find and display the count of the primes up to 100
- Find and display the count of the primes up to 1,000
- Find and display the count of the primes up to 10,000
- Find and display the count of the primes up to 100,000 (optional)
Show all output here.
- Related task
- Also see
11l
F f(n)
R sum((1..n).filter(k -> gcd(@n, k) == 1).map(k -> 1))
F is_prime(n)
R f(n) == n - 1
L(n) 1..25
print(‘ f(#.) == #.’.format(n, f(n))‘’(I is_prime(n) {‘, is prime’} E ‘’))
V count = 0
L(n) 1..10'000
count += is_prime(n)
I n C (100, 1000, 10'000)
print(‘Primes up to #.: #.’.format(n, count))- Output:
f(1) == 1 f(2) == 1, is prime f(3) == 2, is prime f(4) == 2 f(5) == 4, is prime f(6) == 2 f(7) == 6, is prime f(8) == 4 f(9) == 6 f(10) == 4 f(11) == 10, is prime f(12) == 4 f(13) == 12, is prime f(14) == 6 f(15) == 8 f(16) == 8 f(17) == 16, is prime f(18) == 6 f(19) == 18, is prime f(20) == 8 f(21) == 12 f(22) == 10 f(23) == 22, is prime f(24) == 8 f(25) == 20 Primes up to 100: 25 Primes up to 1000: 168 Primes up to 10000: 1229
AArch64 Assembly
/* ARM assembly AARCH64 Raspberry PI 3B */
/* program totient.s */
/************************************/
/* Constantes */
/************************************/
.include "../includeConstantesARM64.inc"
.equ MAXI, 25
/*********************************/
/* Initialized data */
/*********************************/
.data
szMessNumber: .asciz " number @ totient @ @ \n"
szCarriageReturn: .asciz "\n"
szMessPrime: .asciz " is prime."
szMessSpace: .asciz " "
szMessCounterPrime: .asciz "Number of primes to @ : @ \n"
szMessOverflow: .asciz "Overflow function isPrime.\n"
/*********************************/
/* UnInitialized data */
/*********************************/
.bss
sZoneConv: .skip 24
/*********************************/
/* code section */
/*********************************/
.text
.global main
main:
mov x4,#1 // start number
1:
mov x0,x4
bl totient // compute totient
mov x5,x0
mov x0,x4
bl isPrime // control if number is prime
mov x6,x0
mov x0,x4 // display result
ldr x1,qAdrsZoneConv
bl conversion10 // call décimal conversion
ldr x0,qAdrszMessNumber
ldr x1,qAdrsZoneConv // insert conversion in message
bl strInsertAtCharInc
mov x7,x0
mov x0,x5
ldr x1,qAdrsZoneConv
bl conversion10 // call décimal conversion
mov x0,x7
ldr x1,qAdrsZoneConv // insert conversion in message
bl strInsertAtCharInc
mov x7,x0
cmp x6,#1
ldr x1,qAdrszMessPrime
ldr x8,qAdrszMessSpace
csel x1,x1,x8,eq
mov x0,x7
bl strInsertAtCharInc
bl affichageMess // display message
add x4,x4,#1 // increment number
cmp x4,#MAXI // maxi ?
ble 1b // and loop
mov x4,#2 // first number
mov x5,#0 // prime counter
ldr x6,iCst1000 // load constantes
ldr x7,iCst10000
ldr x8,iCst100000
2:
mov x0,x4
bl isPrime
cmp x0,#0
beq 3f
add x5,x5,#1
3:
add x4,x4,#1
cmp x4,#100
bne 4f
mov x0,#100
mov x1,x5
bl displayCounter
b 7f
4:
cmp x4,x6 // 1000
bne 5f
mov x0,x6
mov x1,x5
bl displayCounter
b 7f
5:
cmp x4,x7 // 10000
bne 6f
mov x0,x7
mov x1,x5
bl displayCounter
b 7f
6:
cmp x4,x8 // 100000
bne 7f
mov x0,x8
mov x1,x5
bl displayCounter
7:
cmp x4,x8
ble 2b // and loop
100: // standard end of the program
mov x0, #0 // return code
mov x8,EXIT
svc #0 // perform the system call
qAdrszCarriageReturn: .quad szCarriageReturn
qAdrsZoneConv: .quad sZoneConv
qAdrszMessNumber: .quad szMessNumber
qAdrszMessCounterPrime: .quad szMessCounterPrime
qAdrszMessPrime: .quad szMessPrime
qAdrszMessSpace: .quad szMessSpace
iCst1000: .quad 1000
iCst10000: .quad 10000
iCst100000: .quad 100000
/******************************************************************/
/* display 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,qAdrszMessCounterPrime
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
100:
ldp x2,x3,[sp],16 // restaur registers
ldp x1,lr,[sp],16 // restaur registers
ret
/******************************************************************/
/* compute totient of number */
/******************************************************************/
/* x0 contains number */
totient:
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
/***************************************************/
/* 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"number 1 totient 1 is prime. number 2 totient 1 is prime. number 3 totient 2 is prime. number 4 totient 2 number 5 totient 4 is prime. number 6 totient 2 number 7 totient 6 is prime. number 8 totient 4 number 9 totient 6 number 10 totient 4 number 11 totient 10 is prime. number 12 totient 4 number 13 totient 12 is prime. number 14 totient 6 number 15 totient 8 number 16 totient 8 number 17 totient 16 is prime. number 18 totient 6 number 19 totient 18 is prime. number 20 totient 8 number 21 totient 12 number 22 totient 10 number 23 totient 22 is prime. number 24 totient 8 number 25 totient 20 Number of primes to 100 : 25 Number of primes to 1000 : 168 Number of primes to 10000 : 1229 Number of primes to 100000 : 9592
Ada
with Ada.Text_IO; with Ada.Integer_Text_IO;
with Ada.Strings.Fixed;
procedure Totient is
function Totient (N : in Integer) return Integer is
Tot : Integer := N;
I : Integer;
N2 : Integer := N;
begin
I := 2;
while I * I <= N2 loop
if N2 mod I = 0 then
while N2 mod I = 0 loop
N2 := N2 / I;
end loop;
Tot := Tot - Tot / I;
end if;
if I = 2 then
I := 1;
end if;
I := I + 2;
end loop;
if N2 > 1 then
Tot := Tot - Tot / N2;
end if;
return Tot;
end Totient;
Count : Integer := 0;
Tot : Integer;
Placeholder : String := " n Phi Is_Prime";
Image_N : String renames Placeholder ( 1 .. 3);
Image_Phi : String renames Placeholder ( 6 .. 8);
Image_Prime : String renames Placeholder (11 .. Placeholder'Last);
use Ada.Text_IO;
use Ada.Integer_Text_IO;
Last : constant := 25;
begin
Put_Line (Placeholder);
for N in 1 .. Last loop
Tot := Totient (N);
if N - 1 = Tot then
Count := Count + 1;
end if;
Put (Image_N, N);
Put (Image_Phi, Tot);
Image_Prime := Ada.Strings.Fixed.Tail
(Source => (if N - 1 = Tot then True'Image else False'Image),
Count => Image_Prime'Length);
Put_Line (Placeholder);
end loop;
New_Line;
Put_Line ("Number of primes up to " & Last'Image &" =" & Count'Image);
for N in Last + 1 .. 100_000 loop
Tot := Totient (N);
if Tot = N - 1 then
Count := Count + 1;
end if;
if N = 100 or N = 1_000 or N mod 10_000 = 0 then
Put_Line ("Number of primes up to " & N'Image & " =" & Count'Image);
end if;
end loop;
end Totient;
- Output:
n Phi Is_Prime 1 1 FALSE 2 1 TRUE 3 2 TRUE 4 2 FALSE 5 4 TRUE 6 2 FALSE 7 6 TRUE 8 4 FALSE 9 6 FALSE 10 4 FALSE 11 10 TRUE 12 4 FALSE 13 12 TRUE 14 6 FALSE 15 8 FALSE 16 8 FALSE 17 16 TRUE 18 6 FALSE 19 18 TRUE 20 8 FALSE 21 12 FALSE 22 10 FALSE 23 22 TRUE 24 8 FALSE 25 20 FALSE Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 20000 = 2262 Number of primes up to 30000 = 3245 Number of primes up to 40000 = 4203 Number of primes up to 50000 = 5133 Number of primes up to 60000 = 6057 Number of primes up to 70000 = 6935 Number of primes up to 80000 = 7837 Number of primes up to 90000 = 8713 Number of primes up to 100000 = 9592
ALGOL 68
Uses the totient algorithm from the second Go sample.
BEGIN
# returns the number of integers k where 1 <= k <= n that are mutually prime to n #
PROC totient = ( INT n )INT:
IF n < 3 THEN 1
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 # ;
# show the totient function values for the first 25 integers #
print( ( " n phi(n) remarks", newline ) );
FOR n TO 25 DO
INT tn = totient( n );
print( ( whole( n, -2 ), ": ", whole( tn, -5 )
, IF tn = n - 1 AND tn /= 0 THEN " n is prime" ELSE "" FI, newline
)
)
OD;
# use the totient function to count primes #
INT n100 := 0, n1000 := 0, n10000 := 0, n100000 := 0;
FOR n TO 100 000 DO
IF totient( n ) = n - 1 THEN
IF n <= 100 THEN n100 +:= 1 FI;
IF n <= 1 000 THEN n1000 +:= 1 FI;
IF n <= 10 000 THEN n10000 +:= 1 FI;
IF n <= 100 000 THEN n100000 +:= 1 FI
FI
OD;
print( ( "There are ", whole( n100, -6 ), " primes below 100", newline ) );
print( ( "There are ", whole( n1000, -6 ), " primes below 1 000", newline ) );
print( ( "There are ", whole( n10000, -6 ), " primes below 10 000", newline ) );
print( ( "There are ", whole( n100000, -6 ), " primes below 100 000", newline ) )
END- Output:
n phi(n) remarks 1: 1 2: 1 n is prime 3: 2 n is prime 4: 2 5: 4 n is prime 6: 2 7: 6 n is prime 8: 4 9: 6 10: 4 11: 10 n is prime 12: 4 13: 12 n is prime 14: 6 15: 8 16: 8 17: 16 n is prime 18: 6 19: 18 n is prime 20: 8 21: 12 22: 10 23: 22 n is prime 24: 8 25: 20 There are 25 primes below 100 There are 168 primes below 1 000 There are 1229 primes below 10 000 There are 9592 primes below 100 000
ALGOL-M
The GCD approach is straight-forward and easy to follow, but is not particularly efficient.
BEGIN
% RETURN P MOD Q %
INTEGER FUNCTION MOD (P, Q);
INTEGER P, Q;
BEGIN
MOD := P - Q * (P / Q);
END;
% RETURN GREATEST COMMON DIVISOR OF X AND Y %
INTEGER FUNCTION GCD (X, Y);
INTEGER X, Y;
BEGIN
INTEGER R;
IF X < Y THEN
BEGIN
INTEGER TEMP;
TEMP := X;
X := Y;
Y := TEMP;
END;
WHILE (R := MOD(X, Y)) <> 0 DO
BEGIN
X := Y;
Y := R;
END;
GCD := Y;
END;
% RETURN PHI (ALSO CALLED TOTIENT) OF N %
INTEGER FUNCTION PHI(N);
INTEGER N;
BEGIN
INTEGER I, COUNT;
COUNT := 1;
FOR I := 2 STEP 1 UNTIL N DO
BEGIN
IF GCD(N,I) = 1 THEN COUNT := COUNT + 1;
END;
PHI := COUNT;
END;
COMMENT - EXERCISE THE FUNCTION;
INTEGER N, TOTIENT, COUNT;
WRITE(" N PHI(N) PRIME?");
FOR N := 1 STEP 1 UNTIL 25 DO
BEGIN
WRITE(N, (TOTIENT := PHI(N)));
WRITEON(IF TOTIENT = (N-1) THEN " YES" ELSE " NO");
END;
COMMENT - AND USE IT TO COUNT PRIMES;
WRITE("");
COUNT := 0;
FOR N := 1 STEP 1 UNTIL 1000 DO
BEGIN
IF PHI(N) = (N-1) THEN COUNT := COUNT + 1;
IF N = 100 THEN
WRITE("PRIMES UP TO 100 =", COUNT)
ELSE IF N = 1000 THEN
WRITE("PRIMES UP TO 1000 =", COUNT);
END;
END- Output:
Limiting the search to 1000 keeps the running time of the program within reasonable bounds.
N PHI(N) PRIME?
1 1 NO
2 1 YES
3 2 YES
4 2 NO
5 4 YES
6 2 NO
7 6 YES
8 4 NO
9 6 NO
10 4 NO
11 10 YES
12 4 NO
13 12 YES
14 6 NO
15 8 NO
16 8 NO
17 16 YES
18 6 NO
19 18 YES
20 8 NO
21 12 NO
22 10 NO
23 22 YES
24 8 NO
25 20 NO
PRIMES UP TO 100 = 25
PRIMES UP TO 1000 = 168
APL
task←{
totient ← 1+.=⍳∨⊢
prime ← totient=-∘1
⎕←'Index' 'Totient' 'Prime',(⊢⍪totient¨,[÷2]prime¨)⍳25
{⎕←'There are' (+/prime¨⍳⍵) 'primes below' ⍵}¨100 1000 10000
}
- Output:
Index 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 Totient 1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22 8 20 Prime 0 1 1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0 There are 25 primes below 100 There are 168 primes below 1000 There are 1229 primes below 10000
ARM Assembly
/* ARM assembly Raspberry PI or android with termux */
/* program totient.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 MAXI, 25
/*********************************/
/* Initialized data */
/*********************************/
.data
szMessNumber: .asciz " number @ totient @ @ \n"
szCarriageReturn: .asciz "\n"
szMessPrime: .asciz " is prime."
szMessSpace: .asciz " "
szMessCounterPrime: .asciz "Number of primes to @ : @ \n"
/*********************************/
/* UnInitialized data */
/*********************************/
.bss
sZoneConv: .skip 24
/*********************************/
/* code section */
/*********************************/
.text
.global main
main:
mov r4,#1 @ start number
1:
mov r0,r4
bl totient @ compute totient
mov r5,r0
mov r0,r4
bl isPrime @ control if number is prime
mov r6,r0
mov r0,r4 @ display result
ldr r1,iAdrsZoneConv
bl conversion10 @ call décimal conversion
ldr r0,iAdrszMessNumber
ldr r1,iAdrsZoneConv @ insert conversion in message
bl strInsertAtCharInc
mov r7,r0
mov r0,r5
ldr r1,iAdrsZoneConv
bl conversion10 @ call décimal conversion
mov r0,r7
ldr r1,iAdrsZoneConv @ insert conversion in message
bl strInsertAtCharInc
mov r7,r0
cmp r6,#1
ldreq r1,iAdrszMessPrime
ldrne r1,iAdrszMessSpace
mov r0,r7
bl strInsertAtCharInc
bl affichageMess @ display message
add r4,r4,#1 @ increment number
cmp r4,#MAXI @ maxi ?
ble 1b @ and loop
mov r4,#2 @ first number
mov r5,#0 @ prime counter
ldr r6,iCst1000 @ load constantes
ldr r7,iCst10000
ldr r8,iCst100000
2:
mov r0,r4
bl isPrime
cmp r0,#0
beq 3f
add r5,r5,#1
3:
add r4,r4,#1
cmp r4,#100
bne 4f
mov r0,#100
mov r1,r5
bl displayCounter
b 7f
4:
cmp r4,r6 @ 1000
bne 5f
mov r0,r6
mov r1,r5
bl displayCounter
b 7f
5:
cmp r4,r7 @ 10000
bne 6f
mov r0,r7
mov r1,r5
bl displayCounter
b 7f
6:
cmp r4,r8 @ 100000
bne 7f
mov r0,r8
mov r1,r5
bl displayCounter
7:
cmp r4,r8
ble 2b @ and loop
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
iAdrsZoneConv: .int sZoneConv
iAdrszMessNumber: .int szMessNumber
iAdrszMessCounterPrime: .int szMessCounterPrime
iAdrszMessPrime: .int szMessPrime
iAdrszMessSpace: .int szMessSpace
iCst1000: .int 1000
iCst10000: .int 10000
iCst100000: .int 100000
/******************************************************************/
/* display 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,iAdrszMessCounterPrime
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
100:
pop {r1-r3,pc} @ restaur registers
/******************************************************************/
/* compute totient of number */
/******************************************************************/
/* r0 contains number */
totient:
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
/***************************************************/
/* 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"number 1 totient 1 is prime. number 2 totient 1 is prime. number 3 totient 2 is prime. number 4 totient 2 number 5 totient 4 is prime. number 6 totient 2 number 7 totient 6 is prime. number 8 totient 4 number 9 totient 6 number 10 totient 4 number 11 totient 10 is prime. number 12 totient 4 number 13 totient 12 is prime. number 14 totient 6 number 15 totient 8 number 16 totient 8 number 17 totient 16 is prime. number 18 totient 6 number 19 totient 18 is prime. number 20 totient 8 number 21 totient 12 number 22 totient 10 number 23 totient 22 is prime. number 24 totient 8 number 25 totient 20 Number of primes to 100 : 25 Number of primes to 1000 : 168 Number of primes to 10000 : 1229 Number of primes to 100000 : 9592
Arturo
totient: function [n][
tot: new n
i: 2
while -> n >= i * i [
if 0 = n % i [
while -> 0 = n % i -> n: n / i
'tot - tot / i
]
if 2 = i -> i: 1
'i + 2
]
if n > 1 -> 'tot - tot / n
return tot
]
primes: 0
loop 1..100000 'i [
t: totient i
prime?: 1 = i - t
if i < 26 [
prints ~« Φ(|pad.with:`0` to :string i 2|) = |pad to :string t 2|
if prime? -> prints ", prime"
print ""
]
if 50 = i -> print ""
if in? i [100 1000 10000 100000] -> print ~« |primes| primes =< |i|
if prime? -> 'primes + 1
]
- Output:
Φ(01) = 1 Φ(02) = 1, prime Φ(03) = 2, prime Φ(04) = 2 Φ(05) = 4, prime Φ(06) = 2 Φ(07) = 6, prime Φ(08) = 4 Φ(09) = 6 Φ(10) = 4 Φ(11) = 10, prime Φ(12) = 4 Φ(13) = 12, prime Φ(14) = 6 Φ(15) = 8 Φ(16) = 8 Φ(17) = 16, prime Φ(18) = 6 Φ(19) = 18, prime Φ(20) = 8 Φ(21) = 12 Φ(22) = 10 Φ(23) = 22, prime Φ(24) = 8 Φ(25) = 20 25 primes =< 100 168 primes =< 1000 1229 primes =< 10000 9592 primes =< 100000
AutoHotkey
global cptext := a_tab "Nr" a_tab "Phi" a_tab "Prime?`n---------------------------------------------`n"
divisores(num)
{
serie := ""
loop % num
if !mod(num,a_index)
serie .= a_index ","
return serie
}
gcd(serieA,serieB)
{
emComum := 0
loop,parse,serieA,csv
if A_LoopField in %serieB%
emComum += 1
return emComum
}
principal(voltas,phi:=0)
{
loop %voltas%
{
cp := A_Index
cpcount := 0
numA := divisores(cp)
loop % a_index
{
numA := divisores(cp)
numB := divisores(A_Index)
fim := gcd(numA,numB)
if (fim = 1)
cpcount += 1
}
if (cpcount = cp-1)
{
if phi
cptext .= a_tab cp a_tab cpcount a_tab "1`n"
totalPrimes += 1
}
else
cptext .= a_tab cp a_tab cpcount a_tab "0`n"
}
return totalPrimes
}
totalPrimes := principal(25,1)
msgbox % cptext "`n`ntotal primes = " totalPrimes ; Number 1 is a prime number ? If yes, add 1 to totalPrimes
totalPrimes := principal(100)
msgbox % "total primes in 1 .. 100 = " totalPrimes
totalPrimes := principal(1000) ; caution... pure gcd method
msgbox % "total primes in 1 .. 1000 = " totalPrimes ; takes 3 minutes or more
;totalPrimes := principal(10000)
;msgbox % "total primes in 1 .. 10000 = " totalPrimes
ExitApp
- Output:
--------------------------- Totient function.ahk --------------------------- Nr Phi Prime? --------------------------------------------- 1 1 0 2 1 1 3 2 1 4 2 0 5 4 1 6 2 0 7 6 1 8 4 0 9 6 0 10 4 0 11 10 1 12 4 0 13 12 1 14 6 0 15 8 0 16 8 0 17 16 1 18 6 0 19 18 1 20 8 0 21 12 0 22 10 0 23 22 1 24 8 0 25 20 0 total primes = 9 --------------------------- OK --------------------------- total primes in 1 .. 100 = 25 --------------------------- OK --------------------------- total primes in 1 .. 1000 = 168 --------------------------- OK
AWK
# syntax: GAWK -f TOTIENT_FUNCTION.AWK
BEGIN {
print(" N Phi isPrime")
for (n=1; n<=1000000; n++) {
tot = totient(n)
if (n-1 == tot) {
count++
}
if (n <= 25) {
printf("%2d %3d %s\n",n,tot,(n-1==tot)?"true":"false")
if (n == 25) {
printf("\n Limit PrimeCount\n")
printf("%7d %10d\n",n,count)
}
}
else if (n ~ /^100+$/) {
printf("%7d %10d\n",n,count)
}
}
exit(0)
}
function totient(n, i,tot) {
tot = n
for (i=2; i*i<=n; i+=2) {
if (n % i == 0) {
while (n % i == 0) {
n /= i
}
tot -= tot / i
}
if (i == 2) {
i = 1
}
}
if (n > 1) {
tot -= tot / n
}
return(tot)
}
- Output:
N Phi isPrime
1 1 false
2 1 true
3 2 true
4 2 false
5 4 true
6 2 false
7 6 true
8 4 false
9 6 false
10 4 false
11 10 true
12 4 false
13 12 true
14 6 false
15 8 false
16 8 false
17 16 true
18 6 false
19 18 true
20 8 false
21 12 false
22 10 false
23 22 true
24 8 false
25 20 false
Limit PrimeCount
25 9
100 25
1000 168
10000 1229
100000 9592
1000000 78498
BASIC
10 DEFINT A-Z: DIM B$(2): B$(0)="No": B$(1)="Yes"
20 PRINT " N Totient Prime"
30 FOR N=1 TO 25
40 GOSUB 200
50 P=-(T=N-1)
60 C=C+P
70 PRINT USING "## ####### \ \";N;T;B$(P)
80 NEXT N
90 F$="Number of primes up to ######: ####"
100 PRINT USING F$;25;C
110 FOR N=26 TO 10000
120 GOSUB 200
130 C=C-(T=N-1)
140 IF N=100 OR N=1000 OR N=10000 THEN PRINT USING F$;N;C
150 NEXT N
160 END
200 REM T = TOTIENT(N)
210 T=N: Z=N
220 I=2: GOSUB 270
230 I=3
240 IF I*I<=Z THEN GOSUB 270: I=I+2: GOTO 240
250 IF Z>1 THEN T=T-T\Z
260 RETURN
270 IF Z MOD I<>0 THEN RETURN
280 IF Z MOD I=0 THEN Z = Z\I: GOTO 280
290 T = T-T\I
300 RETURN
- Output:
N Totient Prime 1 1 No 2 1 Yes 3 2 Yes 4 2 No 5 4 Yes 6 2 No 7 6 Yes 8 4 No 9 6 No 10 4 No 11 10 Yes 12 4 No 13 12 Yes 14 6 No 15 8 No 16 8 No 17 16 Yes 18 6 No 19 18 Yes 20 8 No 21 12 No 22 10 No 23 22 Yes 24 8 No 25 20 No Number of primes up to 25: 9 Number of primes up to 100: 25 Number of primes up to 1000: 168 Number of primes up to 10000: 1229
BCPL
get "libhdr"
let totient(n) = valof
$( let tot = n and i = 2
while i*i <= n
$( if n rem i = 0
$( while n rem i = 0 do n := n / i
tot := tot - tot/i
$)
if i=2 then i:=1
i := i+2
$)
resultis n>1 -> tot-tot/n, tot
$)
let start() be
$( let count = 0
writes(" N Totient Prime*N")
for n = 1 to 25
$( let tot = totient(n)
let prime = n-1 = tot
if prime then count := count + 1
writef("%I2 %I7 %S*N", n, tot, prime -> " Yes", " No")
$)
writef("Number of primes up to %I6: %I4*N", 25, count)
for n = 26 to 10000
$( if totient(n) = n-1 then count := count + 1
if n = 100 | n = 1000 | n = 10000 then
writef("Number of primes up to %I6: %I4*N", n, count)
$)
$)- Output:
N Totient Prime 1 1 No 2 1 Yes 3 2 Yes 4 2 No 5 4 Yes 6 2 No 7 6 Yes 8 4 No 9 6 No 10 4 No 11 10 Yes 12 4 No 13 12 Yes 14 6 No 15 8 No 16 8 No 17 16 Yes 18 6 No 19 18 Yes 20 8 No 21 12 No 22 10 No 23 22 Yes 24 8 No 25 20 No Number of primes up to 25: 9 Number of primes up to 100: 25 Number of primes up to 1000: 168 Number of primes up to 10000: 1229
BQN
GCD function is taken from BQNcrate.
The totient function is similar to APL and J, except it is made as a train. An explicit version of Totient is {+´1=𝕩GCD¨1+↕𝕩}
GCD ← {𝕨(|𝕊⍟(>⟜0)⊣)𝕩}
Totient ← +´1=⊢GCD¨1+↕
- Usage:
Totient¨1+↕25
⟨ 1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22 8 20 ⟩
"Number"‿"Totient"‿"Prime?"∾˘{𝕩∾(⊢≍(𝕩-1)⊸=)Totient¨𝕩}1+↕25
┌─
╵ "Number" 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
"Totient" 1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22 8 20
"Prime?" 0 1 1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0
┘
The primes library from bqn-libs includes a Totient function based on factoring that's much faster for numbers in the hundreds and above.
C
Translation of the second Go example
/*Abhishek Ghosh, 7th December 2018*/
#include<stdio.h>
int totient(int n){
int tot = n,i;
for(i=2;i*i<=n;i+=2){
if(n%i==0){
while(n%i==0)
n/=i;
tot-=tot/i;
}
if(i==2)
i=1;
}
if(n>1)
tot-=tot/n;
return tot;
}
int main()
{
int count = 0,n,tot;
printf(" n %c prime",237);
printf("\n---------------\n");
for(n=1;n<=25;n++){
tot = totient(n);
if(n-1 == tot)
count++;
printf("%2d %2d %s\n", n, tot, n-1 == tot?"True":"False");
}
printf("\nNumber of primes up to %6d =%4d\n", 25,count);
for(n = 26; n <= 100000; n++){
tot = totient(n);
if(tot == n-1)
count++;
if(n == 100 || n == 1000 || n%10000 == 0){
printf("\nNumber of primes up to %6d = %4d\n", n, count);
}
}
return 0;
}
Output :
n φ prime --------------- 1 1 False 2 1 True 3 2 True 4 2 False 5 4 True 6 2 False 7 6 True 8 4 False 9 6 False 10 4 False 11 10 True 12 4 False 13 12 True 14 6 False 15 8 False 16 8 False 17 16 True 18 6 False 19 18 True 20 8 False 21 12 False 22 10 False 23 22 True 24 8 False 25 20 False Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 20000 = 2262 Number of primes up to 30000 = 3245 Number of primes up to 40000 = 4203 Number of primes up to 50000 = 5133 Number of primes up to 60000 = 6057 Number of primes up to 70000 = 6935 Number of primes up to 80000 = 7837 Number of primes up to 90000 = 8713 Number of primes up to 100000 = 9592
C#
using static System.Console;
using static System.Linq.Enumerable;
public class Program
{
static void Main()
{
for (int i = 1; i <= 25; i++) {
int t = Totient(i);
WriteLine(i + "\t" + t + (t == i - 1 ? "\tprime" : ""));
}
WriteLine();
for (int i = 100; i <= 100_000; i *= 10) {
WriteLine($"{Range(1, i).Count(x => Totient(x) + 1 == x):n0} primes below {i:n0}");
}
}
static int Totient(int n) {
if (n < 3) return 1;
if (n == 3) return 2;
int totient = n;
if ((n & 1) == 0) {
totient >>= 1;
while (((n >>= 1) & 1) == 0) ;
}
for (int i = 3; i * i <= n; i += 2) {
if (n % i == 0) {
totient -= totient / i;
while ((n /= i) % i == 0) ;
}
}
if (n > 1) totient -= totient / n;
return totient;
}
}
- Output:
1 1 2 1 prime 3 2 prime 4 2 5 4 prime 6 2 7 6 prime 8 4 9 6 10 4 11 10 prime 12 4 13 12 prime 14 6 15 8 16 8 17 16 prime 18 6 19 18 prime 20 8 21 12 22 10 23 22 prime 24 8 25 20 25 primes below 100 168 primes below 1,000 1,229 primes below 10,000 9,592 primes below 100,000
C++
#include <cassert>
#include <iomanip>
#include <iostream>
#include <vector>
class totient_calculator {
public:
explicit totient_calculator(int max) : totient_(max + 1) {
for (int i = 1; i <= max; ++i)
totient_[i] = i;
for (int i = 2; i <= max; ++i) {
if (totient_[i] < i)
continue;
for (int j = i; j <= max; j += i)
totient_[j] -= totient_[j] / i;
}
}
int totient(int n) const {
assert (n >= 1 && n < totient_.size());
return totient_[n];
}
bool is_prime(int n) const {
return totient(n) == n - 1;
}
private:
std::vector<int> totient_;
};
int count_primes(const totient_calculator& tc, int min, int max) {
int count = 0;
for (int i = min; i <= max; ++i) {
if (tc.is_prime(i))
++count;
}
return count;
}
int main() {
const int max = 10000000;
totient_calculator tc(max);
std::cout << " n totient prime?\n";
for (int i = 1; i <= 25; ++i) {
std::cout << std::setw(2) << i
<< std::setw(9) << tc.totient(i)
<< std::setw(8) << (tc.is_prime(i) ? "yes" : "no") << '\n';
}
for (int n = 100; n <= max; n *= 10) {
std::cout << "Count of primes up to " << n << ": "
<< count_primes(tc, 1, n) << '\n';
}
return 0;
}
- Output:
n totient prime? 1 1 no 2 1 yes 3 2 yes 4 2 no 5 4 yes 6 2 no 7 6 yes 8 4 no 9 6 no 10 4 no 11 10 yes 12 4 no 13 12 yes 14 6 no 15 8 no 16 8 no 17 16 yes 18 6 no 19 18 yes 20 8 no 21 12 no 22 10 no 23 22 yes 24 8 no 25 20 no Count of primes up to 100: 25 Count of primes up to 1000: 168 Count of primes up to 10000: 1229 Count of primes up to 100000: 9592 Count of primes up to 1000000: 78498 Count of primes up to 10000000: 664579
CLU
totient = proc (n: int) returns (int)
tot: int := n
i: int := 2
while i*i <= n do
if n//i = 0 then
while n//i = 0 do
n := n/i
end
tot := tot-tot/i
end
if i=2 then i:=1 end
i := i+2
end
if n>1 then
tot := tot-tot/n
end
return(tot)
end totient
start_up = proc ()
po: stream := stream$primary_output()
count: int := 0
stream$putl(po, " N Totient Prime")
for n: int in int$from_to(1, 25) do
tot: int := totient(n)
stream$putright(po, int$unparse(n), 2)
stream$putright(po, int$unparse(tot), 9)
if n-1 = tot then
stream$putright(po, "Yes", 7)
count := count + 1
else
stream$putright(po, "No", 7)
end
stream$putl(po, "")
end
stream$putl(po, "Number of primes up to 25:\t" || int$unparse(count))
for n: int in int$from_to(26, 100000) do
if totient(n) = n-1 then
count := count + 1
end
if n = 100 cor n = 1000 cor n // 10000 = 0 then
stream$putl(po, "Number of primes up to "
|| int$unparse(n) || ":\t"
|| int$unparse(count))
end
end
end start_up- Output:
N Totient Prime 1 1 No 2 1 Yes 3 2 Yes 4 2 No 5 4 Yes 6 2 No 7 6 Yes 8 4 No 9 6 No 10 4 No 11 10 Yes 12 4 No 13 12 Yes 14 6 No 15 8 No 16 8 No 17 16 Yes 18 6 No 19 18 Yes 20 8 No 21 12 No 22 10 No 23 22 Yes 24 8 No 25 20 No Number of primes up to 25: 9 Number of primes up to 100: 25 Number of primes up to 1000: 168 Number of primes up to 10000: 1229 Number of primes up to 20000: 2262 Number of primes up to 30000: 3245 Number of primes up to 40000: 4203 Number of primes up to 50000: 5133 Number of primes up to 60000: 6057 Number of primes up to 70000: 6935 Number of primes up to 80000: 7837 Number of primes up to 90000: 8713 Number of primes up to 100000: 9592
COBOL
IDENTIFICATION DIVISION.
PROGRAM-ID. TOTIENT-FUNCTION.
DATA DIVISION.
WORKING-STORAGE SECTION.
01 TOTIENT-CALCULATION.
03 N PIC 9(6).
03 TOTIENT PIC 9(6).
03 DIVISOR PIC 9(6).
03 DIV-SQUARE PIC 9(6).
03 MODULUS PIC 9(6).
01 LOOP-VARS.
03 PRIME-COUNT PIC 9(4).
03 CURRENT-N PIC 9(6).
03 PRIME-TOTIENT PIC 9(6).
01 OUTPUTFORMAT.
03 HEADER PIC X(20) VALUE " N Totient Prime?".
03 TOTIENT-ROW.
05 OUT-N PIC Z9.
05 FILLER PIC XX VALUE SPACES.
05 OUT-TOTIENT PIC Z(6)9.
05 FILLER PIC XX VALUE SPACES.
05 OUT-PRIME PIC X(5).
03 PRIME-COUNT-ROW.
05 FILLER PIC X(22) VALUE "Number of primes up to".
05 OUT-PRMTO PIC Z(5)9.
05 FILLER PIC XX VALUE ": ".
05 OUT-PRMCNT PIC ZZZ9.
PROCEDURE DIVISION.
BEGIN.
MOVE ZERO TO PRIME-COUNT.
DISPLAY HEADER.
PERFORM SHOW-SMALL-TOTIENT
VARYING CURRENT-N FROM 1 BY 1
UNTIL CURRENT-N IS GREATER THAN 25.
MOVE 25 TO OUT-PRMTO.
MOVE PRIME-COUNT TO OUT-PRMCNT.
DISPLAY PRIME-COUNT-ROW.
PERFORM COUNT-PRIMES
VARYING CURRENT-N FROM 26 BY 1
UNTIL CURRENT-N IS GREATER THAN 10000.
STOP RUN.
SHOW-SMALL-TOTIENT.
MOVE CURRENT-N TO N.
PERFORM CALCULATE-TOTIENT.
MOVE CURRENT-N TO OUT-N.
MOVE TOTIENT TO OUT-TOTIENT.
MOVE "No" TO OUT-PRIME.
SUBTRACT 1 FROM CURRENT-N GIVING PRIME-TOTIENT.
IF TOTIENT IS EQUAL TO PRIME-TOTIENT,
MOVE "Yes" TO OUT-PRIME,
ADD 1 TO PRIME-COUNT.
DISPLAY TOTIENT-ROW.
COUNT-PRIMES.
MOVE CURRENT-N TO N.
PERFORM CALCULATE-TOTIENT.
SUBTRACT 1 FROM CURRENT-N GIVING PRIME-TOTIENT.
IF TOTIENT IS EQUAL TO PRIME-TOTIENT, ADD 1 TO PRIME-COUNT.
IF CURRENT-N IS EQUAL TO 100 OR 1000 OR 10000,
MOVE CURRENT-N TO OUT-PRMTO,
MOVE PRIME-COUNT TO OUT-PRMCNT,
DISPLAY PRIME-COUNT-ROW.
CALCULATE-TOTIENT.
MOVE N TO TOTIENT.
MOVE 2 TO DIVISOR.
MOVE 4 TO DIV-SQUARE.
PERFORM DIVIDE-STEP.
PERFORM DIVIDE-STEP
VARYING DIVISOR FROM 3 BY 2
UNTIL DIV-SQUARE IS GREATER THAN N.
IF N IS GREATER THAN 1,
COMPUTE TOTIENT = TOTIENT - TOTIENT / N.
DIVIDE-STEP.
MULTIPLY DIVISOR BY DIVISOR GIVING DIV-SQUARE.
DIVIDE N BY DIVISOR GIVING MODULUS.
MULTIPLY DIVISOR BY MODULUS.
SUBTRACT MODULUS FROM N GIVING MODULUS.
IF MODULUS IS ZERO,
PERFORM DIVIDE-OUT UNTIL MODULUS IS NOT ZERO,
COMPUTE TOTIENT = TOTIENT - TOTIENT / DIVISOR.
DIVIDE-OUT.
DIVIDE DIVISOR INTO N.
DIVIDE N BY DIVISOR GIVING MODULUS.
MULTIPLY DIVISOR BY MODULUS.
SUBTRACT MODULUS FROM N GIVING MODULUS.
- Output:
N Totient Prime? 1 1 No 2 1 Yes 3 2 Yes 4 2 No 5 4 Yes 6 2 No 7 6 Yes 8 4 No 9 6 No 10 4 No 11 10 Yes 12 4 No 13 12 Yes 14 6 No 15 8 No 16 8 No 17 16 Yes 18 6 No 19 18 Yes 20 8 No 21 12 No 22 10 No 23 22 Yes 24 8 No 25 20 No Number of primes up to 25: 9 Number of primes up to 100: 25 Number of primes up to 1000: 168 Number of primes up to 10000: 1229
Cowgol
include "cowgol.coh";
sub totient(n: uint32): (tot: uint32) is
tot := n;
var i: uint32 := 2;
while i*i <= n loop
if n%i == 0 then
while n%i == 0 loop
n := n/i;
end loop;
tot := tot - tot/i;
end if;
if i == 2 then
i := 1;
end if;
i := i + 2;
end loop;
if n > 1 then
tot := tot - tot/n;
end if;
end sub;
var count: uint16 := 0;
print("N\tTotient\tPrime\n");
var n: uint32 := 1;
while n <= 25 loop
var tot := totient(n);
print_i32(n);
print_char('\t');
print_i32(tot);
print_char('\t');
if n-1 == tot then
count := count + 1;
print("Yes\n");
else
print("No\n");
end if;
n := n + 1;
end loop;
print("Number of primes up to 25:\t");
print_i16(count);
print_nl();
while n <= 100000 loop
tot := totient(n);
if n-1 == tot then
count := count + 1;
end if;
if n == 100 or n == 1000 or n % 10000 == 0 then
print("Number of primes up to ");
print_i32(n);
print(":\t");
print_i16(count);
print_nl();
end if;
n := n + 1;
end loop;- Output:
N Totient Prime 1 1 No 2 1 Yes 3 2 Yes 4 2 No 5 4 Yes 6 2 No 7 6 Yes 8 4 No 9 6 No 10 4 No 11 10 Yes 12 4 No 13 12 Yes 14 6 No 15 8 No 16 8 No 17 16 Yes 18 6 No 19 18 Yes 20 8 No 21 12 No 22 10 No 23 22 Yes 24 8 No 25 20 No Number of primes up to 25: 9 Number of primes up to 100: 25 Number of primes up to 1000: 168 Number of primes up to 10000: 1229 Number of primes up to 20000: 2262 Number of primes up to 30000: 3245 Number of primes up to 40000: 4203 Number of primes up to 50000: 5133 Number of primes up to 60000: 6057 Number of primes up to 70000: 6935 Number of primes up to 80000: 7837 Number of primes up to 90000: 8713 Number of primes up to 100000: 9592
D
import std.stdio;
int totient(int n) {
int tot = n;
for (int i = 2; i * i <= n; i += 2) {
if (n % i == 0) {
while (n % i == 0) {
n /= i;
}
tot -= tot / i;
}
if (i==2) {
i = 1;
}
}
if (n > 1) {
tot -= tot / n;
}
return tot;
}
void main() {
writeln(" n φ prime");
writeln("---------------");
int count = 0;
for (int n = 1; n <= 25; n++) {
int tot = totient(n);
if (n - 1 == tot) {
count++;
}
writefln("%2d %2d %s", n,tot, n - 1 == tot);
}
writeln;
writefln("Number of primes up to %6d = %4d", 25, count);
for (int n = 26; n <= 100_000; n++) {
int tot = totient(n);
if (n - 1 == tot) {
count++;
}
if (n == 100 || n == 1_000 || n % 10_000 == 0) {
writefln("Number of primes up to %6d = %4d", n, count);
}
}
}
- Output:
n φ prime --------------- 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 20000 = 2262 Number of primes up to 30000 = 3245 Number of primes up to 40000 = 4203 Number of primes up to 50000 = 5133 Number of primes up to 60000 = 6057 Number of primes up to 70000 = 6935 Number of primes up to 80000 = 7837 Number of primes up to 90000 = 8713 Number of primes up to 100000 = 9592
Delphi
See Pascal.
Draco
proc totient(word n) word:
word tot, i;
tot := n;
i := 2;
while i*i <= n do
if n%i = 0 then
while n%i = 0 do n := n/i od;
tot := tot - tot/i
fi;
if i=2 then i:=1 fi;
i := i+2
od;
if n>1 then
tot - tot/n
else
tot
fi
corp
proc main() void:
word count, n, tot;
bool prime;
count := 0;
writeln(" N Totient Prime");
for n from 1 upto 25 do
tot := totient(n);
prime := n-1 = tot;
if prime then count := count+1 fi;
writeln(n:2, " ", tot:7, " ", if prime then " Yes" else " No" fi)
od;
writeln("Number of primes up to ",25:6,": ",count:4);
for n from 25 upto 10000 do
if totient(n) = n-1 then count := count+1 fi;
if n=100 or n=1000 or n=10000 then
writeln("Number of primes up to ",n:6,": ",count:4)
fi
od
corp- Output:
N Totient Prime 1 1 No 2 1 Yes 3 2 Yes 4 2 No 5 4 Yes 6 2 No 7 6 Yes 8 4 No 9 6 No 10 4 No 11 10 Yes 12 4 No 13 12 Yes 14 6 No 15 8 No 16 8 No 17 16 Yes 18 6 No 19 18 Yes 20 8 No 21 12 No 22 10 No 23 22 Yes 24 8 No 25 20 No Number of primes up to 25: 9 Number of primes up to 100: 25 Number of primes up to 1000: 168 Number of primes up to 10000: 1229
Dyalect
func totient(n) {
var tot = n
var 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
}
print("n\tphi\tprime")
var count = 0
for n in 1..25 {
var tot = totient(n)
var isPrime = n - 1 == tot
if isPrime {
count += 1
}
print("\(n)\t\(tot)\t\(isPrime)")
}
print("\nNumber of primes up to 25 \t= \(count)")
for n in 26..100000 {
var tot = totient(n)
if tot == n - 1 {
count += 1
}
if n == 100 || n == 1000 || n % 10000 == 0 {
print("Number of primes up to \(n) \t= \(count)")
}
}- Output:
n phi prime 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 20000 = 2262 Number of primes up to 30000 = 3245 Number of primes up to 40000 = 4203 Number of primes up to 50000 = 5133 Number of primes up to 60000 = 6057 Number of primes up to 70000 = 6935 Number of primes up to 80000 = 7837 Number of primes up to 90000 = 8713 Number of primes up to 100000 = 9592
EasyLang
func totient n .
tot = n
i = 2
while i <= sqrt n
if n mod i = 0
while n mod i = 0
n = n div i
.
tot -= tot div i
.
if i = 2
i = 1
.
i += 2
.
if n > 1
tot -= tot div n
.
return tot
.
numfmt 0 3
print " N Prim Phi"
for n = 1 to 25
tot = totient n
x$ = " "
if n - 1 = tot
x$ = " x "
.
print n & x$ & tot
.
print ""
for n = 1 to 100000
tot = totient n
if n - 1 = tot
cnt += 1
.
if n = 100 or n = 1000 or n = 10000 or n = 100000
print n & " - " & cnt & " primes"
.
.
Factor
USING: combinators formatting io kernel math math.primes.factors
math.ranges sequences ;
IN: rosetta-code.totient-function
: Φ ( n -- m )
{
{ [ dup 1 < ] [ drop 0 ] }
{ [ dup 1 = ] [ drop 1 ] }
[
dup unique-factors
[ 1 [ 1 - * ] reduce ] [ product ] bi / *
]
} cond ;
: show-info ( n -- )
[ Φ ] [ swap 2dup - 1 = ] bi ", prime" "" ?
"Φ(%2d) = %2d%s\n" printf ;
: totient-demo ( -- )
25 [1,b] [ show-info ] each nl 0 100,000 [1,b] [
[ dup Φ - 1 = [ 1 + ] when ]
[ dup { 100 1,000 10,000 100,000 } member? ] bi
[ dupd "%4d primes <= %d\n" printf ] [ drop ] if
] each drop ;
MAIN: totient-demo
- Output:
Φ( 1) = 1 Φ( 2) = 1, prime Φ( 3) = 2, prime Φ( 4) = 2 Φ( 5) = 4, prime Φ( 6) = 2 Φ( 7) = 6, prime Φ( 8) = 4 Φ( 9) = 6 Φ(10) = 4 Φ(11) = 10, prime Φ(12) = 4 Φ(13) = 12, prime Φ(14) = 6 Φ(15) = 8 Φ(16) = 8 Φ(17) = 16, prime Φ(18) = 6 Φ(19) = 18, prime Φ(20) = 8 Φ(21) = 12 Φ(22) = 10 Φ(23) = 22, prime Φ(24) = 8 Φ(25) = 20 25 primes <= 100 168 primes <= 1000 1229 primes <= 10000 9592 primes <= 100000
FreeBASIC
#define esPar(n) (((n) And 1) = 0)
#define esImpar(n) (esPar(n) = 0)
Function Totient(n As Integer) As Integer
'delta son números no divisibles por 2,3,5
Dim delta(7) As Integer = {6,4,2,4,2,4,6,2}
Dim As Integer i, quot, idx, result
' div mod por constante es rápido.
'i = 2
result = n
If (2*2 <= n) Then
If Not(esImpar(n)) Then
' eliminar números con factor 2,4,8,16,...
While Not(esImpar(n))
n = n \ 2
Wend
'eliminar los múltiplos de 2
result -= result \ 2
End If
End If
'i = 3
If (3*3 <= n) And (n Mod 3 = 0) Then
Do
quot = n \ 3
If n <> quot*3 Then
Exit Do
Else
n = quot
End If
Loop Until false
result -= result \ 3
End If
'i = 5
If (5*5 <= n) And (n Mod 5 = 0) Then
Do
quot = n \ 5
If n <> quot*5 Then
Exit Do
Else
n = quot
End If
Loop Until false
result -= result \ 5
End If
i = 7
idx = 1
'i = 7,11,13,17,19,23,29,...,49 ..
While i*i <= n
quot = n \ i
If n = quot*i Then
Do
If n <> quot*i Then
Exit Do
Else
n = quot
End If
quot = n \ i
Loop Until false
result -= result \ i
End If
i = i + delta(idx)
idx = (idx+1) And 7
Wend
If n > 1 Then result -= result \ n
Totient = result
End Function
Sub ContandoPrimos(n As Integer)
Dim As Integer i, cnt = 0
For i = 1 To n
If Totient(i) = (i-1) Then cnt += 1
Next i
Print Using " ####### ######"; i-1; cnt
End Sub
Function esPrimo(n As Ulongint) As String
esPrimo = "False"
If n = 1 then Return "False"
If (n=2) Or (n=3) Then Return "True"
If n Mod 2=0 Then Exit Function
If n Mod 3=0 Then Exit Function
Dim As Ulongint limite = Sqr(N)+1
For i As Ulongint = 6 To limite Step 6
If N Mod (i-1)=0 Then Exit Function
If N Mod (i+1)=0 Then Exit Function
Next i
Return "True"
End Function
Sub display(n As Integer)
Dim As Integer idx, phi
If n = 0 Then Exit Sub
Print " n phi(n) esPrimo"
For idx = 1 To n
phi = Totient(idx)
Print Using "### ### \ \"; idx; phi; esPrimo(idx)
Next idx
End Sub
Dim l As Integer
display(25)
Print Chr(10) & " Limite Son primos"
ContandoPrimos(25)
l = 100
Do
ContandoPrimos(l)
l = l*10
Loop Until l > 1000000
End
- Output:
n phi(n) esPrimo
1 1 False
2 1 True
3 2 True
4 2 False
5 4 True
6 2 False
7 6 True
8 4 False
9 6 False
10 4 False
11 10 True
12 4 False
13 12 True
14 6 False
15 8 False
16 8 False
17 16 True
18 6 False
19 18 True
20 8 False
21 12 False
22 10 False
23 22 True
24 8 False
25 20 False
Limite Son primos
25 9
100 25
1000 168
10000 1229
100000 9592
1000000 78498
Fōrmulæ
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for storage and transfer purposes more than visualization and edition.
Programs in Fōrmulæ are created/edited online in its website.
In this page you can see and run the program(s) related to this task and their results. You can also change either the programs or the parameters they are called with, for experimentation, but remember that these programs were created with the main purpose of showing a clear solution of the task, and they generally lack any kind of validation.
Solution
Case 1
Case 2
Forth
Translation from C
: totient \ n -- n' ;
DUP DUP 2 ?DO ( tot n )
DUP I DUP * < IF LEAVE THEN \ for(i=2;i*i<=n;i+=2){
DUP I MOD 0= IF \ if(n%i==0){
BEGIN DUP I /MOD SWAP 0= WHILE ( tot n n/i ) \ while(n%i==0);
NIP ( tot n/i ) \ n/=i;
REPEAT
DROP ( tot n ) \ Remove the new n on exit from loop
SWAP DUP I / - SWAP ( tot' n ) \ tot-=tot/i;
THEN
2 I 2 = + +LOOP \ If I = 2 add 1 else add 2 to loop index.
DUP 1 > IF OVER SWAP / - ELSE DROP THEN ;
: bool. \ f -- ;
IF ." True " ELSE ." False" THEN ;
: count-primes \ n -- n' ;
0 SWAP 2 ?DO I DUP totient 1+ = - LOOP ;
: challenge \ -- ;
CR ." n φ prime" CR
26 1 DO
I 3 .r
I totient DUP 4 .r 4 SPACES
1+ I = bool. CR
LOOP CR
100001 100 DO
." Number of primes up to " I 6 .R ." is " I count-primes 4 .R CR
I 9 * +LOOP ;
- Output:
challenge n φ prime 1 1 False 2 1 True 3 2 True 4 2 False 5 4 True 6 2 False 7 6 True 8 4 False 9 6 False 10 4 False 11 10 True 12 4 False 13 12 True 14 6 False 15 8 False 16 8 False 17 16 True 18 6 False 19 18 True 20 8 False 21 12 False 22 10 False 23 22 True 24 8 False 25 20 False Number of primes up to 100 is 25 Number of primes up to 1000 is 168 Number of primes up to 10000 is 1229 Number of primes up to 100000 is 9592
Go
Results for the larger values of n are very slow to emerge.
package main
import "fmt"
func gcd(n, k int) int {
if n < k || k < 1 {
panic("Need n >= k and k >= 1")
}
s := 1
for n&1 == 0 && k&1 == 0 {
n >>= 1
k >>= 1
s <<= 1
}
t := n
if n&1 != 0 {
t = -k
}
for t != 0 {
for t&1 == 0 {
t >>= 1
}
if t > 0 {
n = t
} else {
k = -t
}
t = n - k
}
return n * s
}
func totient(n int) int {
tot := 0
for k := 1; k <= n; k++ {
if gcd(n, k) == 1 {
tot++
}
}
return tot
}
func main() {
fmt.Println(" n phi prime")
fmt.Println("---------------")
count := 0
for n := 1; n <= 25; n++ {
tot := totient(n)
isPrime := n-1 == tot
if isPrime {
count++
}
fmt.Printf("%2d %2d %t\n", n, tot, isPrime)
}
fmt.Println("\nNumber of primes up to 25 =", count)
for n := 26; n <= 100000; n++ {
tot := totient(n)
if tot == n-1 {
count++
}
if n == 100 || n == 1000 || n%10000 == 0 {
fmt.Printf("\nNumber of primes up to %-6d = %d\n", n, count)
}
}
}
- Output:
n phi prime --------------- 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 20000 = 2262 Number of primes up to 30000 = 3245 Number of primes up to 40000 = 4203 Number of primes up to 50000 = 5133 Number of primes up to 60000 = 6057 Number of primes up to 70000 = 6935 Number of primes up to 80000 = 7837 Number of primes up to 90000 = 8713 Number of primes up to 100000 = 9592
The following much quicker version (runs in less than 150 ms on my machine) uses Euler's product formula rather than repeated invocation of the gcd function to calculate the totient:
package main
import "fmt"
func totient(n int) int {
tot := n
for i := 2; i*i <= n; i += 2 {
if n%i == 0 {
for n%i == 0 {
n /= i
}
tot -= tot / i
}
if i == 2 {
i = 1
}
}
if n > 1 {
tot -= tot / n
}
return tot
}
func main() {
fmt.Println(" n phi prime")
fmt.Println("---------------")
count := 0
for n := 1; n <= 25; n++ {
tot := totient(n)
isPrime := n-1 == tot
if isPrime {
count++
}
fmt.Printf("%2d %2d %t\n", n, tot, isPrime)
}
fmt.Println("\nNumber of primes up to 25 =", count)
for n := 26; n <= 100000; n++ {
tot := totient(n)
if tot == n-1 {
count++
}
if n == 100 || n == 1000 || n%10000 == 0 {
fmt.Printf("\nNumber of primes up to %-6d = %d\n", n, count)
}
}
}
The output is the same as before.
Haskell
{-# LANGUAGE BangPatterns #-}
import Control.Monad (when)
import Data.Bool (bool)
totient
:: (Integral a)
=> a -> a
totient n
| n == 0 = 1 -- by definition phi(0) = 1
| n < 0 = totient (-n) -- phi(-n) is taken to be equal to phi(n)
| otherwise = loop n n 2 --
where
loop !m !tot !i
| i * i > m = bool tot (tot - (tot `div` m)) (1 < m)
| m `mod` i == 0 = loop m_ tot_ i_
| otherwise = loop m tot i_
where
i_
| i == 2 = 3
| otherwise = 2 + i
m_ = nextM m
tot_ = tot - (tot `div` i)
nextM !x
| x `mod` i == 0 = nextM $ x `div` i
| otherwise = x
main :: IO ()
main = do
putStrLn "n\tphi\tprime\n---------------------"
let loop !i !count
| i >= 10 ^ 6 = return ()
| otherwise = do
let i_ = succ i
tot = totient i_
isPrime = tot == pred i_
count_
| isPrime = succ count
| otherwise = count
when (25 >= i_) $
putStrLn $ show i_ ++ "\t" ++ show tot ++ "\t" ++ show isPrime
when
(i_ `elem`
25 :
[ 10 ^ k
| k <- [2 .. 6] ]) $
putStrLn $ "Number of primes up to " ++ show i_ ++ " = " ++ show count_
loop (i + 1) count_
loop 0 0
- Output:
n phi prime --------------------- 1 1 False 2 1 True 3 2 True 4 2 False 5 4 True 6 2 False 7 6 True 8 4 False 9 6 False 10 4 False 11 10 True 12 4 False 13 12 True 14 6 False 15 8 False 16 8 False 17 16 True 18 6 False 19 18 True 20 8 False 21 12 False 22 10 False 23 22 True 24 8 False 25 20 False Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 100000 = 9592 Number of primes up to 1000000 = 78498
J
nth_prime =: p: NB. 2 is the zeroth prime
totient =: 5&p:
primeQ =: 1&p:
NB. first row contains the integer
NB. second row totient
NB. third 1 iff prime
(, totient ,: primeQ) >: i. 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22 8 20
0 1 1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0
NB. primes first exceeding the limits
[&.:(p:inv) 10 ^ 2 + i. 4
101 1009 10007 100003
p:inv 101 1009 10007 100003
25 168 1229 9592
NB. limit and prime count
(,. p:inv) 10 ^ 2 + i. 5
100 25
1000 168
10000 1229
100000 9592
1e6 78498
Java
Efficiently pre-compute all needed values of phi up to 10,000,000 in .344 sec.
public class TotientFunction {
public static void main(String[] args) {
computePhi();
System.out.println("Compute and display phi for the first 25 integers.");
System.out.printf("n Phi IsPrime%n");
for ( int n = 1 ; n <= 25 ; n++ ) {
System.out.printf("%2d %2d %b%n", n, phi[n], (phi[n] == n-1));
}
for ( int i = 2 ; i < 8 ; i++ ) {
int max = (int) Math.pow(10, i);
System.out.printf("The count of the primes up to %,10d = %d%n", max, countPrimes(1, max));
}
}
private static int countPrimes(int min, int max) {
int count = 0;
for ( int i = min ; i <= max ; i++ ) {
if ( phi[i] == i-1 ) {
count++;
}
}
return count;
}
private static final int max = 10000000;
private static final int[] phi = new int[max+1];
private static final void computePhi() {
for ( int i = 1 ; i <= max ; i++ ) {
phi[i] = i;
}
for ( int i = 2 ; i <= max ; i++ ) {
if (phi[i] < i) continue;
for ( int j = i ; j <= max ; j += i ) {
phi[j] -= phi[j] / i;
}
}
}
}
- Output:
Compute and display phi for the first 25 integers. n Phi IsPrime 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false The count of the primes up to 100 = 25 The count of the primes up to 1,000 = 168 The count of the primes up to 10,000 = 1229 The count of the primes up to 100,000 = 9592 The count of the primes up to 1,000,000 = 78498 The count of the primes up to 10,000,000 = 664579
jq
# jq optimizes the recursive call of _gcd in the following:
def gcd(a;b):
def _gcd:
if .[1] != 0 then [.[1], .[0] % .[1]] | _gcd else .[0] end;
[a,b] | _gcd ;
def count(s): reduce s as $x (0; .+1);
def totient:
. as $n
| count( range(0; .) | select( gcd($n; .) == 1) );
# input: determines the maximum via range(0; .)
# and so may be `infinite`
def primes_via_totient:
range(0; .) | select(totient == . - 1);The tasks
def task:
def task1($n):
range(1;$n)
| totient as $totient
| {i: ., $totient, isprime: ($totient == ( . - 1 ))};
task1(26);
def onepass:
reduce (10000 | primes_via_totient) as $p ({};
if $p < 10000
then .["10^4"] += 1
| if $p < 1000
then .["10^3"] += 1
| if $p < 100
then .["10^2"] += 1
else . end else . end else . end) ;
task, "\nCounts of primes up to the given limits:", onepass- Output:
{"i":1,"totient":1,"isprime":false}
{"i":2,"totient":1,"isprime":true}
{"i":3,"totient":2,"isprime":true}
{"i":4,"totient":2,"isprime":false}
{"i":5,"totient":4,"isprime":true}
{"i":6,"totient":2,"isprime":false}
{"i":7,"totient":6,"isprime":true}
{"i":8,"totient":4,"isprime":false}
{"i":9,"totient":6,"isprime":false}
{"i":10,"totient":4,"isprime":false}
{"i":11,"totient":10,"isprime":true}
{"i":12,"totient":4,"isprime":false}
{"i":13,"totient":12,"isprime":true}
{"i":14,"totient":6,"isprime":false}
{"i":15,"totient":8,"isprime":false}
{"i":16,"totient":8,"isprime":false}
{"i":17,"totient":16,"isprime":true}
{"i":18,"totient":6,"isprime":false}
{"i":19,"totient":18,"isprime":true}
{"i":20,"totient":8,"isprime":false}
{"i":21,"totient":12,"isprime":false}
{"i":22,"totient":10,"isprime":false}
{"i":23,"totient":22,"isprime":true}
{"i":24,"totient":8,"isprime":false}
{"i":25,"totient":20,"isprime":false}
Counts of primes up to the given limits:
{"10^4":1229,"10^3":168,"10^2":25}
Julia
φ(n) = sum(1 for k in 1:n if gcd(n, k) == 1)
is_prime(n) = φ(n) == n - 1
function runphitests()
for n in 1:25
println(" φ($n) == $(φ(n))", is_prime(n) ? ", is prime" : "")
end
count = 0
for n in 1:100_000
count += is_prime(n)
if n in [100, 1000, 10_000, 100_000]
println("Primes up to $n: $count")
end
end
end
runphitests()
- Output:
φ(1) == 1 φ(2) == 1, is prime φ(3) == 2, is prime φ(4) == 2 φ(5) == 4, is prime φ(6) == 2 φ(7) == 6, is prime φ(8) == 4 φ(9) == 6 φ(10) == 4 φ(11) == 10, is prime φ(12) == 4 φ(13) == 12, is prime φ(14) == 6 φ(15) == 8 φ(16) == 8 φ(17) == 16, is prime φ(18) == 6 φ(19) == 18, is prime φ(20) == 8 φ(21) == 12 φ(22) == 10 φ(23) == 22, is prime φ(24) == 8 φ(25) == 20 Primes up to 100: 25 Primes up to 1000: 168 Primes up to 10000: 1229 Primes up to 100000: 9592
Kotlin
// Version 1.3.21
fun totient(n: Int): Int {
var tot = n
var nn = n
var i = 2
while (i * i <= nn) {
if (nn % i == 0) {
while (nn % i == 0) nn /= i
tot -= tot / i
}
if (i == 2) i = 1
i += 2
}
if (nn > 1) tot -= tot / nn
return tot
}
fun main() {
println(" n phi prime")
println("---------------")
var count = 0
for (n in 1..25) {
val tot = totient(n)
val isPrime = n - 1 == tot
if (isPrime) count++
System.out.printf("%2d %2d %b\n", n, tot, isPrime)
}
println("\nNumber of primes up to 25 = $count")
for (n in 26..100_000) {
val tot = totient(n)
if (tot == n-1) count++
if (n == 100 || n == 1000 || n % 10_000 == 0) {
System.out.printf("\nNumber of primes up to %-6d = %d\n", n, count)
}
}
}
- Output:
Same as Go example.
Lua
Averages about 7 seconds under LuaJIT
-- Return the greatest common denominator of x and y
function gcd (x, y)
return y == 0 and math.abs(x) or gcd(y, x % y)
end
-- Return the totient number for n
function totient (n)
local count = 0
for k = 1, n do
if gcd(n, k) == 1 then count = count + 1 end
end
return count
end
-- Determine (inefficiently) whether p is prime
function isPrime (p)
return totient(p) == p - 1
end
-- Output totient and primality for the first 25 integers
print("n", string.char(237), "prime")
print(string.rep("-", 21))
for i = 1, 25 do
print(i, totient(i), isPrime(i))
end
-- Count the primes up to 100, 1000 and 10000
local pCount, i, limit = 0, 1
for power = 2, 4 do
limit = 10 ^ power
repeat
i = i + 1
if isPrime(i) then pCount = pCount + 1 end
until i == limit
print("\nThere are " .. pCount .. " primes below " .. limit)
end
- Output:
n φ prime --------------------- 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false There are 25 primes below 100 There are 168 primes below 1000 There are 1229 primes below 10000
Alternative, faster version - around 0.04 seconds using LuaJIT on TIO.RUN.
do
function totient( n ) -- returns the number of integers k where 1 <= k <= n that are mutually prime to n
if n < 3 then return 1
elseif n == 3 then return 2
else
local result, v, i = n, n, 2
while i * i <= v do
if v % i == 0 then
while v % i == 0 do v = math.floor( v / i ) end
result = result - math.floor( result / i )
end
if i == 2 then
i = 1
end
i = i + 2
end
if v > 1 then result = result - math.floor( result / v ) end
return result
end
end
-- show the totient function values for the first 25 integers
io.write( " n phi(n) remarks\n" )
for n = 1,25 do
local tn = totient( n )
io.write( string.format( "%2d", n ), ": ", string.format( "%5d", tn )
, ( tn == n - 1 and tn ~= 0 and " n is prime" or "" )
, "\n"
)
end
-- use the totient function to count primes
local n100, n1000, n10000, n100000 = 0, 0, 0, 0
for n = 1,100000 do
if totient( n ) == n - 1 then
if n <= 100 then n100 = n100 + 1 end
if n <= 1000 then n1000 = n1000 + 1 end
if n <= 10000 then n10000 = n10000 + 1 end
if n <= 100000 then n100000 = n100000 + 1 end
end
end
io.write( "There are ", string.format( "%6d", n100 ), " primes below 100\n" )
io.write( "There are ", string.format( "%6d", n1000 ), " primes below 1 000\n" )
io.write( "There are ", string.format( "%6d", n10000 ), " primes below 10 000\n" )
io.write( "There are ", string.format( "%6d", n100000 ), " primes below 100 000\n" )
end
- Output:
n phi(n) remarks 1: 1 2: 1 n is prime 3: 2 n is prime 4: 2 5: 4 n is prime 6: 2 7: 6 n is prime 8: 4 9: 6 10: 4 11: 10 n is prime 12: 4 13: 12 n is prime 14: 6 15: 8 16: 8 17: 16 n is prime 18: 6 19: 18 n is prime 20: 8 21: 12 22: 10 23: 22 n is prime 24: 8 25: 20 There are 25 primes below 100 There are 168 primes below 1 000 There are 1229 primes below 10 000 There are 9592 primes below 100 000
MAD
NORMAL MODE IS INTEGER
BOOLEAN PRM
INTERNAL FUNCTION(A,B)
ENTRY TO REM.
FUNCTION RETURN A-A/B*B
END OF FUNCTION
INTERNAL FUNCTION(NN)
ENTRY TO TOTENT.
N = NN
TOT = N
THROUGH STEP, FOR I=2, 2, I*I.G.N
WHENEVER REM.(N,I).E.0
THROUGH DIV, FOR N=N, 0, REM.(N,I).NE.0
DIV N = N/I
TOT = TOT-TOT/I
END OF CONDITIONAL
WHENEVER I.E.2, I=1
STEP CONTINUE
WHENEVER N.G.1, TOT = TOT-TOT/N
FUNCTION RETURN TOT
END OF FUNCTION
COUNT = 0
PRINT FORMAT HEADER
THROUGH FRST25, FOR KN=1, 1, KN.G.25
KTOT = TOTENT.(KN)
PRM = KTOT.E.KN-1
WHENEVER PRM, COUNT = COUNT + 1
FRST25 PRINT FORMAT NUMTOT,KN,KTOT,PRM
PRINT FORMAT NUMPRM,25,COUNT
THROUGH CNTPRM, FOR KN=26, 1, KN.G.100000
KTOT = TOTENT.(KN)
WHENEVER KTOT.E.KN-1, COUNT = COUNT+1
WHENEVER KN.E.100 .OR. KN.E.1000 .OR. REM.(KN,10000).E.0,
0 PRINT FORMAT NUMPRM,KN,COUNT
CNTPRM CONTINUE
VECTOR VALUES HEADER = $19H N TOTIENT PRIME*$
VECTOR VALUES NUMTOT = $I2,S2,I7,S2,I5*$
VECTOR VALUES NUMPRM = $22HNUMBER OF PRIMES UP TO,I7,1H:,I6*$
END OF PROGRAM- Output:
N TOTIENT PRIME 1 1 0 2 1 1 3 2 1 4 2 0 5 4 1 6 2 0 7 6 1 8 4 0 9 6 0 10 4 0 11 10 1 12 4 0 13 12 1 14 6 0 15 8 0 16 8 0 17 16 1 18 6 0 19 18 1 20 8 0 21 12 0 22 10 0 23 22 1 24 8 0 25 20 0 NUMBER OF PRIMES UP TO 25: 9 NUMBER OF PRIMES UP TO 100: 25 NUMBER OF PRIMES UP TO 1000: 168 NUMBER OF PRIMES UP TO 10000: 1229 NUMBER OF PRIMES UP TO 20000: 2262 NUMBER OF PRIMES UP TO 30000: 3245 NUMBER OF PRIMES UP TO 40000: 4203 NUMBER OF PRIMES UP TO 50000: 5133 NUMBER OF PRIMES UP TO 60000: 6057 NUMBER OF PRIMES UP TO 70000: 6935 NUMBER OF PRIMES UP TO 80000: 7837 NUMBER OF PRIMES UP TO 90000: 8713 NUMBER OF PRIMES UP TO 100000: 9592
Mathematica / Wolfram Language
Do[
tmp = EulerPhi[i];
If[i - 1 == tmp,
Print["\[CurlyPhi](", i, ")=", tmp, ", is prime"]
,
Print["\[CurlyPhi](", i, ")=", tmp]
]
,
{i, 25}
]
Count[EulerPhi[Range[100]] - Range[100], -1]
Count[EulerPhi[Range[1000]] - Range[1000], -1]
Count[EulerPhi[Range[10000]] - Range[10000], -1]
Count[EulerPhi[Range[100000]] - Range[100000], -1]
(*Alternative much faster way of findings the number primes up to a number*)
(*PrimePi[100]*)
(*PrimePi[1000]*)
(*PrimePi[10000]*)
(*PrimePi[100000]*)
- Output:
φ(1)=1 φ(2)=1, is prime φ(3)=2, is prime φ(4)=2 φ(5)=4, is prime φ(6)=2 φ(7)=6, is prime φ(8)=4 φ(9)=6 φ(10)=4 φ(11)=10, is prime φ(12)=4 φ(13)=12, is prime φ(14)=6 φ(15)=8 φ(16)=8 φ(17)=16, is prime φ(18)=6 φ(19)=18, is prime φ(20)=8 φ(21)=12 φ(22)=10 φ(23)=22, is prime φ(24)=8 φ(25)=20 25 168 1229 9592
Miranda
main :: [sys_message]
main = [Stdout (lay (map showline [1..25])),
Stdout (lay (map countprimes (25:map (10^) [2..5])))]
countprimes :: num->[char]
countprimes n = "There are " ++ show amount ++ " primes up to " ++ show n
where amount = #filter prime [2..n]
showline :: num->[char]
showline n = "phi(" ++ show n ++ ") = " ++ show (totient n) ++ ", " ++ kind
where kind = "prime", if prime n
= "composite", otherwise
prime :: num->bool
prime n = totient n = n - 1
totient :: num->num
totient n = loop n n (2:[3, 5..])
where loop tot n (d:ds)
= tot, if n<=1
= tot - tot div n, if d*d > n
= loop tot n ds, if n mod d ~= 0
= loop (tot - tot div d) (remfac n d) ds, otherwise
remfac n d = n, if n mod d ~= 0
= remfac (n div d) d, otherwise- Output:
phi(1) = 1, composite phi(2) = 1, prime phi(3) = 2, prime phi(4) = 2, composite phi(5) = 4, prime phi(6) = 2, composite phi(7) = 6, prime phi(8) = 4, composite phi(9) = 6, composite phi(10) = 4, composite phi(11) = 10, prime phi(12) = 4, composite phi(13) = 12, prime phi(14) = 6, composite phi(15) = 8, composite phi(16) = 8, composite phi(17) = 16, prime phi(18) = 6, composite phi(19) = 18, prime phi(20) = 8, composite phi(21) = 12, composite phi(22) = 10, composite phi(23) = 22, prime phi(24) = 8, composite phi(25) = 20, composite There are 9 primes up to 25 There are 25 primes up to 100 There are 168 primes up to 1000 There are 1229 primes up to 10000 There are 9592 primes up to 100000
Modula-2
MODULE TotientFunction;
FROM InOut IMPORT WriteString, WriteCard, WriteLn;
VAR count, n, tot: CARDINAL;
PROCEDURE totient(n: CARDINAL): CARDINAL;
VAR tot, i: CARDINAL;
BEGIN
tot := n;
i := 2;
WHILE i*i <= n DO
IF n MOD i = 0 THEN
WHILE n MOD i = 0 DO
n := n DIV i
END;
DEC(tot, tot DIV i)
END;
IF i=2 THEN i := 1 END;
INC(i, 2)
END;
IF n>1 THEN
DEC(tot, tot DIV n)
END;
RETURN tot
END totient;
PROCEDURE ShowPrimeCount(n, count: CARDINAL);
BEGIN
WriteString("Number of primes up to");
WriteCard(n, 6);
WriteString(": ");
WriteCard(count, 4);
WriteLn
END ShowPrimeCount;
BEGIN
count := 0;
WriteString(" N Totient Prime");
WriteLn;
FOR n := 1 TO 25 DO
tot := totient(n);
WriteCard(n, 2);
WriteCard(tot, 9);
IF tot = n-1 THEN
WriteString(" Yes");
INC(count)
ELSE
WriteString(" No")
END;
WriteLn
END;
ShowPrimeCount(25, count);
FOR n := 26 TO 10000 DO
IF totient(n) = n-1 THEN INC(count) END;
IF (n=100) OR (n=1000) OR (n=10000) THEN
ShowPrimeCount(n, count)
END;
END
END TotientFunction.
- Output:
N Totient Prime 1 1 No 2 1 Yes 3 2 Yes 4 2 No 5 4 Yes 6 2 No 7 6 Yes 8 4 No 9 6 No 10 4 No 11 10 Yes 12 4 No 13 12 Yes 14 6 No 15 8 No 16 8 No 17 16 Yes 18 6 No 19 18 Yes 20 8 No 21 12 No 22 10 No 23 22 Yes 24 8 No 25 20 No Number of primes up to 25: 9 Number of primes up to 100: 25 Number of primes up to 1000: 168 Number of primes up to 10000: 1229
Nim
import strformat
func totient(n: int): int =
var tot = n
var nn = n
var i = 2
while i * i <= nn:
if nn mod i == 0:
while nn mod i == 0:
nn = nn div i
dec tot, tot div i
if i == 2:
i = 1
inc i, 2
if nn > 1:
dec tot, tot div nn
tot
echo " n φ prime"
echo "---------------"
var count = 0
for n in 1..25:
let tot = totient(n)
let isPrime = n - 1 == tot
if isPrime:
inc count
echo fmt"{n:2} {tot:2} {isPrime}"
echo ""
echo fmt"Number of primes up to {25:>6} = {count:>4}"
for n in 26..100_000:
let tot = totient(n)
if tot == n - 1:
inc count
if n == 100 or n == 1000 or n mod 10_000 == 0:
echo fmt"Number of primes up to {n:>6} = {count:>4}"
- Output:
n φ prime --------------- 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 20000 = 2262 Number of primes up to 30000 = 3245 Number of primes up to 40000 = 4203 Number of primes up to 50000 = 5133 Number of primes up to 60000 = 6057 Number of primes up to 70000 = 6935 Number of primes up to 80000 = 7837 Number of primes up to 90000 = 8713 Number of primes up to 100000 = 9592
Pascal
Yes, a very slow possibility to check prime
{$IFDEF FPC}
{$MODE DELPHI}
{$IFEND}
function gcd_mod(u, v: NativeUint): NativeUint;inline;
//prerequisites u > v and u,v > 0
var
t: NativeUInt;
begin
repeat
t := u;
u := v;
v := t mod v;
until v = 0;
gcd_mod := u;
end;
function Totient(n:NativeUint):NativeUint;
var
i : NativeUint;
Begin
result := 1;
For i := 2 to n do
inc(result,ORD(GCD_mod(n,i)=1));
end;
function CheckPrimeTotient(n:NativeUint):Boolean;inline;
begin
result := (Totient(n) = (n-1));
end;
procedure OutCountPrimes(n:NativeUInt);
var
i,cnt : NativeUint;
begin
cnt := 0;
For i := 1 to n do
inc(cnt,Ord(CheckPrimeTotient(i)));
writeln(n:10,cnt:8);
end;
procedure display(n:NativeUint);
var
idx,phi : NativeUint;
Begin
if n = 0 then
EXIT;
writeln('number n':5,'Totient(n)':11,'isprime':8);
For idx := 1 to n do
Begin
phi := Totient(idx);
writeln(idx:4,phi:10,(phi=(idx-1)):12);
end
end;
var
i : NativeUint;
Begin
display(25);
writeln('Limit primecount');
i := 100;
repeat
OutCountPrimes(i);
i := i*10;
until i >100000;
end.
- Output
number n Totient(n) isprime
1 1 FALSE
2 1 TRUE
3 2 TRUE
4 2 FALSE
5 4 TRUE
6 2 FALSE
7 6 TRUE
8 4 FALSE
9 6 FALSE
10 4 FALSE
11 10 TRUE
12 4 FALSE
13 12 TRUE
14 6 FALSE
15 8 FALSE
16 8 FALSE
17 16 TRUE
18 6 FALSE
19 18 TRUE
20 8 FALSE
21 12 FALSE
22 10 FALSE
23 22 TRUE
24 8 FALSE
25 20 FALSE
Limit primecount
100 25
1000 168
10000 1229
100000 9592
real 3m39,745s
alternative
changing Totient-funtion in program atop to Computing Euler's totient function on wikipedia, like GO and C.
Impressive speedup.Checking with only primes would be even faster.
function totient(n:NativeUInt):NativeUInt;
const
//delta of numbers not divisible by 2,3,5 (0_1+6->7+4->11 ..+6->29+2->3_1
delta : array[0..7] of NativeUint = (6,4,2,4,2,4,6,2);
var
i, quot,idx: NativeUint;
Begin
// div mod by constant is fast.
//i = 2
result := n;
if (2*2 <= n) then
Begin
IF not(ODD(n)) then
Begin
// remove numbers with factor 2,4,8,16, ...
while not(ODD(n)) do
n := n DIV 2;
//remove count of multiples of 2
dec(result,result DIV 2);
end;
end;
//i = 3
If (3*3 <= n) AND (n mod 3 = 0) then
Begin
repeat
quot := n DIV 3;
IF n <> quot*3 then
BREAK
else
n := quot;
until false;
dec(result,result DIV 3);
end;
//i = 5
If (5*5 <= n) AND (n mod 5 = 0) then
Begin
repeat
quot := n DIV 5;
IF n <> quot*5 then
BREAK
else
n := quot;
until false;
dec(result,result DIV 5);
end;
i := 7;
idx := 1;
//i = 7,11,13,17,19,23,29, ...49 ..
while i*i <= n do
Begin
quot := n DIV i;
if n = quot*i then
Begin
repeat
IF n <> quot*i then
BREAK
else
n := quot;
quot := n DIV i;
until false;
dec(result,result DIV i);
end;
i := i + delta[idx];
idx := (idx+1) AND 7;
end;
if n> 1 then
dec(result,result div n);
end;
- Output
number n Totient(n) isprime
1 1 FALSE
2 1 TRUE
3 2 TRUE
4 2 FALSE
5 4 TRUE
6 2 FALSE
7 6 TRUE
8 4 FALSE
9 6 FALSE
10 4 FALSE
11 10 TRUE
12 4 FALSE
13 12 TRUE
14 6 FALSE
15 8 FALSE
16 8 FALSE
17 16 TRUE
18 6 FALSE
19 18 TRUE
20 8 FALSE
21 12 FALSE
22 10 FALSE
23 22 TRUE
24 8 FALSE
25 20 FALSE
Limit primecount
100 25
1000 168
10000 1229
100000 9592
1000000 78498
10000000 664579
real 0m7,369s
Perl
Direct calculation of 𝜑
use utf8;
binmode STDOUT, ":utf8";
sub gcd {
my ($u, $v) = @_;
while ($v) {
($u, $v) = ($v, $u % $v);
}
return abs($u);
}
push @𝜑, 0;
for $t (1..10000) {
push @𝜑, scalar grep { 1 == gcd($_,$t) } 1..$t;
}
printf "𝜑(%2d) = %3d%s\n", $_, $𝜑[$_], $_ - $𝜑[$_] - 1 ? '' : ' Prime' for 1 .. 25;
print "\n";
for $limit (100, 1000, 10000) {
printf "Count of primes <= $limit: %d\n", scalar grep {$_ == $𝜑[$_] + 1} 0..$limit;
}
- Output:
𝜑( 1) = 1 𝜑( 2) = 1 Prime 𝜑( 3) = 2 Prime 𝜑( 4) = 2 𝜑( 5) = 4 Prime 𝜑( 6) = 2 𝜑( 7) = 6 Prime 𝜑( 8) = 4 𝜑( 9) = 6 𝜑(10) = 4 𝜑(11) = 10 Prime 𝜑(12) = 4 𝜑(13) = 12 Prime 𝜑(14) = 6 𝜑(15) = 8 𝜑(16) = 8 𝜑(17) = 16 Prime 𝜑(18) = 6 𝜑(19) = 18 Prime 𝜑(20) = 8 𝜑(21) = 12 𝜑(22) = 10 𝜑(23) = 22 Prime 𝜑(24) = 8 𝜑(25) = 20 Count of primes <= 100: 25 Count of primes <= 1000: 168 Count of primes <= 10000: 1229
Using 'ntheory' library
Much faster. Output is the same as above.
use utf8;
binmode STDOUT, ":utf8";
use ntheory qw(euler_phi);
my @𝜑 = euler_phi(0,10000); # Returns list of all values in range
printf "𝜑(%2d) = %3d%s\n", $_, $𝜑[$_], $_ - $𝜑[$_] - 1 ? '' : ' Prime' for 1 .. 25;
print "\n";
for $limit (100, 1000, 10000) {
printf "Count of primes <= $limit: %d\n", scalar grep {$_ == $𝜑[$_] + 1} 0..$limit;
}
Phix
As of 1.0.2 there is a builtin phi(), which is memoised and hence much faster.
with javascript_semantics
atom t0 = time()
function totient(integer n)
integer tot = n, i = 2
while i*i<=n do
if mod(n,i)=0 then
while true do
n /= i
if mod(n,i)!=0 then exit end if
end while
tot -= tot/i
end if
i += iff(i=2?1:2)
end while
if n>1 then
tot -= tot/n
end if
return tot
end function
printf(1," n phi prime\n")
printf(1," --------------\n")
for n=1 to 25 do
integer tot = totient(n)
bool prime = (n-1=tot)
printf(1,"%2d %2d %t\n",{n,tot,prime})
end for
printf(1,"\n")
integer count = 0
for n=1 to 1e6 do
-- count += (totient(n)=n-1)
count += (phi(n)=n-1)
if find(n,{25,1e2,1e3,1e4,1e5,1e6,1e7}) then
printf(1,"Number of primes up to %d = %d\n",{n,count})
end if
end for
?elapsed(time()-t0)
- Output:
Using the builtin phi() and with limit upped to 1e7 (above finishes in 0.8s).
n phi prime -------------- 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 100000 = 9592 Number of primes up to 1000000 = 78498 Number of primes up to 10000000 = 664579 "49.3s"
alternate, much faster
with javascript_semantics
atom t0 = time()
constant maxphi = 1e7
sequence s_phi
procedure computePhi()
s_phi = tagset(maxphi)
for i=2 to maxphi do
if s_phi[i]=i then
for j=i to maxphi by i do
s_phi[j] -= s_phi[j] / i;
end for
end if
end for
end procedure
function countPrimes(int lo, hi)
int count = 0;
for i=lo to hi do
if s_phi[i] == i-1 then
count += 1
end if
end for
return count;
end function
computePhi()
printf(1," n phi prime\n")
printf(1," --------------\n")
for n=1 to 25 do
integer tot = s_phi[n]
bool prime = (n-1=tot)
printf(1,"%2d %2d %t\n",{n,tot,prime})
end for
printf(1,"\n")
for n in {25,1e2,1e3,1e4,1e5,1e6,1e7} do
printf(1,"Number of primes up to %d = %d\n",{n,countPrimes(1,n)})
end for
?elapsed(time()-t0)
output as above, except it finishes in 0.8s, which is better than 60x faster. See also builtins/gcd.e in 1.0.6+ for a currently-commented-out get_phi(), which is an extensible version of this, and a planned replacement for the current phi(), once adequately tested.
PL/I
totientFunction: procedure options(main);
totient: procedure(nn) returns(fixed);
declare (nn, n, i, tot) fixed;
n = nn;
tot = n;
do i=2 repeat(i+2) while(i*i <= n);
if mod(n,i) = 0 then do;
do while(mod(n,i) = 0);
n = n / i;
end;
tot = tot - tot / i;
end;
if i=2 then i=1;
end;
if n>1 then tot = tot - tot/n;
return(tot);
end totient;
showPrimeCount: procedure(n, primeCount);
declare (n, primeCount) fixed;
put skip edit('There are', primeCount, ' primes up to', n) (A,F(5),A,F(6));
end showPrimeCount;
declare (n, primeCount, tot) fixed;
do n = 1 to 25;
tot = totient(n);
put edit('phi(', n, ') = ', tot) (A,F(2),A,F(2));
if tot = n-1 then do;
put list('; prime');
primeCount = primeCount + 1;
end;
put skip;
end;
call showPrimeCount(25, primeCount);
do n = 26 to 10000;
if totient(n) = n-1 then primeCount = primeCount + 1;
if n=100 | n=1000 | n=10000 then
call showPrimeCount(n, primeCount);
end;
end totientFunction;- Output:
phi( 1) = 1 phi( 2) = 1 ; prime phi( 3) = 2 ; prime phi( 4) = 2 phi( 5) = 4 ; prime phi( 6) = 2 phi( 7) = 6 ; prime phi( 8) = 4 phi( 9) = 6 phi(10) = 4 phi(11) = 10 ; prime phi(12) = 4 phi(13) = 12 ; prime phi(14) = 6 phi(15) = 8 phi(16) = 8 phi(17) = 16 ; prime phi(18) = 6 phi(19) = 18 ; prime phi(20) = 8 phi(21) = 12 phi(22) = 10 phi(23) = 22 ; prime phi(24) = 8 phi(25) = 20 There are 9 primes up to 25 There are 25 primes up to 100 There are 168 primes up to 1000 There are 1229 primes up to 10000
PL/M
100H:
BDOS: PROCEDURE (F,A); DECLARE F BYTE, A ADDRESS; GO TO 5; END BDOS;
EXIT: PROCEDURE; GO TO 0; END EXIT;
PRINT: PROCEDURE (S); DECLARE S ADDRESS; CALL BDOS(9,S); END PRINT;
PRINT$NUM: PROCEDURE (N);
DECLARE (N, P) ADDRESS, CH BASED P BYTE;
DECLARE S (6) BYTE INITIAL ('.....$');
P = .S(5);
DIGIT:
P = P-1;
CH = '0' + N MOD 10;
IF (N := N/10) > 0 THEN GO TO DIGIT;
CALL PRINT(P);
END PRINT$NUM;
TOTIENT: PROCEDURE (N) ADDRESS;
DECLARE (TOT, N, I) ADDRESS;
TOT = N;
I = 2;
DO WHILE I*I <= N;
IF N MOD I = 0 THEN DO;
DO WHILE (N := N / I) MOD I = 0; END;
TOT = TOT - TOT / I;
END;
IF I=2 THEN I=1;
I = I+2;
END;
IF N > 1 THEN TOT = TOT - TOT / N;
RETURN TOT;
END TOTIENT;
SHOW$PRIME$COUNT: PROCEDURE (N, COUNT);
DECLARE (N, COUNT) ADDRESS;
CALL PRINT(.'THERE ARE $');
CALL PRINT$NUM(COUNT);
CALL PRINT(.' PRIMES UP TO $');
CALL PRINT$NUM(N);
CALL PRINT(.(13,10,'$'));
END SHOW$PRIME$COUNT;
DECLARE (N, TOT) ADDRESS;
DECLARE PRIME$COUNT ADDRESS INITIAL (0);
DO N = 1 TO 25;
CALL PRINT(.'PHI($');
CALL PRINT$NUM(N);
CALL PRINT(.') = $');
CALL PRINT$NUM(TOT := TOTIENT(N));
IF TOT = N-1 THEN DO;
CALL PRINT(.'; PRIME$');
PRIME$COUNT = PRIME$COUNT + 1;
END;
CALL PRINT(.(13,10,'$'));
END;
CALL SHOW$PRIME$COUNT(25, PRIME$COUNT);
DO N = 26 TO 10000;
IF TOTIENT(N) = N-1 THEN
PRIME$COUNT = PRIME$COUNT + 1;
IF N=100 OR N=1000 OR N=10000 THEN
CALL SHOW$PRIME$COUNT(N, PRIME$COUNT);
END;
CALL EXIT;
EOF- Output:
PHI(1) = 1 PHI(2) = 1; PRIME PHI(3) = 2; PRIME PHI(4) = 2 PHI(5) = 4; PRIME PHI(6) = 2 PHI(7) = 6; PRIME PHI(8) = 4 PHI(9) = 6 PHI(10) = 4 PHI(11) = 10; PRIME PHI(12) = 4 PHI(13) = 12; PRIME PHI(14) = 6 PHI(15) = 8 PHI(16) = 8 PHI(17) = 16; PRIME PHI(18) = 6 PHI(19) = 18; PRIME PHI(20) = 8 PHI(21) = 12 PHI(22) = 10 PHI(23) = 22; PRIME PHI(24) = 8 PHI(25) = 20 THERE ARE 9 PRIMES UP TO 25 THERE ARE 25 PRIMES UP TO 100 THERE ARE 168 PRIMES UP TO 1000 THERE ARE 1229 PRIMES UP TO 10000
PicoLisp
(gc 32)
(de gcd (A B)
(until (=0 B)
(let M (% A B)
(setq A B B M) ) )
(abs A) )
(de totient (N)
(let C 0
(for I N
(and (=1 (gcd N I)) (inc 'C)) )
(cons C (= C (dec N))) ) )
(de p? (N)
(let C 0
(for A N
(and
(cdr (totient A))
(inc 'C) ) )
C ) )
(let Fmt (3 7 10)
(tab Fmt "N" "Phi" "Prime?")
(tab Fmt "-" "---" "------")
(for N 25
(tab Fmt
N
(car (setq @ (totient N)))
(cdr @) ) ) )
(println
(mapcar p? (25 100 1000 10000 100000)) )- Output:
N Phi Prime? - --- ------ 1 1 2 1 T 3 2 T 4 2 5 4 T 6 2 7 6 T 8 4 9 6 10 4 11 10 T 12 4 13 12 T 14 6 15 8 16 8 17 16 T 18 6 19 18 T 20 8 21 12 22 10 23 22 T 24 8 25 20 (9 25 168 1229 9592)
Python
from math import gcd
def φ(n):
return sum(1 for k in range(1, n + 1) if gcd(n, k) == 1)
if __name__ == '__main__':
def is_prime(n):
return φ(n) == n - 1
for n in range(1, 26):
print(f" φ({n}) == {φ(n)}{', is prime' if is_prime(n) else ''}")
count = 0
for n in range(1, 10_000 + 1):
count += is_prime(n)
if n in {100, 1000, 10_000}:
print(f"Primes up to {n}: {count}")
- Output:
φ(1) == 1 φ(2) == 1, is prime φ(3) == 2, is prime φ(4) == 2 φ(5) == 4, is prime φ(6) == 2 φ(7) == 6, is prime φ(8) == 4 φ(9) == 6 φ(10) == 4 φ(11) == 10, is prime φ(12) == 4 φ(13) == 12, is prime φ(14) == 6 φ(15) == 8 φ(16) == 8 φ(17) == 16, is prime φ(18) == 6 φ(19) == 18, is prime φ(20) == 8 φ(21) == 12 φ(22) == 10 φ(23) == 22, is prime φ(24) == 8 φ(25) == 20 Primes up to 100: 25 Primes up to 1000: 168 Primes up to 10000: 1229
Quackery
gcd is defined at Greatest common divisor#Quackery.
[ 0 swap dup times
[ i over gcd
1 = rot + swap ]
drop ] is totient ( n --> n )
[ 0 temp put
times
[ i dup 1+ totient
= temp tally ]
temp take ] is primecount ( n --> n )
25 times
[ say "The totient of "
i^ 1+ dup echo
say " is "
dup totient dup echo
say ", so it is "
1+ != if say "not "
say "prime." cr ]
cr
' [ 100 1000 10000 100000 ]
witheach
[ say "There are "
dup primecount echo
say " primes up to " echo
say "." cr ]((Out}}
The totient of 1 is 1, so it is not prime. The totient of 2 is 1, so it is prime. The totient of 3 is 2, so it is prime. The totient of 4 is 2, so it is not prime. The totient of 5 is 4, so it is prime. The totient of 6 is 2, so it is not prime. The totient of 7 is 6, so it is prime. The totient of 8 is 4, so it is not prime. The totient of 9 is 6, so it is not prime. The totient of 10 is 4, so it is not prime. The totient of 11 is 10, so it is prime. The totient of 12 is 4, so it is not prime. The totient of 13 is 12, so it is prime. The totient of 14 is 6, so it is not prime. The totient of 15 is 8, so it is not prime. The totient of 16 is 8, so it is not prime. The totient of 17 is 16, so it is prime. The totient of 18 is 6, so it is not prime. The totient of 19 is 18, so it is prime. The totient of 20 is 8, so it is not prime. The totient of 21 is 12, so it is not prime. The totient of 22 is 10, so it is not prime. The totient of 23 is 22, so it is prime. The totient of 24 is 8, so it is not prime. The totient of 25 is 20, so it is not prime. There are 25 primes up to 100. There are 168 primes up to 1000. There are 1229 primes up to 10000. There are 9592 primes up to 100000.
Racket
#lang racket
(require math/number-theory)
(define (prime*? n) (= (totient n) (sub1 n)))
(for ([n (in-range 1 26)])
(printf "φ(~a) = ~a~a~a\n"
n
(totient n)
(if (prime*? n) " is prime" "")
(if (prime? n) " (confirmed)" "")))
(for/fold ([count 0] #:result (void)) ([n (in-range 1 10001)])
(define new-count (if (prime*? n) (add1 count) count))
(when (member n '(100 1000 10000))
(printf "Primes up to ~a: ~a\n" n new-count))
new-count)
- Output:
φ(1) = 1 φ(2) = 1 is prime (confirmed) φ(3) = 2 is prime (confirmed) φ(4) = 2 φ(5) = 4 is prime (confirmed) φ(6) = 2 φ(7) = 6 is prime (confirmed) φ(8) = 4 φ(9) = 6 φ(10) = 4 φ(11) = 10 is prime (confirmed) φ(12) = 4 φ(13) = 12 is prime (confirmed) φ(14) = 6 φ(15) = 8 φ(16) = 8 φ(17) = 16 is prime (confirmed) φ(18) = 6 φ(19) = 18 is prime (confirmed) φ(20) = 8 φ(21) = 12 φ(22) = 10 φ(23) = 22 is prime (confirmed) φ(24) = 8 φ(25) = 20 Primes up to 100: 25 Primes up to 1000: 168 Primes up to 10000: 1229
Raku
(formerly Perl 6)
This is an incredibly inefficient way of finding prime numbers.
use Prime::Factor;
my \𝜑 = 0, |(1..*).hyper.map: -> \t { t * [*] t.&prime-factors.squish.map: { 1 - 1/$_ } }
printf "𝜑(%2d) = %3d %s\n", $_, 𝜑[$_], $_ - 𝜑[$_] - 1 ?? '' !! 'Prime' for 1 .. 25;
(1e2, 1e3, 1e4, 1e5).map: -> $limit {
say "\nCount of primes <= $limit: " ~ +(^$limit).grep: {$_ == 𝜑[$_] + 1}
}
- Output:
𝜑( 1) = 1 𝜑( 2) = 1 Prime 𝜑( 3) = 2 Prime 𝜑( 4) = 2 𝜑( 5) = 4 Prime 𝜑( 6) = 2 𝜑( 7) = 6 Prime 𝜑( 8) = 4 𝜑( 9) = 6 𝜑(10) = 4 𝜑(11) = 10 Prime 𝜑(12) = 4 𝜑(13) = 12 Prime 𝜑(14) = 6 𝜑(15) = 8 𝜑(16) = 8 𝜑(17) = 16 Prime 𝜑(18) = 6 𝜑(19) = 18 Prime 𝜑(20) = 8 𝜑(21) = 12 𝜑(22) = 10 𝜑(23) = 22 Prime 𝜑(24) = 8 𝜑(25) = 20 Count of primes <= 100: 25 Count of primes <= 1000: 168 Count of primes <= 10000: 1229 Count of primes <= 100000: 9592
REXX
Version 1 Idiomatic
/*REXX program calculates the totient numbers for a range of numbers, and count primes. */
parse arg N . /*obtain optional argument from the CL.*/
if N=='' | N=="," then N= 25 /*Not specified? Then use the default.*/
tell= N>0 /*N positive>? Then display them all. */
N= abs(N) /*use the absolute value of N for loop.*/
w= length(N) /*W: is used in aligning the output. */
primes= 0 /*the number of primes found (so far).*/
/*if N was negative, only count primes.*/
do j=1 for N; T= phi(j) /*obtain totient number for a number. */
prime= word('(prime)', 1 + (T \== j-1 ) ) /*determine if J is a prime number. */
if prime\=='' then primes= primes + 1 /*if a prime, then bump the prime count*/
if tell then say 'totient number for ' right(j, w) " ──► " right(T, w) ' ' prime
end /*j*/
say
say right(primes, w) ' primes detected for numbers up to and including ' N
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
gcd: parse arg x,y; do until y==0; parse value x//y y with y x
end /*until*/; return x
/*──────────────────────────────────────────────────────────────────────────────────────*/
phi: procedure; parse arg z; if z==1 then return 1
#= 1
do m=2 for z-2; if gcd(m, z)==1 then #= # + 1
end /*m*/; return #
- output when using the default input of: 25
totient number for 1 ──► 1 totient number for 2 ──► 1 (prime) totient number for 3 ──► 2 (prime) totient number for 4 ──► 2 totient number for 5 ──► 4 (prime) totient number for 6 ──► 2 totient number for 7 ──► 6 (prime) totient number for 8 ──► 4 totient number for 9 ──► 6 totient number for 10 ──► 4 totient number for 11 ──► 10 (prime) totient number for 12 ──► 4 totient number for 13 ──► 12 (prime) totient number for 14 ──► 6 totient number for 15 ──► 8 totient number for 16 ──► 8 totient number for 17 ──► 16 (prime) totient number for 18 ──► 6 totient number for 19 ──► 18 (prime) totient number for 20 ──► 8 totient number for 21 ──► 12 totient number for 22 ──► 10 totient number for 23 ──► 22 (prime) totient number for 24 ──► 8 totient number for 25 ──► 20 9 primes detected for numbers up to and including 25
- output when using the input of: -100
25 primes detected for numbers up to and including 100
- output when using the input of: -1000
168 primes detected for numbers up to and including 1000
- output when using the input of: -10000
1229 primes detected for numbers up to and including 10000
- output when using the input of: -1000000
9592 primes detected for numbers up to and including 100000 >> 1000 seconds
Version 2 Optimized
This REXX version moved the gcd function in-line.
It is about 7% faster than the 1st version.
/*REXX program calculates the totient numbers for a range of numbers, and count primes. */
parse arg N . /*obtain optional argument from the CL.*/
if N=='' | N=="," then N= 25 /*Not specified? Then use the default.*/
tell= N>0 /*N positive>? Then display them all. */
N= abs(N) /*use the absolute value of N for loop.*/
w= length(N) /*W: is used in aligning the output. */
primes= 0 /*the number of primes found (so far).*/
/*if N was negative, only count primes.*/
do j=1 for N; T= phi(j) /*obtain totient number for a number. */
prime= word('(prime)', 1 + (T \== j-1 ) ) /*determine if J is a prime number. */
if prime\=='' then primes= primes + 1 /*if a prime, then bump the prime count*/
if tell then say 'totient number for ' right(j, w) " ──► " right(T, w) ' ' prime
end /*j*/
say
say right(primes, w) ' primes detected for numbers up to and including ' N
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
phi: procedure; parse arg z; if z==1 then return 1
#= 1
do m=2 for z-2; parse value m z with x y
do until y==0; parse value x//y y with y x
end /*until*/
if x==1 then #= # + 1
end /*m*/
return #
- output is identical to the 1st REXX version.
Version 3 Highly optimized
Libraries: How to use
Library: Functions
Library: Numbers
Library: Settings
Library: Abend
Version 1 shows that the Gcd() approach is not very efficient. Is there a better algorithm available? Yes, see procedure Totient in below program. Several other entries also implement this. This is the A-part of the program.
The task requires the calculation of all Totient(x) values up to some maximum N. Almost all entries run thru the numbers 1...N and calculate Totient(x) for each x, using Gcd() or an other method. Maybe there's an algorithm that recurring calculates Totient(x) from previous values? Yes, at the cost of some memory, see procedure Totients in below program. I found this in de Java-entry only. This is the B-part in the program.
include Settings
say version; say 'Totient function (Phi)'; say
call Time('r')
numeric digits 10; cycl. = 0; m = 7
call First25A
call PrimeCountA m
say Format(Time('e'),,3) 'seconds'; say
call Time('r')
call Totients 10**m
call First25B
call PrimeCountB m
say Format(Time('e'),,3) 'seconds'
exit
First25A:
procedure
say 'A: using calls to an optimized Totient()'; say
say ' N Phi(N) Prime?'
say Copies('-',16)
n = 0
do i = 1 to 25
p = Totient(i)
if p = i-1 then do
n = n+1; pr = 'Yes'
end
else
pr = ''
say Right(i,2) Right(p,6) pr
end
say Copies('-',16); say
say 'Found' Right(n,6) 'primes <' Right(25,8)
return
PrimeCountA:
procedure
arg x
n = 0; d = 1
do i = 1
p = Totient(i)
if p = i-1 then
n = n+1
e = Xpon(i)
if e > d then do
say 'Found' Right(n,6) 'primes <' Right(10**e,8)
if e > x-1 then
leave i
d = e
end
end
say
return
First25B:
procedure expose toti.
say 'B: generate and save all Totients, and use the stored values'; say
say ' N Phi(N) Prime?'
say Copies('-',16)
n = 0
do i = 1 to 25
p = toti.totient.i
if p = i-1 then do
n = n+1; pr = 'Yes'
end
else
pr = ''
say Right(i,2) Right(p,6) pr
end
say Copies('-',16)
say
say 'Found' Right(n,6) 'primes <' Right(25,8)
return
PrimeCountB:
procedure expose toti.
arg x
n = 0; d = 1
do i = 1
p = toti.totient.i
if p = i-1 then
n = n+1
e = Xpon(i)
if e > d then do
say 'Found' Right(n,6) 'primes <' Right(10**e,8)
if e > x-1 then
leave i
d = e
end
end
say
return
Totient:
/* Euler's totient function */
procedure expose toti.
arg x
/* Fast values */
if x < 3 then
return 1
if x < 5 then
return 2
/* Multiplicative property using Factors */
f = Factors(x)+1; fact.factor.f = 0
y = 1; v = fact.factor.1; n = 1
do i = 2 to f
a = fact.factor.i
if a = v then
n = n+1
else do
y = y*v**(n-1)*(v-1)
v = a; n = 1
end
end
return y
Totients:
/* Euler's totient numbers */
procedure expose toti.
arg x
/* Recurring sequence */
do i = 1 to x
toti.totient.i = i
end
do i = 2 to x
if toti.totient.i < i then
iterate i
do j = i by i to x
toti.totient.j = toti.totient.j-toti.totient.j/i
end
end
toti.0 = x
return x
include Functions
include Numbers
include Abend
Running with max = 100000 as in Version 1, yields
- Output:
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022 Totient function (Phi) A: using calls to an optimized Totient() N Phi(N) Prime? ---------------- 1 1 2 1 Yes 3 2 Yes 4 2 5 4 Yes 6 2 7 6 Yes 8 4 9 6 10 4 11 10 Yes 12 4 13 12 Yes 14 6 15 8 16 8 17 16 Yes 18 6 19 18 Yes 20 8 21 12 22 10 23 22 Yes 24 8 25 20 ---------------- A: Found 9 primes < 25 Found 25 primes < 100 Found 168 primes < 1000 Found 1229 primes < 10000 Found 9592 primes < 100000 1.718 seconds B: generate and save all Totients, and use the stored values N Phi(N) Prime? ---------------- 1 1 2 1 Yes 3 2 Yes 4 2 5 4 Yes 6 2 7 6 Yes 8 4 9 6 10 4 11 10 Yes 12 4 13 12 Yes 14 6 15 8 16 8 17 16 Yes 18 6 19 18 Yes 20 8 21 12 22 10 23 22 Yes 24 8 25 20 ---------------- B: Found 9 primes < 25 Found 25 primes < 100 Found 168 primes < 1000 Found 1229 primes < 10000 Found 9592 primes < 100000 0.563 seconds
That's much MUCH better compared to Version 1. Also B is considerable faster than A.
Taking max = 1 million for only Primecount gives:
- Output:
A: Found 25 primes < 100 Found 168 primes < 1000 Found 1229 primes < 10000 Found 9592 primes < 100000 Found 78498 primes < 1000000 25.337 seconds B: Found 25 primes < 100 Found 168 primes < 1000 Found 1229 primes < 10000 Found 9592 primes < 100000 Found 78498 primes < 1000000 7.423 seconds
And finally max = 10 million, only for solution B (A would probable take > 1000 seconds):
- Output:
B: Found 25 primes < 100 Found 168 primes < 1000 Found 1229 primes < 10000 Found 9592 primes < 100000 Found 78498 primes < 1000000 Found 664579 primes < 10000000 96.495 seconds
RPL
GCD approach
Very compact code, but very slow algorithm.
≪ → n
≪ 0 1 n FOR j
n j
WHILE DUP REPEAT SWAP OVER MOD END
DROP 1 == +
NEXT
≫ ≫
‘PHI’ STO
≪ 1 1 ROT FOR j
j DUP PHI - 1 == +
2 STEP
≫
'CNTPR' STO
25 PHI 10000 CNTPR
- Output:
2: 20 1: 1229
Faster version
Calculator simulator's timedog unfortunately prevents from running the code to count up to 100,000.
≪ DUP 2 OVER √
FOR j
IF DUP j MOD NOT THEN
WHILE DUP j MOD NOT REPEAT
j /
END
SWAP DUP j / - SWAP
END
IF j 2 == THEN 1 'j' STO END
2 STEP
IF DUP 1 > THEN OVER SWAP / - ELSE DROP END
≫
‘PHI’ STO
Ruby
require "prime"
def 𝜑(n)
n.prime_division.inject(1){|res, (pr, exp)| res * (pr-1) * pr**(exp-1) }
end
1.upto 25 do |n|
tot = 𝜑(n)
puts "#{n}\t #{tot}\t #{"prime" if n-tot==1}"
end
[100, 1_000, 10_000, 100_000].each do |u|
puts "Number of primes up to #{u}: #{(1..u).count{|n| n-𝜑(n) == 1} }"
end
- Output:
1 1 2 1 prime 3 2 prime 4 2 5 4 prime 6 2 7 6 prime 8 4 9 6 10 4 11 10 prime 12 4 13 12 prime 14 6 15 8 16 8 17 16 prime 18 6 19 18 prime 20 8 21 12 22 10 23 22 prime 24 8 25 20 Number of primes up to 100: 25 Number of primes up to 1000: 168 Number of primes up to 10000: 1229 Number of primes up to 100000: 9592
Rust
use num::integer::gcd;
fn main() {
// Compute the totient of the first 25 natural integers
println!("N\t phi(n)\t Prime");
for n in 1..26 {
let phi_n = phi(n);
println!("{}\t {}\t {:?}", n, phi_n, phi_n == n - 1);
}
// Compute the number of prime numbers for various steps
[1, 100, 1000, 10000, 100000]
.windows(2)
.scan(0, |acc, tuple| {
*acc += (tuple[0]..=tuple[1]).filter(is_prime).count();
Some((tuple[1], *acc))
})
.for_each(|x| println!("Until {}: {} prime numbers", x.0, x.1));
}
fn is_prime(n: &usize) -> bool {
phi(*n) == *n - 1
}
fn phi(n: usize) -> usize {
(1..=n).filter(|&x| gcd(n, x) == 1).count()
}
Output is:
N phi(n) Prime 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false Until 100: 25 prime numbers Until 1000: 168 prime numbers Until 10000: 1229 prime numbers Until 100000: 9592 prime numbers
S-BASIC
$lines
rem - return p mod q
function mod(p, q = integer) = integer
end = p - q * (p / q)
rem - return greatest common divisor of x and y
function gcd(x, y = integer) = integer
var r, temp = integer
if x < y then
begin
temp = x
x = y
y = temp
end
r = mod(x, y)
while r <> 0 do
begin
x = y
y = r
r = mod(x, y)
end
end = y
rem - return phi (also called totient) of n
function phi(n = integer) = integer
var i, count = integer
count = 1
for i = 2 to n
if gcd(n, i) = 1 then count = count + 1
next i
end = count
rem - exercise the function
var n, totient, count = integer
print " n Phi(n) Prime?"
for n = 1 to 25
totient = phi(n)
print using "#### #### ";n, totient;
if totient + 1 = n then
print "yes"
else
print "no"
next n
rem - and further test it by counting primes
print
count = 0
for n = 1 to 1000
if phi(n) = n - 1 then count = count + 1
if n = 100 then
print "Primes up to 100 = ";count
else if n = 1000 then
print "Primes up to 1000 = ";count
next n
end
- Output:
n Phi(n) Prime? 1 1 no 2 1 yes 3 2 yes 4 2 no 5 4 yes 6 2 no 7 6 yes 8 4 no 9 6 no 10 4 no 11 10 yes 12 4 no 13 12 yes 14 6 no 15 8 no 16 8 no 17 16 yes 18 6 no 19 18 yes 20 8 no 21 12 no 22 10 no 23 22 yes 24 8 no 25 20 no Primes up to 100 = 25 Primes up to 1000 = 168
Scala
The most concise way to write the totient function in Scala is using a naive lazy list:
@tailrec
def gcd(a: Int, b: Int): Int = if(b == 0) a else gcd(b, a%b)
def totientLaz(num: Int): Int = LazyList.range(2, num).count(gcd(num, _) == 1) + 1
The above approach, while concise, is not very performant. It must check the GCD with every number between 2 and num, giving it an O(n*log(n)) time complexity. A better solution is to use the product definition of the totient, repeatedly extracting prime factors from num:
def totientPrd(num: Int): Int = {
@tailrec
def dTrec(f: Int, n: Int): Int = if(n%f == 0) dTrec(f, n/f) else n
@tailrec
def tTrec(ac: Int, i: Int, n: Int): Int = if(n != 1){
if(n%i == 0) tTrec(ac*(i - 1)/i, i + 1, dTrec(i, n))
else tTrec(ac, i + 1, n)
}else{
ac
}
tTrec(num, 2, num)
}
This version is significantly faster, but the introduction of multiple recursive methods makes it far less concise. We can, however, embed the recursion into a lazy list to obtain a function as fast as the second example yet almost as concise as the first, at the cost of some readability:
@tailrec
def scrub(f: Long, num: Long): Long = if(num%f == 0) scrub(f, num/f) else num
def totientLazPrd(num: Long): Long = LazyList.iterate((num, 2: Long, num)){case (ac, i, n) => if(n%i == 0) (ac*(i - 1)/i, i + 1, scrub(i, n)) else (ac, i + 1, n)}.find(_._3 == 1).get._1
To generate the output up to 100000, Longs are necessary.
- Output:
1 (Not Prime): 1 2 (Prime) : 1 3 (Prime) : 2 4 (Not Prime): 2 5 (Prime) : 4 6 (Not Prime): 2 7 (Prime) : 6 8 (Not Prime): 4 9 (Not Prime): 6 10 (Not Prime): 4 11 (Prime) : 10 12 (Not Prime): 4 13 (Prime) : 12 14 (Not Prime): 6 15 (Not Prime): 8 16 (Not Prime): 8 17 (Prime) : 16 18 (Not Prime): 6 19 (Prime) : 18 20 (Not Prime): 8 21 (Not Prime): 12 22 (Not Prime): 10 23 (Prime) : 22 24 (Not Prime): 8 25 (Not Prime): 20 Prime Count <= N... 100: 25 1000: 168 10000: 1229 100000: 9592
SETL
program totient;
loop for n in [1..1000000] do
tot := totient(n);
if tot = n-1 then prime +:= 1; end if;
if n <= 25 then
print(lpad(str n, 2), " ",
lpad(str tot, 2), " ",
if tot = n-1 then "prime" else "" end if);
end if;
if n in [1000,10000,100000,1000000] then
print(lpad(str prime,8), "primes up to" + lpad(str n,8));
end if;
end loop;
proc totient(n);
tot := n;
i := 2;
loop while i*i <= n do
if n mod i = 0 then
loop while n mod i = 0 do
n div:= i;
end loop;
tot -:= tot div i;
end if;
if i=2 then i:=3;
else i+:=2;
end if;
end loop;
if n>1 then
tot -:= tot div n;
end if;
return tot;
end proc;
end program;- Output:
1 1
2 1 prime
3 2 prime
4 2
5 4 prime
6 2
7 6 prime
8 4
9 6
10 4
11 10 prime
12 4
13 12 prime
14 6
15 8
16 8
17 16 prime
18 6
19 18 prime
20 8
21 12
22 10
23 22 prime
24 8
25 20
168 primes up to 1000
1229 primes up to 10000
9592 primes up to 100000
78498 primes up to 1000000
Shale
#!/usr/local/bin/shale
primes library
init dup var {
pl sieve type primes::()
100000 0 pl generate primes::()
} =
go dup var {
i var
c0 var
c1 var
c2 var
c3 var
p var
" N Phi Is prime" println
i 1 =
{ i 25 <= } {
i pl isprime primes::() { "True" } { "False" } if i pl phi primes::() i "%2d %3d %s\n" printf
i++
} while
c0 0 =
c1 0 =
c2 0 =
c3 0 =
i 0 =
{ i count pl:: < } {
p i pl get primes::() =
p 100 < { c0++ } ifthen
p 1000 < { c1++ } ifthen
p 10000 < { c2 ++ } ifthen
p 100000 < { c3 ++ } ifthen
i++
} while
"" println
" Limit Count" println
c0 " 100 %5d\n" printf
c1 " 1,000 %5d\n" printf
c2 " 10,000 %5d\n" printf
c3 "100,000 %5d\n" printf
} =
init()
go()
- Output:
N Phi Is prime
1 1 False
2 1 True
3 2 True
4 2 False
5 4 True
6 2 False
7 6 True
8 4 False
9 6 False
10 4 False
11 10 True
12 4 False
13 12 True
14 6 False
15 8 False
16 8 False
17 16 True
18 6 False
19 18 True
20 8 False
21 12 False
22 10 False
23 22 True
24 8 False
25 20 False
Limit Count
100 25
1,000 168
10,000 1229
100,000 9592
Sidef
The Euler totient function is built-in as Number.euler_phi(), but we can easily re-implement it using its multiplicative property: phi(p^k) = (p-1)*p^(k-1).
func 𝜑(n) {
n.factor_exp.prod {|p|
(p[0]-1) * p[0]**(p[1]-1)
}
}
for n in (1..25) {
var totient = 𝜑(n)
printf("𝜑(%2s) = %3s%s\n", n, totient, totient==(n-1) ? ' - prime' : '')
}
- Output:
𝜑( 1) = 1 𝜑( 2) = 1 - prime 𝜑( 3) = 2 - prime 𝜑( 4) = 2 𝜑( 5) = 4 - prime 𝜑( 6) = 2 𝜑( 7) = 6 - prime 𝜑( 8) = 4 𝜑( 9) = 6 𝜑(10) = 4 𝜑(11) = 10 - prime 𝜑(12) = 4 𝜑(13) = 12 - prime 𝜑(14) = 6 𝜑(15) = 8 𝜑(16) = 8 𝜑(17) = 16 - prime 𝜑(18) = 6 𝜑(19) = 18 - prime 𝜑(20) = 8 𝜑(21) = 12 𝜑(22) = 10 𝜑(23) = 22 - prime 𝜑(24) = 8 𝜑(25) = 20
[100, 1_000, 10_000, 100_000].each {|limit|
var pi = (1..limit -> count_by {|n| 𝜑(n) == (n-1) })
say "Number of primes <= #{limit}: #{pi}"
}
- Output:
Number of primes <= 100: 25 Number of primes <= 1000: 168 Number of primes <= 10000: 1229 Number of primes <= 100000: 9592
Tiny BASIC
1 indicates prime, 0 indicates not prime.
REM PRINT THE DATA FOR N=1 TO 25
LET N = 0
10 LET N = N + 1
IF N = 26 THEN GOTO 100
GOSUB 1000
IF T = N - 1 THEN LET P = 1
IF T <> N - 1 THEN LET P = 0
PRINT N," ", T," ",P
GOTO 10
100 REM COUNT PRIMES BELOW 10000
LET C = 0
LET N = 2
110 GOSUB 1000
IF T = N - 1 THEN LET C = C + 1
IF N = 100 THEN PRINT "There are ", C, " primes below 100."
IF N = 1000 THEN PRINT "There are ", C, " primes below 1000."
IF N = 10000 THEN PRINT "There are ", C, " primes below 10000."
LET N = N + 1
IF N < 10001 THEN GOTO 110
END
1000 REM TOTIENT FUNCTION OF INTEGER N
LET M = 1
LET T = 0
1010 IF M > N THEN RETURN
LET A = N
LET B = M REM NEED TO RENAME THESE SO THEY ARE PRESERVED
GOSUB 2000
IF G = 1 THEN LET T = T + 1
LET M = M + 1
GOTO 1010
2000 REM GCD OF INTEGERS A, B
2010 IF A>B THEN GOTO 2020
LET B = B - A
IF A=0 THEN LET G = B
IF A=0 THEN RETURN
2020 LET S = A
LET A = B
LET B = S
GOTO 2010- Output:
1 1 02 1 1 3 2 1 4 2 0 5 4 1 6 2 0 7 6 1 8 4 0 9 6 0 10 4 0 11 10 1 12 4 0 13 12 1 14 6 0 15 8 0 16 8 0 17 16 1 18 6 0 19 18 1 20 8 0 21 12 0 22 10 0 23 22 1 24 8 0 25 20 0 There are 25 primes below 100. There are 168 primes below 1000.
There are 1229 primes below 10000.
VBA
Private Function totient(ByVal n As Long) As Long
Dim tot As Long: tot = n
Dim i As Long: i = 2
Do While i * i <= n
If n Mod i = 0 Then
Do While True
n = n \ i
If n Mod i <> 0 Then Exit Do
Loop
tot = tot - tot \ i
End If
i = i + IIf(i = 2, 1, 2)
Loop
If n > 1 Then
tot = tot - tot \ n
End If
totient = tot
End Function
Public Sub main()
Debug.Print " n phi prime"
Debug.Print " --------------"
Dim count As Long
Dim tot As Integer, n As Long
For n = 1 To 25
tot = totient(n)
prime = (n - 1 = tot)
count = count - prime
Debug.Print Format(n, "@@"); Format(tot, "@@@@@"); Format(prime, "@@@@@@@@")
Next n
Debug.Print
Debug.Print "Number of primes up to 25 = "; Format(count, "@@@@")
For n = 26 To 100000
count = count - (totient(n) = n - 1)
Select Case n
Case 100, 1000, 10000, 100000
Debug.Print "Number of primes up to"; n; String$(6 - Len(CStr(n)), " "); "="; Format(count, "@@@@@")
Case Else
End Select
Next n
End Sub
- Output:
n phi prime -------------- 1 1 False 2 1 True 3 2 True 4 2 False 5 4 True 6 2 False 7 6 True 8 4 False 9 6 False 10 4 False 11 10 True 12 4 False 13 12 True 14 6 False 15 8 False 16 8 False 17 16 True 18 6 False 19 18 True 20 8 False 21 12 False 22 10 False 23 22 True 24 8 False 25 20 False Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 100000 = 9592
Visual Basic .NET
Imports System.Linq.Enumerable
Module Module1
Sub Main()
For i = 1 To 25
Dim t = Totient(i)
Console.WriteLine("{0}{1}{2}{3}", i, vbTab, t, If(t = i - 1, vbTab & "prime", ""))
Next
Console.WriteLine()
Dim j = 100
While j <= 100000
Console.WriteLine($"{Range(1, j).Count(Function(x) Totient(x) + 1 = x):n0} primes below {j:n0}")
j *= 10
End While
End Sub
Function Totient(n As Integer) As Integer
If n < 3 Then
Return 1
End If
If n = 3 Then
Return 2
End If
Dim tot = n
If (n And 1) = 0 Then
tot >>= 1
Do
n >>= 1
Loop While (n And 1) = 0
End If
Dim i = 3
While i * i <= n
If n Mod i = 0 Then
tot -= tot \ i
Do
n \= i
Loop While (n Mod i) = 0
End If
i += 2
End While
If n > 1 Then
tot -= tot \ n
End If
Return tot
End Function
End Module
- Output:
1 1 2 1 prime 3 2 prime 4 2 5 4 prime 6 2 7 6 prime 8 4 9 6 10 4 11 10 prime 12 4 13 12 prime 14 6 15 8 16 8 17 16 prime 18 6 19 18 prime 20 8 21 12 22 10 23 22 prime 24 8 25 20 25 primes below 100 168 primes below 1,000 1,229 primes below 10,000 9,592 primes below 100,000
Wren
import "./fmt" for Fmt
var totient = Fn.new { |n|
var tot = n
var i = 2
while (i*i <= n) {
if (n%i == 0) {
while(n%i == 0) n = (n/i).floor
tot = tot - (tot/i).floor
}
if (i == 2) i = 1
i = i + 2
}
if (n > 1) tot = tot - (tot/n).floor
return tot
}
System.print(" n phi prime")
System.print("---------------")
var count = 0
for (n in 1..25) {
var tot = totient.call(n)
var isPrime = (n - 1) == tot
if (isPrime) count = count + 1
System.print("%(Fmt.d(2, n)) %(Fmt.d(2, tot)) %(isPrime)")
}
System.print("\nNumber of primes up to 25 = %(count)")
for (n in 26..100000) {
var tot = totient.call(n)
if (tot == n - 1) count = count + 1
if (n == 100 || n == 1000 || n%10000 == 0) {
System.print("\nNumber of primes up to %(Fmt.d(-6, n)) = %(count)")
}
}
- Output:
n phi prime --------------- 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false Number of primes up to 25 = 9 Number of primes up to 100 = 25 Number of primes up to 1000 = 168 Number of primes up to 10000 = 1229 Number of primes up to 20000 = 2262 Number of primes up to 30000 = 3245 Number of primes up to 40000 = 4203 Number of primes up to 50000 = 5133 Number of primes up to 60000 = 6057 Number of primes up to 70000 = 6935 Number of primes up to 80000 = 7837 Number of primes up to 90000 = 8713 Number of primes up to 100000 = 9592
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;
];
int N, Phi, Pwr, C, Limit;
[Text(0, "n phi is prime^m^j");
for N:= 1 to 25 do
[IntOut(0, N); ChOut(0, 9\tab\);
Phi:= Totient(N);
IntOut(0, Phi); ChOut(0, 9\tab\);
Text(0, if Phi = N-1 then "true" else "false");
CrLf(0);
];
CrLf(0);
for Pwr:= 2 to 4 do
[C:= 0;
Limit:= fix(Pow(10.0, float(Pwr)));
IntOut(0, Limit); ChOut(0, 9\tab\);
for N:= 1 to Limit do
[Phi:= Totient(N);
if Phi = N-1 then C:= C+1;
];
IntOut(0, C); CrLf(0);
];
]- Output:
n phi is prime 1 1 false 2 1 true 3 2 true 4 2 false 5 4 true 6 2 false 7 6 true 8 4 false 9 6 false 10 4 false 11 10 true 12 4 false 13 12 true 14 6 false 15 8 false 16 8 false 17 16 true 18 6 false 19 18 true 20 8 false 21 12 false 22 10 false 23 22 true 24 8 false 25 20 false 100 25 1000 168 10000 1229
zkl
fcn totient(n){ [1..n].reduce('wrap(p,k){ p + (n.gcd(k)==1) }) }
fcn isPrime(n){ totient(n)==(n - 1) }foreach n in ([1..25]){
println("\u03c6(%2d) ==%3d %s"
.fmt(n,totient(n),isPrime(n) and "is prime" or ""));
}- Output:
φ( 1) == 1 φ( 2) == 1 is prime φ( 3) == 2 is prime φ( 4) == 2 φ( 5) == 4 is prime φ( 6) == 2 φ( 7) == 6 is prime φ( 8) == 4 φ( 9) == 6 φ(10) == 4 φ(11) == 10 is prime φ(12) == 4 φ(13) == 12 is prime φ(14) == 6 φ(15) == 8 φ(16) == 8 φ(17) == 16 is prime φ(18) == 6 φ(19) == 18 is prime φ(20) == 8 φ(21) == 12 φ(22) == 10 φ(23) == 22 is prime φ(24) == 8 φ(25) == 20
count:=0;
foreach n in ([1..10_000]){ // yes, this is sloooow
count+=isPrime(n);
if(n==100 or n==1000 or n==10_000)
println("Primes <= %,6d : %,5d".fmt(n,count));
}- Output:
Primes <= 100 : 25 Primes <= 1,000 : 168 Primes <= 10,000 : 1,229
Categories:
- Programming Tasks
- Prime Numbers
- 11l
- AArch64 Assembly
- Ada
- ALGOL 68
- ALGOL-M
- APL
- ARM Assembly
- Arturo
- AutoHotkey
- AWK
- BASIC
- BCPL
- BQN
- C
- C sharp
- C++
- CLU
- COBOL
- Cowgol
- D
- Delphi
- Draco
- Dyalect
- EasyLang
- Factor
- FreeBASIC
- Fōrmulæ
- Forth
- Go
- Haskell
- J
- Java
- Jq
- Julia
- Kotlin
- Lua
- MAD
- Mathematica
- Wolfram Language
- Miranda
- Modula-2
- Nim
- Pascal
- Perl
- Ntheory
- Phix
- PL/I
- PL/M
- PicoLisp
- Python
- Quackery
- Racket
- Raku
- REXX
- RPL
- Ruby
- Rust
- S-BASIC
- Scala
- SETL
- Shale
- Sidef
- Tiny BASIC
- VBA
- Visual Basic .NET
- Wren
- Wren-fmt
- XPL0
- Zkl