結果
問題 | 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 tatyamconstexpr 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;}