オイラーのトーシェント関数

この記事は ひとりアドベントカレンダーRosettaCodeで楽しむプログラミング Advent Calendar 2025の20日めの記事です。

定義

• 与えられた整数 $n$ までの整数で、$n$ と互いに素なものの個数
• $1 \leq k \leq n$ の整数 $k$ で $\gcd(n,k) = 1$ であるものの個数

ちなみに、$\varphi(N) = N - 1$ ならば $N$ は素数

タスク

• $1 \leq i \leq 25$ について、$i$, $\phi(i)$, $i$ が素数かどうか、を示せ。
• （$\varphi()$ を用いた素数判定により）100, 1,000, 10,000, 100,000（オプション）以下の素数の個数をそれぞれ示せ。

考える

Wikipedia に、$n = \prod p_k^{e_k}$ と素因数分解できるとき、$\displaystyle \varphi(n) = n \prod \left (1 - \frac{1}{p_k} \right)$ とある。
$p_k$ は $n$ の素因数なので一度は割ることができ、それを1小さい値に差し替えた $\prod (p_k - 1) \cdot p_k^{e_k - 1}$ なので $n$ より小さい。

元の方の式では素因数 $p_k$ だけを使いその個数 $e_k$ を使っていないので、こちらをベースにして考える。
より手間を減らすために、$X \cdot (1 - \frac{1}{P}) = X - X/P$ という関係に注意して、

import Data.Numbers.Primes
import Data.List

totient n = foldl' step n $ map head $ group $ primeFactors n
  where
    step x p = x - div x p

-- さらに圧縮
totient n = foldl' (\x p -> x - div x p) n . map head . group . primeFactors $ n

しかしタスクは、これを素数判定の方法として使うことを意図しているので、
Data.Numbers.Primes を使うのはズルだろう。
それがあるなら totient なしで length $ takeWhile (ub >=) primes だけで済んでしまう。

素因数分解は、2、3以降の素数で順に試し割りをしていくことで実現できる。
手持ちのコードはこんな感じ：

primeFactors :: Int -> [Int]
primeFactors n = loop n (2 : [3, 5 ..])
  where
    loop m pps@(p:ps)                -- 擬素数列 pps と、割る数 m について
      | m == 1    = []               -- mは1まで割りきっていたら、もう素因数はない
      | m < p * p = [m]              -- m <= √p まで調べたらそれ以上に素因数があるなら m 自体だけ
      | r == 0    = p : loop q pps   -- mがpで割りきれたらpは素因数 pps は進めずループ
      | otherwise = loop m ps        -- 割りきれないのなら次の擬素数を試しにいく
      where
        (q,r) = divMod m p           -- 試し割り

Haskell的には unfoldr で書くのだろうか。

primeFactors n = unfoldr step (n, 2 : [3, 5 ..])
  where
    step (m, pps@(p:ps))
      | m == 1    = Nothing
      | m < p * p = Just (m, (1, pps)) -- dropWhile (m >) ps かもしれないが、やる意味はない
      | r == 0    = Just (p, (q, pps))
      | otherwise = step (m, ps)       -- 中でループして次の素数を試しに行く
      where
        (q,r) = divMod m p

RosettaCodeのコード

上の totient と primeFactors がfusionしたような構造だった。

{-# LANGUAGE BangPatterns #-}

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                                     -- 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             -- でもこれだと毎回 i == 2 のチェック入るよ。
          | otherwise = 2 + i
        m_ = nextM m                 -- mがiで割り切れたとき↓
        nextM !x                     -- 割れるだけ割ってiをmから無くしてから次に進む
          | x `mod` i == 0 = nextM $ x `div` i
          | otherwise = x
        tot_ = tot - (tot `div` i)   -- 素因数 i について答えを修正する計算

教科書的には、状態が変化した変数は i に対して i' とする感じだと思うけど、アンダースコアを付けるこのスタイルはどこの流儀なのだろう？
自分のブラウザ設定とフォントとあいまって、サイトでめちゃくちゃ見づらい。
自分はプライム付けるのも嫌いなので番号を振るスタイルにしている。

Bang Patternを思考停止で付けてるけど、ガード i * i > m を判定するために、m, i は速攻で評価されるので、tot に付けたもの以外は無意味だと思う。

nextM の x も、ガードを判定するためには値を求める必要があるし。

i の計算で毎回 i == 2 を避けるためには、外に追い出すのが簡単。

という自分の好みで直すとこんな感じ：

totient :: Integral a => a -> a
totient 0 = 1
totient n = loop (abs n) (abs n) (2 : [3, 5 ..])
  where
    loop m tot (p:ps)
      | m < p * p = if m == 1 then tot else tf m
      | divided   = loop m1 (tf p) ps
      | otherwise = loop m  tot    ps
      where
        divided = mod m p == 0
        tf d = tot - div tot d
        m1 = nextM m
        nextM x = case divMod x p of
          (q, 0) -> nextM q
          _      -> x

余談：nextM

nextM を until で書こうとすると、

m1 = (\(q, r) -> p * q + r) $ until ((0 <) . snd) (flip divMod p . fst) (m, 0)

終了条件を「余りが出た」にすると、余りが出たときには割りすぎているので、1ステップ戻すことになってかっこわるい。

m1 = until ((0 <) . flip mod p) (flip div p) m

終了条件を「もう割れない、次割ると余りが出る」にすると、div と mod が分かれてしまう。

余談：flip

flip div p を (`div` p) としろ、と VSCode もごちゃごちゃ言ってくるのだけど、

• 関数をバッククオートで囲むことで演算子と見なす、という構文上の例外措置
• セクションを作ることで演算子を関数と見なす、という構文上の例外措置

という相反するものを同時に使うとか、ツリーにタッチしているとしか思えない。
\x -> div x p なら何の問題もない。

OEIS

OEIS A000010 にもHaskellのコードがあった。

a n = length (filter (==1) (map (gcd n) [1..n])) -- Allan C. Wechsler, Dec 29 2014

そうだけど、何の役にも立たない記述だな。
Juliaのやつの方が面白い。

# Computes the first N terms of the sequence.
function A000010List(N)
    phi = [i for i in 1:N + 1]
    for i in 2:N + 1
        if phi[i] == i
            for j in i:i:N + 1
                phi[j] -= div(phi[j], i)
    end end end
return phi end

エラトステネスのふるいを作って、初期値はインデックスと同じにしておく。
小さい方から見ていって、そのまま手つかずな番号は素数なので、自身を含めその全ての倍数に対して $1 - \frac{1}{p}$ を掛ける。
通り過ぎた枠は、$\varphi()$ の値が完成している。

余談：セクション

可換則の成り立つ演算の演算子について、セクションを (== 1) と書く人が多いが、
それ (\x -> equal x 1) と同じでわざわざ第1引数に差し込んでいる。
(1 ==) なら (\x -> equal 1 x) つまり equal 1 なので、自分はこっちの方が好き。

余談：Project Euler

トーシェント関数と言えば Project Euler 214 が直球でそれ。