結果

問題 No.3038 シャッフルの再現
ユーザー qwewe
提出日時 2025-05-14 13:01:33
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
RE  
実行時間 -
コード長 8,980 bytes
コンパイル時間 1,055 ms
コンパイル使用メモリ 86,548 KB
実行使用メモリ 7,848 KB
最終ジャッジ日時 2025-05-14 13:03:47
合計ジャッジ時間 3,938 ms
ジャッジサーバーID
(参考情報)
judge1 / judge4
このコードへのチャレンジ
(要ログイン)
ファイルパターン 結果
sample RE * 1
other RE * 21
権限があれば一括ダウンロードができます

ソースコード

diff #

#include <iostream>
#include <vector>
#include <map>
#include <utility> // For std::pair

// Use __int128 for intermediate calculations in modular exponentiation and matrix multiplication
// to prevent overflow, especially with large moduli or intermediate values.
// This type is supported by GCC and Clang. Check platform compatibility if using other compilers.
#ifdef __GNUC__
using LargeInt = __int128;
#else
// Fallback type if __int128 is not available. This might lead to incorrect results due to overflow.
using LargeInt = long long; 
#warning "__int128 not supported, using long long as fallback. Potential overflow issues."
#endif

// Modular exponentiation: computes (base^exp) % modulus
long long power_mod(long long base, long long exp, long long modulus) {
    long long res = 1;
    base %= modulus;
    // Ensure base is non-negative
    if (base < 0) base += modulus; 
    while (exp > 0) {
        if (exp % 2 == 1) res = (LargeInt)res * base % modulus;
        base = (LargeInt)base * base % modulus;
        exp /= 2;
    }
    return res;
}

// Computes the Legendre symbol (a/p), assuming p is an odd prime.
// Returns 1 if a is a quadratic residue modulo p, -1 if a quadratic non-residue, 0 if a is divisible by p.
long long legendre(long long a, long long p) {
    a %= p;
    // Ensure a is in the range [0, p-1]
    if (a < 0) a += p; 
    if (a == 0) return 0;
    
    // Use Euler's criterion: (a/p) = a^((p-1)/2) mod p
    long long res = power_mod(a, (p - 1) / 2, p);
    
    // The result must be 1 or p-1 (which represents -1 mod p)
    if (res == p - 1) return -1; 
    // If res is not p-1, it must be 1 due to properties of modular arithmetic.
    return 1; 
}


// Structure to represent a 2x2 matrix
struct Matrix {
    long long mat[2][2];
};

// Multiplies two 2x2 matrices A and B modulo modulus
Matrix multiply(Matrix A, Matrix B, long long modulus) {
    Matrix C;
    for(int i=0; i<2; ++i) {
        for(int j=0; j<2; ++j) {
            C.mat[i][j] = 0;
            for(int k=0; k<2; ++k) {
                 // Use LargeInt for intermediate product to avoid overflow
                 C.mat[i][j] = (C.mat[i][j] + (LargeInt)A.mat[i][k] * B.mat[k][j]) % modulus;
            }
        }
    }
    return C;
}

// Computes matrix A raised to the power exp modulo modulus using binary exponentiation
Matrix matrix_power(Matrix A, long long exp, long long modulus) {
    Matrix res;
    // Initialize result matrix as Identity matrix
    res.mat[0][0] = 1; res.mat[1][1] = 1;
    res.mat[0][1] = 0; res.mat[1][0] = 0; 
    
    // Reduce initial matrix A elements modulo modulus. This ensures values are in range [0, modulus-1].
    A.mat[0][0] %= modulus; if(A.mat[0][0] < 0) A.mat[0][0] += modulus;
    A.mat[0][1] %= modulus; if(A.mat[0][1] < 0) A.mat[0][1] += modulus;
    A.mat[1][0] %= modulus; if(A.mat[1][0] < 0) A.mat[1][0] += modulus;
    A.mat[1][1] %= modulus; if(A.mat[1][1] < 0) A.mat[1][1] += modulus;

    while (exp > 0) {
        if (exp % 2 == 1) res = multiply(res, A, modulus);
        A = multiply(A, A, modulus);
        exp /= 2;
    }
    return res;
}

// Computes the prime factorization of a positive integer n using trial division
// Returns a vector of pairs {prime, exponent}
std::vector<std::pair<long long, int>> prime_factorize(long long n) {
    std::vector<std::pair<long long, int>> factors;
    for (long long i = 2; i * i <= n; ++i) {
        if (n % i == 0) {
            int count = 0;
            while (n % i == 0) {
                n /= i;
                count++;
            }
            factors.push_back({i, count});
        }
    }
    // If n is still greater than 1, it must be a prime factor itself
    if (n > 1) {
        factors.push_back({n, 1});
    }
    return factors;
}

