結果

問題 No.1039 Project Euler でやれ
ユーザー lam6er
提出日時 2025-04-15 23:42:32
言語 PyPy3
(7.3.15)
結果
WA  
実行時間 -
コード長 2,637 bytes
コンパイル時間 225 ms
コンパイル使用メモリ 81,800 KB
実行使用メモリ 85,076 KB
最終ジャッジ日時 2025-04-15 23:45:20
合計ジャッジ時間 3,718 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample -- * 2
other WA * 1 TLE * 1 -- * 16
権限があれば一括ダウンロードができます

ソースコード

diff # 

import sys
from math import isqrt
from collections import defaultdict

MOD = 10**9 + 7

def factorize(n):
    factors = defaultdict(int)
    while n % 2 == 0:
        factors[2] += 1
        n //= 2
    i = 3
    while i <= isqrt(n):
        while n % i == 0:
            factors[i] += 1
            n //= i
        i += 2
    if n > 1:
        factors[n] += 1
    return factors

def partitions(e):
    res = []
    def dfs(path, remaining, start):
        if remaining == 0:
            res.append(path.copy())
            return
        for i in range(start, 0, -1):
            if i > remaining:
                continue
            path.append(i)
            dfs(path, remaining - i, i)
            path.pop()
    dfs([], e, e)
    return res

def factorial(n, mod):
    res = 1
    for i in range(1, n+1):
        res = res * i % mod
    return res

def gl_size(r, p, mod):
    res = 1
    for i in range(r):
        term = pow(p, r, mod) - pow(p, i, mod)
        res = res * term % mod
    return res

def main():
    M = int(sys.stdin.readline())
    if M == 1:
        print(1)
        return

    factors = factorize(M)
    primes = list(factors.keys())
    exponents = [factors[p] for p in primes]

    fact = {}
    inv_fact = {}
    max_fact = max(exponents) if exponents else 0
    current = 1
    for i in range(max_fact + 1):
        fact[i] = current
        current = current * (i + 1) % MOD

    total = 0
    from itertools import product

    parts_list = []
    for p, e in zip(primes, exponents):
        parts = partitions(e)
        parts_list.append( [(p, part) for part in parts] )

    combinations = product(*parts_list)
    for comb in combinations:
        aut = 1
        for (p, part) in comb:
            part = sorted(part, reverse=True)
            same = all(x == part[0] for x in part)
            if same and part[0] == 1:
                r = len(part)
                aut_p = gl_size(r, p, MOD)
            else:
                aut_p = 1
                phi = 1
                for k in part:
                    pk = pow(p, k, MOD)
                    phi = phi * (pk - pk // p) % MOD
                aut_p = phi
                hom = 1
                for i in range(len(part)):
                    for j in range(i):
                        kj = part[i]
                        hom = hom * pow(p, kj, MOD) % MOD
                aut_p = aut_p * hom % MOD
            aut = aut * aut_p % MOD
        m_fact = factorial(M, MOD)
        inv_aut = pow(aut, MOD-2, MOD)
        res = m_fact * inv_aut % MOD
        total = (total + res) % MOD

    print(total)

if __name__ == '__main__':
    main()
0