#991434 (C++17) No.2798 Multiple Chain

提出ソース

結果

問題	No.2798 Multiple Chain
ユーザー	IceKnight1093
提出日時	2024-06-28 21:59:06
言語	C++17 (gcc 12.3.0 + boost 1.83.0)
結果	AC
実行時間	11 ms / 2,000 ms
コード長	11,684 bytes
コンパイル時間	3,066 ms
コンパイル使用メモリ	232,124 KB
実行使用メモリ	6,948 KB
最終ジャッジ日時	2024-06-28 21:59:11
合計ジャッジ時間	4,262 ms
ジャッジサーバーID （参考情報）	judge2 / judge5

このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示

入力	結果	実行時間実行使用メモリ
testcase_00	AC	9 ms 6,816 KB
testcase_01	AC	9 ms 6,944 KB
testcase_02	AC	9 ms 6,940 KB
testcase_03	AC	9 ms 6,944 KB
testcase_04	AC	9 ms 6,944 KB
testcase_05	AC	10 ms 6,944 KB
testcase_06	AC	11 ms 6,944 KB
testcase_07	AC	9 ms 6,940 KB
testcase_08	AC	9 ms 6,940 KB
testcase_09	AC	9 ms 6,944 KB
testcase_10	AC	9 ms 6,948 KB
testcase_11	AC	9 ms 6,940 KB
testcase_12	AC	9 ms 6,944 KB
testcase_13	AC	10 ms 6,940 KB
testcase_14	AC	9 ms 6,944 KB
testcase_15	AC	9 ms 6,940 KB
testcase_16	AC	8 ms 6,940 KB
testcase_17	AC	8 ms 6,940 KB
testcase_18	AC	9 ms 6,940 KB
testcase_19	AC	8 ms 6,944 KB
testcase_20	AC	9 ms 6,944 KB
testcase_21	AC	8 ms 6,940 KB
testcase_22	AC	8 ms 6,940 KB
testcase_23	AC	9 ms 6,944 KB
testcase_24	AC	9 ms 6,944 KB
testcase_25	AC	10 ms 6,944 KB
testcase_26	AC	11 ms 6,944 KB
testcase_27	AC	10 ms 6,944 KB
testcase_28	AC	9 ms 6,940 KB
testcase_29	AC	9 ms 6,940 KB
testcase_30	AC	8 ms 6,940 KB
testcase_31	AC	8 ms 6,940 KB
testcase_32	AC	9 ms 6,944 KB
testcase_33	AC	9 ms 6,940 KB
testcase_34	AC	9 ms 6,940 KB
testcase_35	AC	9 ms 6,940 KB
testcase_36	AC	9 ms 6,944 KB
testcase_37	AC	9 ms 6,940 KB
testcase_38	AC	8 ms 6,940 KB
testcase_39	AC	9 ms 6,940 KB
testcase_40	AC	9 ms 6,940 KB
testcase_41	AC	8 ms 6,944 KB
testcase_42	AC	9 ms 6,944 KB
testcase_43	AC	9 ms 6,944 KB
testcase_44	AC	9 ms 6,940 KB
testcase_45	AC	9 ms 6,944 KB
testcase_46	AC	9 ms 6,940 KB
testcase_47	AC	9 ms 6,944 KB
testcase_48	AC	9 ms 6,940 KB
testcase_49	AC	9 ms 6,940 KB
testcase_50	AC	9 ms 6,940 KB
testcase_51	AC	9 ms 6,940 KB
testcase_52	AC	9 ms 6,940 KB
testcase_53	AC	9 ms 6,940 KB

権限があれば一括ダウンロードができます

ソースコード

raw source code

// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#include "bits/stdc++.h"
using namespace std;
using ll = long long int;
mt19937_64 RNG(chrono::high_resolution_clock::now().time_since_epoch().count());

// Integer factorization in O(N^{1/4}
// uses squfof from msieve https://github.com/radii/msieve
// with fixes to work for n = p^3
// works up to 10^18
// probably fails on 5003^5 which is ~10^{18.5}
 
