結果
問題 | No.2365 Present of good number |
ユーザー |
|
提出日時 | 2023-06-30 22:56:50 |
言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
結果 |
AC
|
実行時間 | 25 ms / 2,000 ms |
コード長 | 4,045 bytes |
コンパイル時間 | 3,179 ms |
コンパイル使用メモリ | 238,644 KB |
最終ジャッジ日時 | 2025-02-15 04:29:01 |
ジャッジサーバーID (参考情報) |
judge4 / judge3 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
sample | AC * 2 |
other | AC * 39 |
ソースコード
#include <bits/stdc++.h> using namespace std; using ll = long long; #define rep(i,n) for (int i=0;i<(int)(n);i++) // https://trap.jp/post/1444/ by tatyam constexpr uint32_t totient(uint32_t x){ uint32_t ans = x; for(uint32_t i = 2; i * i <= x; i++) if(x % i == 0){ ans /= i; ans *= i - 1; do x /= i; while(x % i == 0); } if(x != 1){ ans /= x; ans *= x - 1; } return ans; } template<uint32_t P> struct Modint{ static_assert(P < 0x80000000, "P must be smaller than 2^31"); uint32_t a; Modint<totient(P)> b; private: static uint32_t mod(uint64_t x){ if(x < P * 2) return uint32_t(x); return uint32_t(x % P) + P; } static uint32_t mul(uint32_t a, uint32_t b){ return mod(uint64_t(a) * b); } static uint32_t pow(uint32_t a, uint32_t b){ uint32_t ans = 1; while(b){ if(b & 1) ans = mul(ans, a); a = mul(a, a); b >>= 1; } return ans; } public: Modint(uint64_t x): a(mod(x)), b(x){} Modint(uint32_t a, Modint<totient(P)> b): a(a), b(b){} Modint(): a(0), b(0){} uint32_t val() const { if(a < P) return a; return a - P; } Modint& operator*=(const Modint& other){ a = mul(a, other.a); b *= other.b; return *this; } Modint operator*(const Modint& other) const { return Modint(*this) *= other; } Modint& operator+=(const Modint& other){ a += other.a; if(a >= P * 2) a -= P; if(a >= P * 2) a -= P; b += other.b; return *this; } Modint operator+(const Modint& other) const { return Modint(*this) += other; } Modint pow(const Modint& other) const { return {pow(a, other.b.a), b.pow(other.b)}; }; }; template<> struct Modint<1>{ uint32_t a; Modint(uint64_t x): a(bool(x)){} uint32_t val() const { return 0; } Modint& operator*=(const Modint& other){ a &= other.a; return *this; } Modint operator*(const Modint& other) const { return Modint(*this) *= other; } Modint& operator+=(const Modint& other){ a |= other.a; return *this; } Modint operator+(const Modint& other) const { return Modint(*this) += other; } Modint pow(const Modint& other) const { return {a || !other.a}; } }; using mint = Modint<1000000007>; map<int,mint> div(ll x){ map<int,mint> rt; ll cx=x; for(int i=2;i*i<=cx;i++){ if(x%i==0){ int ct=0; while(x%i==0){ ct++; x/=i; } if(!rt.count(i)){ rt.insert({i,0}); } rt.at(i)+=ct; } if(x==1) break; } if(x!=1){ rt.insert({x,1}); } return rt; } map<int,mint> f(map<int,mint> &nw){ map<int,mint> rt; for(auto [x,n]:nw){ auto mp=div(x+1); for(auto [nx,nn]:mp){ if(!rt.count(nx)){ rt.insert({nx,0}); } rt.at(nx)+=nn*n; } } return rt; } vector<vector<mint>> mmul(vector<vector<mint>> a,vector<vector<mint>> b){ vector<vector<mint>> c={{0,0},{0,0}}; rep(i,2){ rep(k,2){ rep(j,2){ c.at(i).at(j)+=a.at(i).at(k)*b.at(k).at(j); } } } return c; } int main(){ int n; ll k; cin>>n>>k; map<int,mint> nw=div(n); for(ll i=0;i<k;i++){ if(nw.size()<=2){ vector<mint> v={mint(0),mint(0)}; { int flg=0; for(auto [x,n]:nw){ if(x==2) v.at(0)=n; else if(x==3) v.at(1)=n; else{ flg=1; } } if(flg){ nw=f(nw); continue; } } vector<vector<vector<mint>>> mat(60); mat.at(0)={{0,2},{1,0}}; rep(i,59){ mat.at(i+1)=mmul(mat.at(i),mat.at(i)); } ll lp=k-i; vector<vector<mint>> vv={{1,0},{0,1}}; rep(i,60){ if((lp>>i)&1){ vv=mmul(vv,mat.at(i)); } } map<int,mint> tmp={{2,0},{3,0}}; rep(i,2){ rep(j,2){ tmp.at(i+2)+=vv.at(i).at(j)*v.at(j); } } nw=tmp; i+=lp; break; } nw=f(nw); } mint ans=1; for(auto [x,n]:nw){ ans*=(mint(x).pow(n)); } cout<<ans.val()<<endl; }