Problem 95: 社交数（3個以上の友愛数）

• 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

問題 95. 社交数（3個以上の友愛数）

原文 Problem 95: Amicable chains

問題の要約：すべての要素が$10^6$以下の最長の社交数ループの要素の最小値を求めよ

社交数(Wikipedia)とはあまりなじみのない言葉ですがProblem 21: 友愛数で出てきた友愛数の3個以上のものを言うようです。

従ってまずそこで作った自分自身を除いた約数（proper divisors）の和」を求める関数sumpdivを使って$10^6$以下の整数の約数の和を配列に格納します。

その配列からループになっている友愛数のチェーンを探す関数amiChainを使って探し長さが最大の物を見つけます。

これでも1分強で終了するので、まあOKかなと思いますが。やはりさらなる高速化に挑戦します。

import sympy
def sumpdiv(n):
  return sum(sympy.divisors(n))-n

def amiChain(spd, i):         # Create a chain
  chain = []
  while i not in chain:       
    if i > N:  return []      # out of range
    chain, i = chain+[i], spd[i]
  return chain[chain.index(i):]

N = 10**6
spd = [sumpdiv(n) for n in range(0,N)]
#print(spd[:20])

maxc = 1 
for i in range(len(spd)):
  chain = amiChain(spd, i)
  if len(chain) > maxc:
    maxc, maxchain = len(chain), chain

print(min(maxchain), len(maxchain), maxchain)

約数関数の公式を使った高速化

約数関数(Wikipedia)によると約数関数$\sigma (n)$の公式はいくつかあるみたいですが以下の2つを使います。(ここでは$ \sigma(n)=\sigma_1(n)$とします)

\large 乗法定理　\sigma(ab)=\sigma(a) \times \sigma(b)\ \ \ \ \  if \ gcd(a,b) = 1　 \cdots (1)\\
\large n = p^k とすると \sigma(p^k)=1+p+p^2+\cdots +p^k　 \cdots (2)\\ \\

まず整数$n$に含まれる素数$p$から$(2)$を使って$\sigma(p^k)$を求める関数pksumを作ります。要は整数$n$に含まれる素数すべてに関して$\sigma(p^k)$を求めてそれを$(1)$の乗法定理を使って掛ければ$\sigma(n)$が求まるということになります。

def pksum(n, p):   # calculate 1+p+p^2+p^3+,,,
  ret, pk = 1,1
  while n % p == 0:
    n //= p
    pk *= p
    ret += pk
  return ret

次にこのpksumを使って[0,nmax]の約数の和を求めるsumdivrangeを作ります。エラストテネスのふるいの原理で素数の倍数に対してpksumを計算して掛けて行きます。

後で使うかもしれないので、第二引数のproperフラグで整数自身を含むかを指定できるようにしました。

def sumdivrange(nmax, proper=False):
  sdiv = (nmax + 1) * [1]
  for p in range(2, len(sdiv)):
    if sdiv[p] > 1: continue # not prime
    for pi in range(p, len(sdiv), p):  # pi = p x i
      sdiv[pi] *= pksum(pi, p)         # 1+p+p^2+p^3+,,,
  if not proper: return sdiv
  for i in range(len(sdiv)): sdiv[i] -= i   # for proper divisor, subtract the number itself
  return sdiv

これで部品はそろったので、最初のプログラムの"spd = ..."の1行を以下に変更するだけです。実行時間は約24秒とまあまあ改善出来ました。

spd = sumdivrange(N, proper=True)  # sum of proper divisors

（開発環境：Google Colab）