namespace NT{
    template<typename T>
    struct bigger_type{};
    template<typename T> using bigger_type_t = typename bigger_type<T>::type;
    template<> struct bigger_type<int>{using type = long long;};
    template<> struct bigger_type<unsigned int>{using type = unsigned long long;};
    //template<> struct bigger_type<int64_t>{using type = __int128;};
    //template<> struct bigger_type<uint64_t>{using type = unsigned __int128;};
 
    template<typename int_t = unsigned long long>
    struct Mod_Int{
        static inline int_t add(int_t const&a, int_t const&b, int_t const&mod){
            int_t ret = a+b;
            if(ret>=mod) ret-=mod;
            return ret;
        }
        static inline int_t sub(int_t const&a, int_t const&b, int_t const&mod){
            return add(a, mod-b);
        }
        static inline int_t mul(int_t const&a, int_t const&b, int_t const&mod){
            uint64_t ret = a * (uint64_t)b - (uint64_t)((long double)a * b / mod - 1.1) * mod;
			if(-ret < ret){
				ret = mod-1-(-ret-1)%mod;
			} else {
				ret%=mod;
			}
				
			//ret = min(ret, ret+mod);
			int64_t out = ret;
			/*if(out != a*(__int128) b % mod){
				cerr << (long double)a * b / mod << " " << (uint64_t)((long double)a * b / mod - 0.1) << "\n";
				cerr << mod << " " << ret << " " << ret+mod << " " << out << " " << (int64_t)(a*(__int128) b % mod) << "\n";
				assert(0);
			}*/
			return out;
            //return a*static_cast<bigger_type_t<int_t>>(b)%mod;
        }
        static inline int_t pow(int_t const&a, int_t const&b, int_t const&mod){
            int_t ret = 1;
            int_t base = a;
            for(int i=0;b>>i;++i){
                if((b>>i)&1) ret = mul(ret, base, mod);
                base = mul(base, base, mod);
            }
            return ret;
        }
    };
 
    template<typename T>
    typename enable_if<is_integral<T>::value, bool>::type is_prime(T x){
        if(x<T(4)) return x>T(1);
        for(T i=2;i*i<=x;++i){
            if(x%i == 0) return false;
        }
        return true;
    }
 
    template<typename T>
    typename enable_if<is_integral<T>::value, bool>::type miller_rabin_single(T const&x, T base){
        if(x<T(4)) return x>T(1);
        if(x%2 == 0) return false;
        base%=x;
        if(base == 0) return true;
 
        T xm1 = x-1;
        int j = 1;
        T d = xm1/2;
        while(d%2 == 0){ // could use __builtin_ctz
            d/=2;
            ++j;
        }
        T t = Mod_Int<T>::pow(base, d, x);
        if(t==T(1) || t==T(xm1)) return true;
        for(int k=1;k<j;++k){
            t = Mod_Int<T>::mul(t, t, x);
            if(t == xm1) return true;
            if(t<=1) break;
        }
        return false;
    }
 
    template<typename T>
    typename enable_if<is_integral<T>::value, bool>::type miller_rabin_multi(T const&){return true;}
    template<typename T, typename... S>
    typename enable_if<is_integral<T>::value, bool>::type miller_rabin_multi(T const&x, T const&base, S const&...bases){
        if(!miller_rabin_single(x, base)) return false;
        return miller_rabin_multi(x, bases...);
    }
 
    template<typename T>
    typename enable_if<is_integral<T>::value, bool>::type miller_rabin(T const&x){
        if(x < 316349281) return miller_rabin_multi(x, T(11000544), T(31481107));
        if(x < 4759123141ull) return miller_rabin_multi(x, T(2), T(7), T(61));
        return miller_rabin_multi(x, T(2), T(325), T(9375), T(28178), T(450775), T(9780504), T(1795265022));
    }
 