// Computes the Pisano period pi(p) for a prime p
long long compute_pisano_prime(long long p) {
    // Handle base cases for p=2 and p=5
    if (p == 2) return 3;
    if (p == 5) return 20;

    // Determine the upper bound B for pi(p) based on Legendre symbol (5/p)
    long long B;
    if (legendre(5, p) == 1) {
        // If 5 is a quadratic residue mod p, pi(p) divides p-1
        B = p - 1;
    } else { 
        // If 5 is a quadratic non-residue mod p, pi(p) divides 2(p+1)
        B = 2 * (p + 1);
    }

    // Find the prime factorization of the bound B
    std::vector<std::pair<long long, int>> factors = prime_factorize(B);
    
    // Start with pi = B and reduce it by prime factors q until A^(pi/q) is not Identity
    long long pi = B;

    // The Fibonacci matrix M = {{1,1},{1,0}}
    Matrix A_base; 
    A_base.mat[0][0] = 1; A_base.mat[0][1] = 1;
    A_base.mat[1][0] = 1; A_base.mat[1][1] = 0;
    
    // Identity matrix I
    Matrix I; 
    I.mat[0][0] = 1; I.mat[0][1] = 0;
    I.mat[1][0] = 0; I.mat[1][1] = 1;

    // Iterate through each distinct prime factor q of B
    for (const auto& factor_pair : factors) {
        long long q = factor_pair.first;
        // Repeatedly divide pi by q as long as pi is divisible by q and A^(pi/q) mod p is the Identity matrix
        while (pi > 0 && pi % q == 0) { 
            long long exponent_check = pi / q;
            // Prevent checking exponent 0, which should not happen with correct logic but safe check
            if (exponent_check == 0) break; 
            
            // Compute A^(pi/q) mod p
            Matrix A_pow = matrix_power(A_base, exponent_check, p);
            
            // Check if A_pow is the Identity matrix
            bool is_identity = true;
            for(int r=0; r<2; ++r) {
                 for(int c=0; c<2; ++c) {
                     if (A_pow.mat[r][c] != I.mat[r][c]) {
                         is_identity = false;
                         break;
                     }
                 }
                 if (!is_identity) break;
            }

            if (is_identity) {
                // If it is Identity, then the period is smaller. Reduce pi by factor q.
                pi /= q; 
            } else {
                // If it's not Identity, we found the correct power for this prime factor q.
                // Stop dividing pi by q and move to the next prime factor.
                break; 
            }
        }
    }
    // The final value of pi is the smallest period.
    return pi;
}


int main() {
    // Use faster I/O operations
    std::ios_base::sync_with_stdio(false);
    std::cin.tie(NULL);

    int N; // Number of prime factors of M
    std::cin >> N;

    // Read the prime factorization of M: pairs of {prime p_i, exponent k_i}
    std::vector<std::pair<long long, int>> PK(N);
    for (int i = 0; i < N; ++i) {
        std::cin >> PK[i].first >> PK[i].second;
    }

    // Map to store the maximum exponent for each prime factor found in the LCM calculation
    std::map<long long, int> max_powers; 
    // The final modulus for the result
    long long P_mod = 1000000007;

    // Process each prime power factor p_i^k_i of M
    for (int i = 0; i < N; ++i) {
        long long p = PK[i].first;
        int k = PK[i].second;

        // Compute the Pisano period for the prime base p
        long long pi_p = compute_pisano_prime(p);
        
        // Use the property pi(p^k) = p^(k-1) * pi(p) (assuming Wall's conjecture holds, widely believed true)
        // Compute the prime factorization of pi(p^k)
        std::map<long long, int> current_factors_map;
        
        // Add factors from pi(p)
        std::vector<std::pair<long long, int>> factors_pi_p = prime_factorize(pi_p);
        for(const auto& factor : factors_pi_p) {
             current_factors_map[factor.first] += factor.second;
        }

        // Add the factor p^(k-1)
        if (k > 1) {
             current_factors_map[p] += (k - 1);
        }
        
        // Update the global max_powers map with the maximum exponents encountered so far for each prime factor
        for(const auto& pair : current_factors_map) {
            long long prime_factor = pair.first;
            int exponent = pair.second;
            // Update if this prime factor is new or has a larger exponent than previously seen
            if (max_powers.find(prime_factor) == max_powers.end() || exponent > max_powers[prime_factor]) {
                max_powers[prime_factor] = exponent;
            }
        }
    }

    // Compute the final LCM by multiplying together prime factors raised to their maximum required exponents, modulo P_mod
    long long final_lcm = 1;
    for (const auto& pair : max_powers) {
        long long q = pair.first;
        int exponent = pair.second;
        // Calculate q^exponent mod P_mod
        if (exponent > 0) { // Only consider factors with positive exponent
           long long term = power_mod(q, exponent, P_mod);
           // Multiply into the final LCM result modulo P_mod
           final_lcm = (LargeInt)final_lcm * term % P_mod;
        }
    }

    // Output the final result
    std::cout << final_lcm << std::endl;

    return 0;
}
0