結果
| 問題 |
No.950 行列累乗
|
| コンテスト | |
| ユーザー |
maroon_kuri
|
| 提出日時 | 2019-12-14 00:29:19 |
| 言語 | C++14 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
WA
|
| 実行時間 | - |
| コード長 | 16,661 bytes |
| コンパイル時間 | 4,241 ms |
| コンパイル使用メモリ | 242,500 KB |
| 実行使用メモリ | 31,048 KB |
| 最終ジャッジ日時 | 2024-06-27 23:30:13 |
| 合計ジャッジ時間 | 8,596 ms |
|
ジャッジサーバーID (参考情報) |
judge4 / judge5 |
(要ログイン)
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 4 |
| other | AC * 56 WA * 1 |
ソースコード
#include <bits/stdc++.h>
using namespace std;
using ll=long long;
#define int ll
#define rng(i,a,b) for(int i=int(a);i<int(b);i++)
#define rep(i,b) rng(i,0,b)
#define gnr(i,a,b) for(int i=int(b)-1;i>=int(a);i--)
#define per(i,b) gnr(i,0,b)
#define pb push_back
#define eb emplace_back
#define a first
#define b second
#define bg begin()
#define ed end()
#define all(x) x.bg,x.ed
#ifdef LOCAL
#define dmp(x) cerr<<__LINE__<<" "<<#x<<" "<<x<<endl
#else
#define dmp(x) void(0)
#endif
template<class t,class u> void chmax(t&a,u b){if(a<b)a=b;}
template<class t,class u> void chmin(t&a,u b){if(b<a)a=b;}
template<class t> using vc=vector<t>;
template<class t> using vvc=vc<vc<t>>;
using pi=pair<int,int>;
using vi=vc<int>;
template<class t,class u>
ostream& operator<<(ostream& os,const pair<t,u>& p){
return os<<"{"<<p.a<<","<<p.b<<"}";
}
template<class t> ostream& operator<<(ostream& os,const vc<t>& v){
os<<"{";
for(auto e:v)os<<e<<",";
return os<<"}";
}
#define mp make_pair
#define mt make_tuple
#define one(x) memset(x,-1,sizeof(x))
#define zero(x) memset(x,0,sizeof(x))
#ifdef LOCAL
void dmpr(ostream&os){os<<endl;}
template<class T,class... Args>
void dmpr(ostream&os,const T&t,const Args&... args){
os<<t<<" ";
dmpr(os,args...);
}
#define dmp2(...) dmpr(cerr,"Line:",__LINE__,##__VA_ARGS__)
#else
#define dmp2(...) void(0)
#endif
using uint=unsigned;
using ull=unsigned long long;
template<class t,size_t n>
ostream& operator<<(ostream&os,const array<t,n>&a){
return os<<vc<t>(all(a));
}
template<int i,class T>
void print_tuple(ostream&,const T&){
}
template<int i,class T,class H,class ...Args>
void print_tuple(ostream&os,const T&t){
if(i)os<<",";
os<<get<i>(t);
print_tuple<i+1,T,Args...>(os,t);
}
template<class ...Args>
ostream& operator<<(ostream&os,const tuple<Args...>&t){
os<<"{";
print_tuple<0,tuple<Args...>,Args...>(os,t);
return os<<"}";
}
void print(ll x,int suc=1){
cout<<x;
if(suc==1)
cout<<"\n";
if(suc==2)
cout<<" ";
}
ll read(){
ll i;
cin>>i;
return i;
}
vi readvi(int n,int off=0){
vi v(n);
rep(i,n)v[i]=read()+off;
return v;
}
template<class T>
void print(const vector<T>&v,int suc=1){
rep(i,v.size())
print(v[i],i==int(v.size())-1?suc:2);
}
string readString(){
string s;
cin>>s;
return s;
}
template<class T>
T sq(const T& t){
return t*t;
}
//#define CAPITAL
void yes(bool ex=true){
#ifdef CAPITAL
cout<<"YES"<<endl;
#else
cout<<"Yes"<<endl;
#endif
if(ex)exit(0);
}
void no(bool ex=true){
#ifdef CAPITAL
cout<<"NO"<<endl;
#else
cout<<"No"<<endl;
#endif
if(ex)exit(0);
}
void possible(bool ex=true){
#ifdef CAPITAL
cout<<"POSSIBLE"<<endl;
#else
cout<<"Possible"<<endl;
#endif
if(ex)exit(0);
}
void impossible(bool ex=true){
#ifdef CAPITAL
cout<<"IMPOSSIBLE"<<endl;
#else
cout<<"Impossible"<<endl;
#endif
if(ex)exit(0);
}
constexpr ll ten(int n){
return n==0?1:ten(n-1)*10;
}
const ll infLL=LLONG_MAX/3;
#ifdef int
const int inf=infLL;
#else
const int inf=INT_MAX/2-100;
#endif
int topbit(signed t){
return t==0?-1:31-__builtin_clz(t);
}
int topbit(ll t){
return t==0?-1:63-__builtin_clzll(t);
}
int botbit(signed a){
return a==0?32:__builtin_ctz(a);
}
int botbit(ll a){
return a==0?64:__builtin_ctzll(a);
}
int popcount(signed t){
return __builtin_popcount(t);
}
int popcount(ll t){
return __builtin_popcountll(t);
}
bool ispow2(int i){
return i&&(i&-i)==i;
}
int mask(int i){
return (int(1)<<i)-1;
}
bool inc(int a,int b,int c){
return a<=b&&b<=c;
}
template<class t> void mkuni(vc<t>&v){
sort(all(v));
v.erase(unique(all(v)),v.ed);
}
ll rand_int(ll l, ll r) { //[l, r]
#ifdef LOCAL
static mt19937_64 gen;
#else
static random_device rd;
static mt19937_64 gen(rd());
#endif
return uniform_int_distribution<ll>(l, r)(gen);
}
template<class t>
int lwb(const vc<t>&v,const t&a){
return lower_bound(all(v),a)-v.bg;
}
using uint=unsigned;
using ull=unsigned long long;
//initfact();
//const uint mod=998244353;
//const uint mod=1000000007;
uint mod=1;
struct mint{
uint v;
mint(ll vv=0){s(vv%mod+mod);}
mint& s(uint vv){
v=vv<mod?vv:vv-mod;
return *this;
}
mint operator-()const{return mint()-*this;}
mint& operator+=(const mint&rhs){return s(v+rhs.v);}
mint&operator-=(const mint&rhs){return s(v+mod-rhs.v);}
mint&operator*=(const mint&rhs){
v=ull(v)*rhs.v%mod;
return *this;
}
mint&operator/=(const mint&rhs){return *this*=rhs.inv();}
mint operator+(const mint&rhs)const{return mint(*this)+=rhs;}
mint operator-(const mint&rhs)const{return mint(*this)-=rhs;}
mint operator*(const mint&rhs)const{return mint(*this)*=rhs;}
mint operator/(const mint&rhs)const{return mint(*this)/=rhs;}
mint pow(int n)const{
mint res(1),x(*this);
while(n){
if(n&1)res*=x;
x*=x;
n>>=1;
}
return res;
}
mint inv()const{return pow(mod-2);}
/*mint inv()const{
int x,y;
int g=extgcd(v,mod,x,y);
assert(g==1);
if(x<0)x+=mod;
return mint(x);
}*/
friend mint operator+(int x,const mint&y){
return mint(x)+y;
}
friend mint operator-(int x,const mint&y){
return mint(x)-y;
}
friend mint operator*(int x,const mint&y){
return mint(x)*y;
}
friend mint operator/(int x,const mint&y){
return mint(x)/y;
}
friend ostream& operator<<(ostream&os,const mint&m){
return os<<m.v;
}
bool operator<(const mint&r)const{return v<r.v;}
bool operator==(const mint&r)const{return v==r.v;}
bool operator!=(const mint&r)const{return v!=r.v;}
explicit operator bool()const{
return v;
}
};
const int vmax=(1<<21)+10;
mint fact[vmax],finv[vmax],invs[vmax];
void initfact(){
fact[0]=1;
rng(i,1,vmax){
fact[i]=fact[i-1]*i;
}
finv[vmax-1]=fact[vmax-1].inv();
for(int i=vmax-2;i>=0;i--){
finv[i]=finv[i+1]*(i+1);
}
for(int i=vmax-1;i>=1;i--){
invs[i]=finv[i]*fact[i-1];
}
}
mint choose(int n,int k){
return fact[n]*finv[n-k]*finv[k];
}
mint binom(int a,int b){
return fact[a+b]*finv[a]*finv[b];
}
mint catalan(int n){
return binom(n,n)-(n-1>=0?binom(n-1,n+1):0);
}
template<class num>
struct Vector{
vector<num> v;
Vector(int s=0){
v.resize(s,num(0));
}
int size()const{
return v.size();
}
num& operator[](int i){
return v[i];
}
num const& operator[](int i)const{
return v[i];
}
Vector& operator+=(const Vector&rhs){
assert(size()==rhs.size());
rep(i,size())
v[i]+=rhs[i];
return *this;
}
Vector& operator-=(const Vector&rhs){
assert(size()==rhs.size());
rep(i,size())
v[i]-=rhs[i];
return *this;
}
Vector& operator*=(const num&x){
rep(i,size())
v[i]*=x;
return *this;
}
Vector& operator/=(const num&x){
num y=num(1)/x;
rep(i,size())
v[i]*=y;
return *this;
}
Vector operator+(const Vector&rhs)const{
return Vector(*this)+=rhs;
}
Vector operator-(const Vector&rhs)const{
return Vector(*this)-=rhs;
}
Vector operator*(const num&x)const{
return Vector(*this)*=x;
}
Vector operator/(const num&x)const{
return Vector(*this)/=x;
}
bool operator==(const Vector&rhs)const{
return v==rhs.v;
}
};
template<class num>
num dot(const Vector<num>&a,const Vector<num>&b){
assert(a.size()==b.size());
const int s=a.size();
num ans(0);
rep(i,s)
ans+=a[i]*b[i];
return ans;
}
template<class num>
ostream&operator<<(ostream&os,const Vector<num>&v){
return os<<v.v;
}
template<class num>
struct Matrix{
using V=Vector<num>;
vector<V> m;
Matrix(int s=0,num z=num(0)){
m.resize(s,V(s));
rep(i,size())
m[i][i]=z;
}
int size()const{
return m.size();
}
Matrix operator*(const Matrix&rhs)const{
assert(size()==rhs.size());
Matrix res(size());
rep(i,size())rep(j,size())rep(k,size())
res[i][j]+=m[i][k]*rhs[k][j];
return res;
}
Matrix& operator*=(const Matrix&rhs){
return *this=(*this)*rhs;
}
V operator*(const V&rhs)const{
assert(size()==rhs.size());
V res(size());
rep(i,size())
res[i]=dot(m[i],rhs);
return res;
}
V& operator[](int i){
return m[i];
}
V const& operator[](int i) const{
return m[i];
}
Matrix& operator/=(const num&r){
rep(i,m.size())m[i]/=r;
return *this;
}
Matrix operator/(const num&r)const{
return Matrix(*this)/=r;
}
bool operator==(const Matrix&rhs)const{
return m==rhs.m;
}
};
template<class num>
ostream&operator<<(ostream&os,const Matrix<num>&x){
const int s=x.size();
os<<"----------"<<endl;
rep(i,s){
rep(j,s){
os<<x[i][j];
if(j==s-1){
os<<endl;
}else{
os<<" ";
}
}
}
return os<<"----------";
}
mint phi;
struct F{
mint a,b;
F(mint aa=0,mint bb=0):a(aa),b(bb){}
F operator-()const{
return F(-a,-b);
}
F& operator+=(const F&r){
a+=r.a;
b+=r.b;
return *this;
}
F operator+(const F&r)const{
return F(*this)+=r;
}
F& operator-=(const F&r){
a-=r.a;
b-=r.b;
return *this;
}
F operator-(const F&r)const{
return F(*this)-=r;
}
F operator*(const F&r)const{
return F(a*r.a+b*r.b*phi,a*r.b+b*r.a);
}
F& operator*=(const F&r){
return *this=*this*r;
}
F& operator/=(const F&r){
return *this*=r.inv();
}
F operator/(const F&r)const{
return F(*this)/=r;
}
F pow(int n)const{
F res{1,0},x=*this;
while(n){
if(n&1)res=res*x;
x=x*x;
n>>=1;
}
return res;
}
F inv()const{
return pow(ll(mod)*mod-2);
}
bool operator<(const F&r)const{
if(a!=r.a)return a<r.a;
else return b<r.b;
}
bool operator==(const F&r)const{
return a==r.a&&b==r.b;
}
bool operator!=(const F&r)const{
return a!=r.a||b!=r.b;
}
bool is_zero()const{
return a.v==0&&b.v==0;
}
explicit operator bool()const{
return a||b;
}
};
ostream&operator<<(ostream&os,const F&f){
return os<<"{"<<f.a<<","<<f.b<<"}";
}
vi factors(int x){
vi res;
for(int i=2;i*i<=x;i+=(i==2?1:2)){
if(x%i==0){
res.pb(i);
while(x%i==0)
x/=i;
}
}
if(x>1)res.pb(x);
sort(all(res));
return res;
}
//require a^s=1
//minimum n s.t. a^n=b, 0<=n<=s-1
//or -1 if infeasible
template<class t>
int discrete_log(t a,t b,int s){
dmp2(a,b,s);
assert(a.pow(s)==t(1));
const int L=ceil(sqrt(s));
int ans=inf;
map<t,int> v;
{
t x=a.pow(L);
t cur(1);
for(int i=0;i<s;i+=L){
if(!v.count(cur))v[cur]=i;
cur*=x;
}
}
t x=a.inv();
dmp(x);
t cur(1);
rep(i,L){
t y=b*cur;
if(v.count(y))
chmin(ans,v[y]+i);
cur*=x;
}
if(ans>=s)
ans=-1;
else{
dmp(L);
dmp(ans);
assert(inc(0,ans,s-1));
assert(a.pow(ans)==b);
}
return ans;
}
//require a^s=1
//returns minimum p s.t. a^p=1
template<class t>
int mult_period(t a,int s){
vi f=factors(s);
int p=s;
for(auto v:f){
while(p%v==0){
if(a.pow(p/v)==t(1))
p/=v;
else
break;
}
}
return p;
}
//require b^s=1
//pair of
//minimum n s.t. a*b^n=c,0<=n<=s-1
//and period
//or (-1,-1) if infeasible
template<class t>
pi discrete_log_helper(t a,t b,t c,int s){
dmp2(a,b,c,s);
if(!a||!b){
if(!c)return pi(0,1);
else return pi(-1,-1);
}
assert(b.pow(s)==t(1));
int w=discrete_log(b,c*a.inv(),s);
if(w==-1)return pi(-1,-1);
return pi(w,mult_period(b,s));
}
int norm_mod(int a,int p){
a%=p;
if(a<0)a+=p;
return a;
}
//p: odd (not necessarily prime)
//10^18 サイズの入力で動くはず(未検証)
int jacobi(int a,int p){
a=norm_mod(a,p);
auto neg=[](int x){return x%2?-1:1;};
if(a==0) return p==1;
else if(a%2) return neg(((p-1)&(a-1))>>1)*jacobi(p%a,a);
else return neg(((p&15)*(p&15)-1)/8)*jacobi(a/2,p);
}
//yosupo sqrt mod
//LOJ 150
//p: prime
//0<=a<p
//returns minimum solution
int sqrt_mod(int a,int p){
if(p==2)return a%2;
if(a==0)return 0;
if(jacobi(a,p)==-1)return -1;
int b,c;
do{
b=rand_int(0,p-1);
c=norm_mod(b*b-a,p);
}while(jacobi(c,p)!=-1);
auto mul=[&](pi x,pi y){
return pi(norm_mod(x.a*y.a+x.b*y.b%p*c,p),norm_mod(x.a*y.b+x.b*y.a,p));
};
pi x(b,1),cur(1,0);
int n=(p+1)/2;
while(n){
if(n&1)cur=mul(cur,x);
x=mul(x,x);
n>>=1;
}
assert(cur.b==0);
return min(cur.a,p-cur.a);
}
int extgcd(int a,int b,int&x,int&y){
if(b==0){
x=1;
y=0;
return a;
}else{
int g=extgcd(b,a%b,y,x);
y-=a/b*x;
return g;
}
}
//x*y=g mod m
//returns (g,y)
pi modinv(int x,int m){
int a,b;
int g=extgcd(x,m,a,b);
if(a<0)a+=m/g;
return pi(g,a);
}
//x = a mod b
//x = c mod d
// => x = p mod q
//returns (p,q) or (-1,-1) if infeasible
pi crt(int a,int b,int c,int d){
if(b<d){
swap(a,c);
swap(b,d);
}
c%=d;
int g,e;
tie(g,e)=modinv(b,d);
int k=(c-a);
if(k%g)return pi(-1,-1);
k/=g;
int m=b/g*d;
int x=(a+k*e%(d/g)*b)%m;
if(x<0)x+=m;
return pi(x,m);
}
/*
void imp(){
cout<<-1<<endl;
exit(0);
}
void answer(int a){
cout<<a<<endl;
exit(0);
}
*/
template<class t>
t det22(const Matrix<t>&a){
return a[0][0]*a[1][1]-a[0][1]*a[1][0];
}
template<class t>
Matrix<t> inv22(const Matrix<t>&a){
Matrix<t> res(2);
res[0][0]=a[1][1];
res[1][0]=-a[1][0];
res[0][1]=-a[0][1];
res[1][1]=a[0][0];
dmp(res*a);
t d=det22(a);
assert(d);
return res/d;
}
pi mergeans(pi a,pi b){
if(a==pi(-1,-1)||b==pi(-1,-1))return pi(-1,-1);
return crt(a.a,a.b,b.a,b.b);
}
int normans(pi a){
if(a==pi(-1,-1))return -1;
int w=a.a;
if(w==0)w+=a.b;
return w;
}
template<class t>
t general_pow(t a,int n){
if(n>=0)return a.pow(n);
else return a.pow(-n).inv();
}
//solve the problem given two distinct eigenvalues
template<class t>
int solve_sub(Matrix<mint> aa,Matrix<mint> bb,t p,t q,int relka){
Matrix<t> a(2);
rep(i,2)rep(j,2)
a[i][j]=aa[i][j];
Matrix<t> b(2);
rep(i,2)rep(j,2)
b[i][j]=bb[i][j];
auto getvec=[&](t w){
Vector<t> res(2);
if(a[0][0]!=w){
res[0]=-a[0][1];
res[1]=a[0][0]-w;
}else{
res[0]=a[1][1]-w;
res[1]=-a[1][0];
}
return res;
};
Vector<t> v[2];
v[0]=getvec(p);
v[1]=getvec(q);
rep(i,2){
dmp(v[i]);
assert(a*v[i]==v[i]*(i==0?p:q));
}
Matrix<t> w(2);
rep(i,2)rep(j,2)
w[i][j]=v[j][i];
dmp(w);
Matrix<t> winv=inv22(w);
dmp(winv);
dmp(w*winv);
t ad=det22(a);
assert(ad==p*q);
Matrix<t> tar=winv*b*w;
dmp(tar);
if(tar[1][0]||tar[0][1])return -1;
t bd=tar[0][0]*tar[1][1];
int off,z;
tie(off,z)=discrete_log_helper(t(1),ad,bd,mod-1);
if(off==-1)return -1;
auto waf=[&](t g,t h){
h/=g.pow(off);
g=g.pow(z);
int x,per;
tie(x,per)=discrete_log_helper(t(1),g,h,relka);
if(x==-1)return pi(-1,-1);
return pi(off+z*x,z*per);
};
return normans(mergeans(
waf(p,tar[0][0]),
waf(q,tar[1][1])
));
}
int solve_simple(Matrix<mint> a,Matrix<mint> b){
assert(a[0][1]==0);
assert(a[1][0]==0);
if(b[0][1]||b[1][0])return -1;
return normans(mergeans(
discrete_log_helper(mint(1),a[0][0],b[0][0],mod-1),
discrete_log_helper(mint(1),a[1][1],b[1][1],mod-1)
));
}
int solve(Matrix<mint> a,Matrix<mint> b,int brutethreshold=2){
phi=0;
if(mod<=brutethreshold){
Matrix<mint> cur(2,1);
rep(i,mod*mod){
if(i){
if(cur==b){
return i;
}
}
cur*=a;
}
return -1;
}
mint c=-a[0][0]-a[1][1];
mint d=a[0][0]*a[1][1]-a[0][1]*a[1][0];
dmp2(c,d);
if(d==0){
if(a==b)return 1;
pi ans(0,1);
rep(i,2)rep(j,2){
pi x=discrete_log_helper(a[i][j],-c,b[i][j],mod-1);
ans=mergeans(x,ans);
}
int w=normans(ans);
if(w>=0)w++;
return w;
}else{
mint dis=c*c/4-d;
dmp(dis);
if(dis==0){
mint alpha=-c/2;
Matrix<mint> apri=a;
apri[0][0]-=alpha;
apri[1][1]-=alpha;
dmp(alpha);
Vector<mint> v(2);
if(apri[0][0]||apri[1][0])
v[0]=1;
else if(apri[0][1]||apri[1][1])
v[1]=1;
else
return solve_simple(a,b);
Vector<mint> u=apri*v;
Matrix<mint> w(2);
w[0][0]=u[0];
w[1][0]=u[1];
w[0][1]=v[0];
w[1][1]=v[1];
auto winv=inv22(w);
auto tar=winv*b*w;
if(tar[1][0])return -1;
if(tar[0][0]!=tar[1][1])return -1;
pi k=discrete_log_helper(mint(1),alpha,tar[0][0],mod-1);
if(k==pi(-1,-1))return -1;
mint z=tar[0][1]/general_pow(alpha,k.a-1);
mint x=(z-k.a)/k.b;
return k.a+x.v*k.b;
}else{
int s=sqrt_mod(dis.v,mod);
dmp(s);
if(s!=-1){
return solve_sub(a,b,-c/2+s,-c/2-s,mod-1);
}else{
phi=dis;
return solve_sub(a,b,F(-c/2,1),F(-c/2,-1),mod+1);
}
}
}
}
bool isprime(int p){
for(int i=2;i*i<=p;i+=(i==2?1:2))
if(p%i==0)return false;
return true;
}
signed main(){
cin.tie(0);
ios::sync_with_stdio(0);
cout<<fixed<<setprecision(20);
int n;cin>>n;
if(n>0){
mod=n;
Matrix<mint> a(2),b(2);
rep(i,2)rep(j,2)cin>>a[i][j].v;
rep(i,2)rep(j,2)cin>>b[i][j].v;
cout<<solve(a,b)<<endl;
}else{
rep(_,100000){
do{
mod=rand_int(100,200);
}while(!isprime(mod));
Matrix<mint> a(2),b(2);
rep(i,2)rep(j,2)
a[i][j]=rand_int(0,mod-1);
rep(i,2)rep(j,2)
b[i][j]=rand_int(0,mod-1);
{
cout<<mod<<endl;
cout<<a<<endl;
cout<<b<<endl;
}
int x=solve(a,b);
int y=solve(a,b,inf);
if(x!=y){
cout<<mod<<endl;
cout<<a<<endl;
cout<<b<<endl;
exit(1);
}
}
}
}
maroon_kuri