    template<typename T>
    typename enable_if<is_integral<T>::value, T>::type isqrt(T const&x){
        assert(x>=T(0));
        T ret = static_cast<T>(sqrtl(x));
        while(ret>0 && ret*ret>x) --ret;
        while(x-ret*ret>2*ret)
            ++ret;
        return ret;
    }
    template<typename T>
    typename enable_if<is_integral<T>::value, T>::type icbrt(T const&x){
        assert(x>=T(0));
        T ret = static_cast<T>(cbrt(x));
        while(ret>0 && ret*ret*ret>x) --ret;
        while(x-ret*ret*ret>3*ret*(ret+1))
            ++ret;
        return ret;
    }
    /*uint64_t isqrt(unsigned __int128 const&x){
        unsigned __int128 ret = sqrtl(x);
        while(ret>0 && ret*ret>x) --ret;
        while(x-ret*ret>2*ret)
            ++ret;
        return ret;
    }*/
    vector<uint16_t> saved;
    // fast prime factorization from
    // https://github.com/radii/msieve
    uint64_t squfof_iter_better(uint64_t const&x, uint64_t const&k, uint64_t const&it_max, uint32_t cutoff_div){
        if(__gcd((uint64_t)k, x)!=1) return __gcd((uint64_t)k, x);
        //cerr << "try: " << x << " " << k << "\n";
        saved.clear();
        uint64_t scaledn = k*x;
        if(scaledn>>62) return 1;
        uint32_t sqrtn = isqrt(scaledn);
        uint32_t cutoff = isqrt(2*sqrtn)/cutoff_div;
        uint32_t q0 = 1;
        uint32_t p1 = sqrtn;
        uint32_t q1 = scaledn-p1*p1;
 
        if(q1 == 0){
            uint64_t factor = __gcd(x, (uint64_t)p1);
            return factor==x ? 1:factor;
        }
 
        uint32_t multiplier = 2*k;
        uint32_t coarse_cutoff = cutoff * multiplier;
        //cerr << "at: " << multiplier << "\n";
 
        uint32_t i, j;
        uint32_t p0 = 0;
        uint32_t sqrtq = 0;
 
        for(i=0;i<it_max;++i){
            uint32_t q, bits, tmp;
 
            tmp = sqrtn + p1 - q1;
            q = 1;
            if (tmp >= q1)
                q += tmp / q1;
 
            p0 = q * q1 - p1;
            q0 = q0 + (p1 - p0) * q;
 
            if (q1 < coarse_cutoff) {
                tmp = q1 / __gcd(q1, multiplier);
 
                if (tmp < cutoff) {
                    saved.push_back((uint16_t)tmp);
                }
            }
 
            bits = 0;
            tmp = q0;
            while(!(tmp & 1)) {
                bits++;
                tmp >>= 1;
            }
            if (!(bits & 1) && ((tmp & 7) == 1)) {
 
                sqrtq = (uint32_t)isqrt(q0);
 
                if (sqrtq * sqrtq == q0) {
                    for(j=0;j<saved.size();++j){
                        if(saved[j] == sqrtq) break;
                    }
                    if(j == saved.size()) break;
                    //else cerr << "skip " << i << "\n";;
                }
            }
            tmp = sqrtn + p0 - q0;
            q = 1;
            if (tmp >= q0)
                q += tmp / q0;
 
            p1 = q * q0 - p0;
            q1 = q1 + (p0 - p1) * q;
 
            if (q0 < coarse_cutoff) {
                tmp = q0 / __gcd(q0, multiplier);
 
                if (tmp < cutoff) {
                    saved.push_back((uint16_t) tmp);
                }
            }
        }
 
        if(sqrtq == 1) { return 1;}
        if(i == it_max) { return 1;}
 
        q0 = sqrtq;
        p1 = p0 + sqrtq * ((sqrtn - p0) / sqrtq);
        q1 = (scaledn - (uint64_t)p1 * (uint64_t)p1) / (uint64_t)q0;
 
        for(j=0;j<it_max;++j) {
            uint32_t q, tmp;
 
            tmp = sqrtn + p1 - q1;
            q = 1;
            if (tmp >= q1)
                q += tmp / q1;
 
            p0 = q * q1 - p1;
            q0 = q0 + (p1 - p0) * q;
 
            if (p0 == p1) {
                q0 = q1;
                break;
            }
 
            tmp = sqrtn + p0 - q0;
            q = 1;
            if (tmp >= q0)
                q += tmp / q0;
 
            p1 = q * q0 - p0;
            q1 = q1 + (p0 - p1) * q;
 
            if (p0 == p1)
                break;
        }
        if(j==it_max) {cerr << "RNG\n"; return 1;} // random fail
 
        uint64_t factor = __gcd((uint64_t)q0, x);
        if(factor == x) factor=1;
        return factor;
    }
    uint64_t squfof(uint64_t const&x){
        static array<uint32_t, 16> multipliers{1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11};
 
        uint64_t cbrt_x = icbrt(x);
        if(cbrt_x*cbrt_x*cbrt_x == x) return cbrt_x;
 
        //uint32_t iter_lim = isqrt(isqrt(x))+10;
        uint32_t iter_lim = 300;
        for(uint32_t iter_fact = 1;iter_fact<20000;iter_fact*=4){
            for(uint32_t const&k : multipliers){
                if(numeric_limits<uint64_t>::max()/k<=x) continue; //would overflow
                uint32_t const it_max = iter_fact*iter_lim;
                uint64_t factor = squfof_iter_better(x, k, it_max, 1);
                if(factor==1 || factor==x) continue;
                return factor;
            }
        }
        cerr << "failed to factor: " << x << "\n";
        assert(0);
        assert(0);
        return 1;
    }
 
