#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; typedef long long int LL; typedef pair P; typedef pair LP; const int INF=1<<30; const LL MAX=1e9+7; void array_show(int *array,int array_n,char middle=' '){ for(int i=0;i &vec_s,int vec_n=-1,char middle=' '){ if(vec_n==-1)vec_n=vec_s.size(); for(int i=0;i &vec_s,int vec_n=-1,char middle=' '){ if(vec_n==-1)vec_n=vec_s.size(); for(int i=0;i ostream& operator<<(ostream& os,const vector& v1){ int n=v1.size(); for(int i=0;i ostream& operator<<(ostream& os,const pair& p){ os< istream& operator>>(istream& is,vector& v1){ int n=v1.size(); for(int i=0;i>v1[i]; return is; } template istream& operator>>(istream& is,pair& p){ is>>p.first>>p.second; return is; } templateT ext_gcd(T a,T b,T& x,T& y){ //ax+by=gcd(a,b) if(ab)y+=(x/b)*a,x%=b; if(y>a)x+=(y/a)*b,y%=a; return gcd_val; } templateclass modint{ private: typedef long long int ll; ll val; ll gcd(ll a,ll b){ if(amodint(T a){ val=(ll)a%mod; if(val<0)val+=mod; } ll value()const{return val;} ll get_mod()const{return mod;} modint& operator++(){ val++; if(val==mod)val=0; return *this; } modint operator++(int){ modint ans=*this; ++*this; return ans; } modint& operator--(){ if(val==0)val=mod; val--; return *this; } modint operator--(int){ modint ans=*this; --*this; return ans; } modint& operator+=(const modint& a){ val+=a.value(); if(val>=mod)val-=mod; return *this; } modint& operator-=(const modint& a){ val-=a.value(); if(val<0)val+=mod; return *this; } modint& operator*=(const modint& a){ val*=a.value(); if(val>=mod)val%=mod; return *this; } modint pow(ll index)const{ assert(index>=0); if(prime && index>=mod-1)index%=mod-1; modint a=*this,ans=1; for(ll i=1;i>=0 && i<=index;i<<=1){ if(index&i)ans*=a; a*=a; } return ans; } modint inverse()const{ if(prime){ assert(val!=0); return pow(mod-2); } ll x,y; ll g=ext_gcd(val,mod,x,y); assert(g==1); return x; } modint& operator/=(const modint& a){ if(prime){ *this=(*this)*a.inverse(); return *this; } ll g=gcd(val,a.value()); modint a2=a.value()/g; val/=g; *this=(*this)*a2.inverse(); return *this; } friend modint operator-(const modint& a,const modint& b){return (modint)a-=b;} friend modint operator+(const modint& a,const modint& b){return (modint)a+=b;} friend modint operator*(const modint& a,const modint& b){return (modint)a*=b;} friend modint operator/(const modint& a,const modint& b){return (modint)a/=b;} friend bool operator==(const modint& a,const modint& b){return a.value()==b.value();} friend bool operator!=(const modint& a,const modint& b){return a.value()!=b.value();} friend modint pow(const modint& a,const ll b){return a.pow(b);} modint operator+() const{return *this;} modint operator-() const{return modint()-*this;} friend ostream& operator<<(ostream& os,const modint& a){ os<>(istream& is,modint& a){ ll val; is>>val; a=val; return is; } }; using mint=modint<1'000'000'007,true>; using modint109=modint<1'000'000'009,true>; using modint998=modint<998'244'353,true>; template class Matrix{ private: typedef long long int ll; T zero=0,e=1; vector> vec_matrix; void init(int n,int m); function same=[](T a,T b){return a==b;}; //function same=[](T a,T b){return (a-b)<=1e-9;}; Matrix gauss_and_inverse(int mode)const; public: Matrix(int n); Matrix(int n,int m); Matrix(vector>& ma); int sizeX() const; int sizeY() const; bool valid() const; T& element(int a,int b); const vector>& get_vec()const; void set_same_function(function _same){same=_same;} void E(); void E(int n); Matrix& operator += (const Matrix& mat_a); Matrix& operator -= (const Matrix& mat_a); Matrix& operator *= (const Matrix& mat_a); friend Matrix operator +(const Matrix& mat_a,const Matrix& mat_b){return (Matrix)mat_a+=mat_b;} friend Matrix operator -(const Matrix& mat_a,const Matrix& mat_b){return (Matrix)mat_a-=mat_b;} friend Matrix operator *(const Matrix& mat_a,const Matrix& mat_b){return (Matrix)mat_a*=mat_b;} Matrix operator+() const{return *this;} Matrix operator-() const{return Matrix()-*this;} Matrix& operator %= (const T mod); Matrix pow_mod(ll m,ll mod=1e9+7)const; Matrix pow(ll m)const; T tr()const; Matrix gauss()const; T det()const; Matrix inverse()const; Matrix submatrix(int a,int b,int c,int d)const;//[a,c)*[b,d) }; template void Matrix::init(int n,int m){ vec_matrix.assign(n,vector(m,zero)); } template bool Matrix::valid() const{ int n=sizeX(),m=sizeY(); for(int i=0;i Matrix::Matrix(int n){ init(n,n); } template Matrix::Matrix(int n,int m){ init(n,m); } template Matrix::Matrix(vector>& ma):vec_matrix(ma){ assert(valid()); }; template int Matrix::sizeX() const{ return vec_matrix.size(); } template int Matrix::sizeY() const{ if(vec_matrix.empty())return 0; return vec_matrix[0].size(); } template void Matrix::E(){ assert(sizeX()==sizeY()); int n=sizeX(); int i,j; for(i=0;i void Matrix::E(int n){ init(n,n); E(); } template Matrix& Matrix::operator += (const Matrix& mat){ const vector>& vec_mat=mat.get_vec; assert(sizeX()==mat.sizeX() && sizeY()==mat.sizeY()); int i,j; int n=sizeX(),m=sizeY(); for(i=0;i Matrix& Matrix::operator -= (const Matrix& mat){ const vector>& vec_mat=mat.get_vec; assert(sizeX()==mat.sizeX() && sizeY()==mat.sizeY()); int i,j; int n=sizeX(),m=sizeY(); for(i=0;i Matrix& Matrix::operator *= (const Matrix& mat){ const vector>& vec_mat=mat.get_vec(); assert(mat.sizeX()==sizeY()); int i,j,k; int n=sizeX(),m=mat.sizeY(),p=sizeY(); vector> vec_mats(n,vector(m,zero)); for(i=0;i Matrix& Matrix::operator %= (const T mod){ for(auto& vm:vec_matrix){ for(auto& num:vm){ if(num>=mod)num%=mod; } } return *this; } template vector operator *(const Matrix& mat,const vector& v1){ const vector>& vec_mat=mat.get_vec(); assert(!vec_mat.empty() && vec_mat[0].size()==v1.size()); int n=vec_mat.size(),m=vec_mat[0].size(); vector vs(n); int i,j; for(i=0;i istream& operator >> (istream& is,Matrix& a){ assert(a.valid()); int n=a.sizeX(),m=a.sizeY(); int i,j; for(i=0;i>a.element(i,j); } } return is; } template ostream& operator << (ostream& os,Matrix& a){ assert(a.valid()); int n=a.sizeX(),m=a.sizeY(); int i,j; for(i=0;i Matrix Matrix::pow_mod(ll m,ll mod)const{ assert(sizeX()==sizeY()); ll p_b=1; Matrix ms(sizeX()); ms.E(); for(;p_b<=m;p_b<<=1); for(p_b>>=1;p_b>0;p_b>>=1){ ms*=ms; ms%=mod; if(m&p_b)ms*=(*this); ms%=mod; } return ms; } template Matrix Matrix::pow(ll m)const{ assert(sizeX()==sizeY()); Matrix ms(sizeX()),ma=*this; ms.E(); for(int i=0;(1LL< T& Matrix::element(int a,int b){ assert(0<=a && a const vector>& Matrix::get_vec() const{ return vec_matrix; } template T Matrix::tr()const{ assert(sizeX()==sizeY()); int n=sizeX(); T a=e; for(int i=0;i Matrix Matrix::gauss_and_inverse(int mode)const{ //mode: 0: gauss triangle, 1: inverse int n=sizeX(),m=sizeY(); int i,j,k,l; vector> v=vec_matrix,inv(n,vector(m)); if(mode==1){ assert(n==m); for(i=0;i=0;i--){ if(v[i][i]==zero)assert(false); for(j=0;j Matrix Matrix::gauss()const{ return gauss_and_inverse(0); } template T Matrix::det()const{ assert(sizeX()==sizeY()); Matrix a=gauss(); T s=e; for(int i=0;i Matrix Matrix::inverse()const{ return gauss_and_inverse(1); } template Matrix Matrix::submatrix(int a,int b,int c,int d)const{ //[a,c)*[b,d) assert(0<=a && a s(c-a,d-b); int i,j; for(i=0;i>p>>n>>m; Matrix ma(4); ma.element(0,0)=(mint)1/p; ma.element(2,0)=(mint)(p-2)/p; ma.element(3,0)=(mint)1/p; ma.element(1,1)=(mint)1/p; ma.element(2,1)=(mint)(p-1)/p; ma.element(0,2)=1; ma.element(1,3)=1; ma=ma.pow(n); mint m1=1-ma.element(1,1); cout<<1-m1.pow(m)<