オイラーの関数をpythonで計算する

自然数nに対して、$\phi$(n)はn以下の自然数でnと互に素なものの個数を与える関数をオイラーの関数といいます。Pyhonを使って、オイラーの関数を計算してみます。

nとある数が互いに素ということは最小公倍数が1ということです。ということは以前に作ったgcd関数を使うことができます。次のgcd関数はユークリッドの互除法を使って最小公倍数を求めることができます。
参考：Pythonで最小公倍数、最大公約数を計算する

#1 ユークリッドの互除法による最大公約数を求める関数
def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

ここで力技でオイラーの関数の計算をしてみます。

#2 gcd関数を使ってオイラーの関数を計算する
def phi(num):
    cnt=0
    for i in range (1,num):
        if gcd(num,i)==1:
            cnt+=1
    return cnt
    
phi(72)    

答えは24になります。

答えが合っているかは数えればよいのですが、大変なのであらかじめ用意されている関数を使って計算します。SymPyには整sympy.ntheory数論論モジュールにtotient関数があります。

#3 sympy.ntheoryモジュールのtotient関数でオイラーの関数を計算する
n=72
import sympy.ntheory.factor_ 
totient_n =  sympy.ntheory.factor_.totient(n)  
print(totient_n)     
print("phi({}) =  {} ".format(n, totient_n)) # 1 5 7 11 13 17 19 23 

phi(72) =  24 

どうやら正しいようです。この関数は英語でEuler’s totient functionというそうです。totientとは見慣れぬ言葉ですが、辞書で引いてもいきなり「互いに素」という訳語が出てきます。totとはラテン語でたくさんという意味らしいです。

自然数nの素因数分解を、$p_1^{k_1}p_2^{k_2}\cdots p_m^{k_m}$ とすると

$\phi(n)=(p_1^{k_1}-p_1^{k_1-1})(p_2^{k_2}-p_2^{k_2-1})\cdots(p_m^{k_m}-p_m^{k_m-1})$
$\displaystyle=n\left(1-\frac{1}{p_1}\right)\left(1-\frac{1}{p_2}\right)\cdots\left(1-\frac{1}{p_m}\right)$

素因数分解をする関数は以前に作成しました。
参考：Pythonによる約数の計算と素因数分解

#4 素因数分解を辞書形式で求める関数 
def prime_f(num):
    divisors={}
    i=2
    while(i <= num):
        k=0
        while(num % i == 0):
            k+=1
            num //= i
        if k != 0:
            divisors[i]=k
        i += 1
    return divisors

{2: 3, 3: 2}

結果は、$2^3×3^2$を辞書形式で計算されます。これを式にあてはめるプログラムは次の通りです。なお、この関数は前のprime_fが実行されている必要があります。

#5 素因数分解の結果の辞書からオイラーの関数を求める関数
def totient_n(num):
    phi=num
    for i in prime_f(num).keys():
        phi*=(1-1/i)
    return int(phi)

totient_n(72)

答えは同じ24です。この考え方は中国の剰余定理を使っています。
参考：中国の剰余定理をPythonで計算する

このほかにeulerlibというのがあるようです。これを使うといろいろできそうですが、JupyterNotebookにはあらかじめ用意されていないようです。