Totient function

From Rosetta Code
Task
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


C[edit]

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

Factor[edit]

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

Go[edit]

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.

Pascal[edit]

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[edit]

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 6[edit]

Works with: Rakudo version 2018.11

This is an incredibly inefficient way of finding prime numbers.


my \𝜑 = 0, |(1..*).hyper(:8degree).map: -> $t { +(^$t).grep: * gcd $t == 1 };
 
printf "𝜑(%2d) = %3d %s\n", $_, 𝜑[$_], $_ - 𝜑[$_] - 1 ?? '' !! 'Prime' for 1 .. 25;
 
(100, 1000, 10000).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

Python[edit]

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

REXX[edit]

/*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 to aligning the output. */
primes= 0 /*the number of primes found (so far).*/
/*if N was negative, only count primes.*/
do j=1 for N; tn= phi(j) /*obtain totient number for a number. */
prime= word('(prime)', 1 + (tn \== 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(tn, w) ' ' prime
end
say
say right(primes, w) ' primes detected for numbers up to and including ' N
exit /*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; #= z==1
do m=1 for z-1; 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

Ruby[edit]

 
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

Sidef[edit]

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

zkl[edit]

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