結果

問題 No.1039 Project Euler でやれ
ユーザー qwewe
提出日時 2025-05-14 13:01:25
言語 PyPy3
(7.3.15)
結果
WA  
実行時間 -
コード長 9,227 bytes
コンパイル時間 337 ms
コンパイル使用メモリ 81,896 KB
実行使用メモリ 86,360 KB
最終ジャッジ日時 2025-05-14 13:03:40
合計ジャッジ時間 2,185 ms
ジャッジサーバーID
(参考情報)
judge1 / judge4
このコードへのチャレンジ
(要ログイン)
ファイルパターン 結果
sample AC * 2
other AC * 15 WA * 3
権限があれば一括ダウンロードができます

ソースコード

diff #

import sys
import itertools

# Set a higher recursion depth limit for potentially deep recursion in partition generation
# Although max exponent K is relatively small (~20 for M < 10^6), the structure of recursion might need it.
sys.setrecursionlimit(2000) 

# Define the modulus P = 10^9 + 7
MOD = 10**9 + 7

# Function for modular exponentiation (base^power % MOD)
# Uses binary exponentiation (also known as exponentiation by squaring) for efficiency
def fast_pow(base, power):
    result = 1
    while power > 0:
        # If power is odd, multiply result with base
        if power % 2 == 1:
            result = (result * base) % MOD
        # Square the base and halve the power
        base = (base * base) % MOD
        power //= 2
    return result

# Function for modular inverse using Fermat's Little Theorem
# Calculates a^(MOD-2) % MOD, which is the modular inverse of a modulo MOD (since MOD is prime)
def inverse(a):
    # Check if a is 0 modulo MOD. If so, its inverse doesn't exist.
    # The size of an automorphism group |Aut(G)| is always >= 1 for |G| >= 1, so a should not be 0.
    if a == 0: 
         raise ValueError("Modular inverse of 0 requested")
    return fast_pow(a, MOD - 2)

# Precompute factorials modulo MOD up to the maximum possible value of M
# M < 10^6, so we compute up to 10^6.
MAX_M_limit = 10**6 
fact = [1] * (MAX_M_limit + 1)
for i in range(1, MAX_M_limit + 1):
    fact[i] = (fact[i - 1] * i) % MOD

# Function to calculate the size of the automorphism group of an abelian p-group
# G = Z_{p^{e_1}} x ... x Z_{p^{e_n}}
# The input 'partition' is a list of exponents [e1, e2, ..., en]
# The formula used requires the partition structure represented as:
# G = (Z_{p^{lambda_1}})^{m_1} x ... x (Z_{p^{lambda_k}})^{m_k} where lambda_1 > lambda_2 > ... > lambda_k > 0
# This function handles converting the input partition list to this structure internally.
def calculate_aut_p_group_size(p, partition):
    # Count occurrences of each exponent value in the input partition
    counts = {}
    for x in partition:
        if x > 0: # Partition parts must be positive integers
          counts[x] = counts.get(x, 0) + 1
    
    # If the partition is empty (this corresponds to k=0 exponent, for M=1), 
    # the group is the trivial group {e}, and |Aut({e})| = 1.
    if not counts: 
        return 1

    # Get the distinct exponents (lambda_i) and sort them in descending order
    distinct_exponents = sorted(counts.keys(), reverse=True)
    
    lambdas = distinct_exponents
    # Get the multiplicities (m_i) corresponding to the sorted distinct exponents
    ms = [counts[lam] for lam in lambdas]
    k = len(lambdas) # Number of distinct exponent types

    # The formula for |Aut(G)| has three main product terms. Calculate them modularly.
    
    # Term 1: Product related to GL(m_i, F_p) factors
    term1 = 1
    for i in range(k): # Loop over distinct exponent types lambda_i
        m_i = ms[i]
        # Calculate p^{m_i} mod MOD
        p_pow_m_i = fast_pow(p, m_i)
        
        current_prod = 1
        # Inner product runs from j=0 to m_i-1
        for j in range(m_i): 
            # Calculate p^j mod MOD
            p_pow_j = fast_pow(p, j)
            # Calculate (p^{m_i} - p^j) mod MOD. Add MOD before taking modulo to handle potential negative result.
            val = (p_pow_m_i - p_pow_j + MOD) % MOD 
            current_prod = (current_prod * val) % MOD
        term1 = (term1 * current_prod) % MOD

    # Term 2: Product over pairs of distinct exponents (1 <= i < j <= k)
    term2 = 1
    for i in range(k):
        for j in range(i + 1, k): # j > i
            lambda_j = lambdas[j]
            m_i = ms[i]
            m_j = ms[j]
            
            # Calculate (p^{lambda_j})^(m_i * m_j) mod MOD
            factor1 = fast_pow(p, lambda_j) 
            factor1_pow = m_i * m_j # Exponent is relatively small, direct computation is fine
            term2_part1 = fast_pow(factor1, factor1_pow) 
            
            # Calculate (p^{m_j})^(m_i * m_j) mod MOD
            factor2 = fast_pow(p, m_j) 
            factor2_pow = m_i * m_j
            term2_part2 = fast_pow(factor2, factor2_pow) 
            
            # Accumulate product terms modulo MOD
            term2 = (term2 * term2_part1) % MOD
            term2 = (term2 * term2_part2) % MOD

    # Term 3: Product related to powers p^{lambda_i - 1}
    term3 = 1
    for i in range(k): # Loop over distinct exponent types lambda_i
        lambda_i = lambdas[i]
        m_i = ms[i]
        
        # Base is p^{lambda_i - 1}
        factor = fast_pow(p, lambda_i - 1) 
        # Exponent is m_i^2
        exponent = m_i * m_i
        # Calculate (p^{lambda_i - 1})^(m_i^2) mod MOD
        term3_part = fast_pow(factor, exponent) 
        
        # Accumulate product term modulo MOD
        term3 = (term3 * term3_part) % MOD

    # Combine all three terms to get the final size of the automorphism group
    total_aut_size = (term1 * term2) % MOD
    total_aut_size = (total_aut_size * term3) % MOD
    
    return total_aut_size

