0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

無平方数(Squarefree number)と メビウス関数(その1)

Last updated at Posted at 2024-08-30

メビウス関数

メビウス関数(Wikipedia)は以下のように定義されます。

\begin{equation}
  \mu(n)=
  \begin{cases}
    0       & \text{$nは平方因子を持つとき$,} \\
    1       & \text{$nは偶数個の相異なる素数の積のとき$,} \\
    -1      & \text{$nは奇数個の相異なる素数の積のとき$。}
  \end{cases}
\end{equation}

$n \le 12$では

n 1 2 3 4 5 6 7 8 9 10 11 12
$\mu(n)$ 1 -1 -1 0 -1 1 -1 0 0 1 -1 0

名前は聞いたことがあるけど何に使うの?と思っている人は多いのではないでしょうか? 今回はなるほどと思えるような問題を解いていきたいと思います。

無平方数(Squarefree number)を数える

無平方数(Squarefree number)(Wikipedia)とは、$2^2=4$や$3^2=9$などの平方数で割ることのできない正の整数のことで、例えば$2 \cdot 3 \cdot 5 = 30$が無平方数となります。

【例題】N以下の無平方数の数をsf(N)と表すときsf($10^8$)を求めよ

もちろんメビウス関数を使って$\mu(n) \neq 0$のnをすべて数えても良いのですが、もう少し賢い方法があります。それをまずsf(40)の場合で考えてみます

  1. まず平方因子を持つ数を数えて、それをnから引く
  2. $n\le 40$以下の平方因子を持つ数は$2^2=4、3^2=9 または 5^2=25$の倍数
  3. この4か9か25の倍数の数はまさに「包除原理(Inclusion-exclusion principle)でProject Euler Problem 1を解く(その2)」と同様に考えることができるので、30以下のkの倍数の数を$m(k)$と表すと
\begin{align}
&sf(40) = 40- (m(4)+m(9)+m(25))\\
&+ (m(4 \times 9)+m(9 \times 25)+m(25 \times 4))\\
&- m(4 \times 9 \times 25)\\
\end{align}

と表されます。要は掛け合わせる個数が奇数なら$-$, 偶数なら$+$となるということですね。40はm(1)と表し、40を超えるものを消し、さらに$m(k^2)=m2(k)$とすると、

\begin{align}
&sf(40) = m(1)- (m(2^2)+m(3^2)+m(5^2))\\
&+ (m((2 \times 3)^2))\\
&= m2(1) - m2(2) - m2(3) - m2(5) + m2(6) \cdots (1)\\
\end{align}

この$(1)$と上のメビウス関数の値を見比べると係数が一致しているのが分かります。すなわち無平方数だけを抜き出しその素因数の個数で$\pm$が決まるわけですね。これを$\mu(n)$を使って書き、$m2(k) = \lfloor n/k^2 \rfloor$ なので

\begin{align}
&sf(40) = \mu(1)m(1) + \mu(2) m2(2) + \mu(3)  m2(3) + \mu(5)  m2(5) + \mu(6)  m2(6)\\
&= \sum_{k=1}^{\lfloor \sqrt{40} \rfloor} \mu(k) m2(k) \\
&= \sum_{k=1}^{\lfloor \sqrt{40} \rfloor} \mu(k)  \lfloor 40/k^2 \rfloor \\
\end{align}

これを一般化すると以下の式が得られます。

\begin{align}
&\large \mathrm{sf(n)} = \sum_{k=1}^{\lfloor \sqrt{n} \rfloor} \mu(k)  \lfloor n/k^2 \rfloor \\
\end{align}

これをpythonのコードにするとsympyのmobius関数を使えば1行で書くことができます。

import sympy as sp
N = 10**8
print(f"SF({N})={sum(sp.mobius(k)*(N//k**2) for k in range(1,int(N**(0.5))+1))}")
# SF(100000000)=60792694

今回は平方因子ない0の無平方数の数を考えましたが次回(その2)は、平方因子が$r$個ある場合の数を考えます。

【別解】メビウス関数を使わない解法

メビウス関数を使わない解法も見つけましたので紹介します。平方因子を含む数の個数を求めてnから引くという考え方です

\begin{align}
&sq(n,i) : n以下のiが最大平方因子の数の個数\\
&例えばsq(16,2)=|\{4,8,12\}|=3, (16は最大平方因子が4)\\ \\
&この最大平方因子iを使って平方因子を含む数は\\
&i^2j (jは無平方数) と表せる\\
&するとn以下の最大平方因子がiの個数は\\
&jの個数すなわち\lfloor n/i^2 \rfloor 以下の無平方数の個数になる\\
&式にすると\\
&sq(n,i) = sf(\lfloor n/i^2 \rfloor)\\
&sf(n)はすべての最大平方因子iで合計してnから引けばよいので\\
&sf(n) = n - \sum_{i=2}^{\sqrt{n}}sq(n,i) = n - \sum_{i=2}^{\sqrt{n}}sf(\lfloor n/i^2 \rfloor)\\
&という漸化式が得られます\\
\end{align}

これは素直にpythonのコードにすると以下のようなシンプルなものになります。

from functools import lru_cache
@lru_cache(maxsize = None)
def SqFree(n):
  return n - sum(SqFree(n//i**2) for i in range(2,int(n**(0.5))+1))
print(SqFree(10**8))
# 60792694

(開発環境:Google Colab)

この考え方はProject Euler、Problem 193: Squarefree Numbersを解くのに役に立ちます。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?