    template<typename T>
    typename enable_if<is_integral<T>::value, vector<T>>::type factorize_brute(T x){
        vector<T> ret;
        while(x%2 == 0){
            x/=2;
            ret.push_back(2);
        }
        for(uint32_t i=3;i*(T)i <= x;i+=2){
            while(x%i == 0){
                x/=i;
                ret.push_back(i);
            }
        }
        if(x>1) ret.push_back(x);
        return ret;
    }
 
    template<typename T>
    typename enable_if<is_integral<T>::value, vector<T>>::type factorize(T x){
        //cerr << "factor: " << x << "\n";
        vector<T> ret;
        const uint32_t trial_limit = 5000;
		auto trial = [&](uint32_t const&i){
			while(x%i == 0){
                x/=i;
                ret.push_back(i);
            }
		};
		trial(2);
		trial(3);
        for(uint32_t i=5, j=2;i<trial_limit && i*i <= x;i+=j, j=6-j){
            trial(i);
        }
        if(x>1){
            static stack<T> s;
            s.push(x);
            while(!s.empty()){
                x = s.top(); s.pop();
                if(!miller_rabin(x)){
                    T factor = squfof(x);
                    if(factor == 1 || factor == x){assert(0); return ret;}
                    //cerr << x << " -> " << factor << "\n";
                    s.push(factor);
                    s.push(x/factor);
                } else {
                    ret.push_back(x);
                }
            }
        }
        sort(ret.begin(), ret.end());
        return ret;
    }
}

int main()
{
    ios::sync_with_stdio(false); cin.tie(0);

    ll n; cin >> n;
    
    const int K = 70;
    vector dp(K+1, vector(K+1, vector(K+1, 0ll)));
    // dp[i][j][k] = number of shapes of length i, area j, last column k
    for (int i = 1; i <= K; ++i) dp[1][i][i] = 1;

    int ops = 0;
    for (int i = 1; i <= K; ++i) {
        for (int j = 1; j <= K; ++j) {
            for (int k = 1; k <= j; ++k) {
                for (int x = k; x <= K; ++x) {
                    dp[i][j][k] += dp[i-1][j-k][x];
                    ++ops;
                }
            }
        }
    }

    auto prm = NT::factorize(n);
    map<ll, int> freq;
    for (auto p : prm) ++freq[p];

    vector<vector<ll>> v;
    for (auto [p, x] : freq) {
        v.push_back(vector(K+1, 0ll));
        for (int i = 1; i <= K; ++i) {
            ll tot = 0;
            for (int y = 1; y <= K; ++y) tot += dp[i][x][y];
            v.back()[i] = tot + v.back()[i-1];
        }
    }
    ll ans = 0, prv = 0;
    for (int i = 1; i <= K; ++i) {
        ll ways = 1;
        for (auto & vec : v) ways *= vec[i];
        ans += ways - prv;
        prv = ways;
    }
    cout << ans << '\n';
}

yukicoder

結果

テストケース

ソースコード