# Function for prime factorization of an integer n
def prime_factorize(n):
    factors = {}
    d = 2
    # Use trial division up to sqrt(n)
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    # If n is still greater than 1 after the loop, it must be a prime factor itself
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors

# Function to generate integer partitions of k with memoization
# Generates partitions as lists of integers in non-decreasing order, e.g., [1, 1, 2] for k=4.
memo_partitions = {} 
def generate_partitions_memo(k): 
    # Use k as the key for memoization
    state = k 
    if state in memo_partitions:
        return memo_partitions[state]
    
    # Base case: Partition of 0 is represented by a list containing an empty list
    if k == 0:
        return [[]]
    
    partitions = []
    # 'i' represents the smallest part in the partition being constructed
    for i in range(1, k + 1):
        # Recursively find partitions for the remaining sum k-i
        results_k_minus_i = generate_partitions_memo(k - i)
        for result in results_k_minus_i:
            # Ensure non-decreasing order:
            # The current part 'i' must be less than or equal to the smallest part in the 'result' list (if 'result' is not empty)
            if not result or result[0] >= i:
                 # If the condition is met, prepend 'i' to the sub-partition 'result'
                 partitions.append([i] + result)

    # Store the generated partitions in the memoization table and return them
    memo_partitions[state] = partitions
    return partitions

# Main program logic
# Read input M
M = int(sys.stdin.readline())

# Handle the trivial case M=1. The only group is the trivial group on {0}. There is 1 such structure.
if M == 1:
    print(1) 
    sys.exit() # Exit after printing the result

# Compute the prime factorization of M. Example: M=12 -> {2: 2, 3: 1}
prime_factors = prime_factorize(M)
# Get list of prime factors p_i
primes = list(prime_factors.keys())
# Get list of corresponding exponents k_i
exponents = [prime_factors[p] for p in primes]

# Generate all partitions for each exponent k_i. Store these lists of partitions.
partitions_for_exponents = []
for k_i in exponents:
    partitions_for_exponents.append(generate_partitions_memo(k_i))

# Calculate M! mod MOD using the precomputed table
M_factorial = fact[M]

# Initialize the total sum of terms M! / |Aut(G)|
total_sum = 0

# An abelian group G of order M is isomorphic to a direct product of its Sylow p-subgroups G_p.
# The structure of each G_p is determined by a partition of the exponent k_i of p in M.
# We need to iterate through all combinations of partitions, one for each prime factor p_i.
# `itertools.product` generates the Cartesian product of the lists of partitions.
list_of_partition_lists = partitions_for_exponents
combinations = list(itertools.product(*list_of_partition_lists))

# Each 'combo' is a tuple of partitions, e.g., (partition_for_p1, partition_for_p2, ...).
# This tuple represents one isomorphism class of an abelian group G of order M.
for combo in combinations:
    # Calculate the size of the automorphism group |Aut(G)| for this isomorphism class.
    # |Aut(G)| = product of |Aut(G_p)| for each Sylow p-subgroup G_p.
    current_aut_size = 1
    for i in range(len(primes)):
        p = primes[i] # The prime factor
        partition = combo[i] # The partition defining the structure of G_p
        # Calculate |Aut(G_p)| using the dedicated function
        aut_p_group_size = calculate_aut_p_group_size(p, partition)
        # Multiply into the total |Aut(G)| size, modulo MOD
        current_aut_size = (current_aut_size * aut_p_group_size) % MOD
    
    # Calculate the term M! / |Aut(G)| modulo MOD.
    # This is equivalent to M! * (|Aut(G)|)^(-1) mod MOD.
    inv_aut_size = inverse(current_aut_size)
    term = (M_factorial * inv_aut_size) % MOD
    # Add this term to the total sum, modulo MOD
    total_sum = (total_sum + term) % MOD

# Print the final total sum modulo MOD
print(total_sum)
0