library
This documentation is automatically generated by online-judge-tools/verification-helper
P-recursive
(FPS/p-recursive.hpp)
- View this file on GitHub
- Last update: 2024-01-12 04:16:01+09:00
- Include:
#include "FPS/p-recursive.hpp"
Depends on
Shift of Sampling Points of Polynomial
(FPS/samplepointshift.hpp)
- Linear Equation
(Math/linearequation.hpp)
- Matrix
(Math/matrix.hpp)
Code
#pragma once
#include "FPS/samplepointshift.hpp"
#include "Math/matrix.hpp"
#include "Math/linearequation.hpp"
template<typename T>Matrix<T> PrefixProdOfPolyMatrix(Matrix<Poly<T>>& m,ll K){
using Mat=Matrix<T>;
int n=m.val.size();
int deg=1;
rep(i,0,n)rep(j,0,n)chmax(deg,(int)m[i][j].size()-1);
ll SQ=1;
while(SQ*SQ*deg<K)SQ<<=1;
T iSQ=T(SQ).inv();
vector<Mat> G(deg+1);
rep(k,0,deg+1){
G[k]=Mat(n,n);
rep(i,0,n)rep(j,0,n)G[k][i][j]=m[i][j].eval(SQ*k);
}
auto process=[&](vector<Mat>& base,T x)->vector<Mat>{
int D=base.size();
vector ret(D,Mat(n,n));
rep(i,0,n)rep(j,0,n){
vector<T> val(D);
rep(k,0,D)val[k]=base[k][i][j];
auto add=SamplePointsShift<T>(val,x);
rep(k,0,D)ret[k][i][j]=add[k];
}
return ret;
};
for(ll w=1;w<SQ;w<<=1){
auto G1=process(G,iSQ*w);
auto G2=process(G,w*deg+1);
auto G3=process(G,iSQ*w+w*deg+1);
rep(i,0,w*deg+1)G1[i]*=G[i],G3[i]*=G2[i];
G1.insert(G1.end(),ALL(G3));
G1.pop_back();
swap(G,G1);
}
Mat ret(n,n);
rep(i,0,n)ret[i][i]=1;
ll k=0;
while(k*SQ+SQ<=K)ret=G[k++]*ret;
k*=SQ;
while(k<K){
Mat mul(n,n);
rep(i,0,n)rep(j,0,n)mul[i][j]=m[i][j].eval(k);
ret=mul*ret;
k++;
}
return ret;
}
// a_{n+i}*f_n(i)+...+a_i*f_0(i)=0
// {f_r}:dec order!!!
template<typename T>vector<Poly<T>> FindPRecursive(vector<T>& a,int d){
int n=a.size();
int k=(n+2)/(d+2)-1;
if(k<=0)return {};
int m=(d+1)*(k+1);
Matrix<T> mat(m-1,m);
rep(i,0,m-1)rep(j,0,k+1){
T base=1;
rep(deg,0,d+1){
mat[i][(d+1)*j+deg]=a[i+j]*base;
base*=(i+j);
}
}
auto basis=LinearEquation(mat,vector<T>(m-1)).second;
if(basis.val.empty())return {};
auto c=basis[0];
vector<Poly<T>> ret;
for(int i=0;i*(d+1)<(int)c.size();i++){
Poly<T> add,base({T(i),T(1)});
for(int j=d;j>=0;j--){
add*=base;
if(c[i*(d+1)+j]!=0)add+=c[i*(d+1)+j];
}
ret.push_back(add);
}
while(ret.back().empty())ret.pop_back();
reverse(ALL(ret));
return ret;
}
template<typename T>T KthtermOfPRecursive(vector<T>& a,vector<Poly<T>>& fs,ll k){
int n=fs.size()-1;
assert(int(a.size())>=n);
if(k<int(a.size()))return a[k];
Matrix<Poly<T>> m(n),den(1);
Matrix<T> base(n);
rep(i,0,n)m[0][i]=-fs[i+1];
rep(i,1,n)m[i][i-1]=fs[0];
den[0][0]=fs[0];
rep(i,0,n)base[i][0]=a[n-1-i];
T ret=(PrefixProdOfPolyMatrix(m,k-n+1)*base)[0][0];
ret/=PrefixProdOfPolyMatrix(den,k-n+1)[0][0];
return ret;
}
template<typename T>T KthtermEsper(vector<T>& a,ll k){
if(k<(int)a.size())return a[k];
int n=a.size()-1;
vector<Fp> b=a;
b.pop_back();
for(int d=0;;d++){
if((n+2)/(d+2)<=1)break;
auto fs=FindPRecursive(b,d);
if(KthtermOfPRecursive(b,fs,n)==a.back()){
return KthtermOfPRecursive(a,fs,k);
}
}
cerr<<"esper Failed"<<'\n';
assert(0);
}
/**
* @brief P-recursive
*/#line 2 "FPS/samplepointshift.hpp"
template<typename T>Poly<T> SamplePointsShift(vector<T>& ys,T c,int m=-1){
ll n=ys.size()-1,C=c.v%T::get_mod();
if(m==-1)m=n+1;
if(C<=n){
Poly<T> res;
rep(i,C,n+1)res.push_back(ys[i]);
if(int(res.size())>=m){
res.resize(m);
return res;
}
auto add=SamplePointsShift<T>(ys,n+1,m-res.size());
for(int i=0;int(res.size())<m;i++){
res.push_back(add[i]);
}
return res;
}
if(C+m>T::get_mod()){
auto res=SamplePointsShift<T>(ys,c,T::get_mod()-c.v);
auto add=SamplePointsShift<T>(ys,0,m-res.size());
rep(i,0,add.size())res.push_back(add[i]);
return res;
}
Poly<T> A(n+1),B(m+n);
rep(i,0,n+1){
A[i]=ys[i]*Fact<T>(i,1)*Fact<T>(n-i,1);
if((n-i)&1)A[i]=-A[i];
}
rep(i,0,m+n)B[i]=Fp(1)/(c-n+i);
auto AB=A*B;
vector<T> res(m);
Fp base=1;
rep(x,0,n+1)base*=(c-x);
rep(i,0,m){
res[i]=AB[n+i]*base;
base*=(c+i+1);
base*=B[i];
}
return res;
}
/**
* @brief Shift of Sampling Points of Polynomial
*/
#line 2 "Math/matrix.hpp"
template<class T>struct Matrix{
int h,w; vector<vector<T>> val; T det;
Matrix(){}
Matrix(int n):h(n),w(n),val(vector<vector<T>>(n,vector<T>(n))){}
Matrix(int n,int m):h(n),w(m),val(vector<vector<T>>(n,vector<T>(m))){}
vector<T>& operator[](const int i){return val[i];}
Matrix& operator+=(const Matrix& m){
assert(h==m.h and w==m.w);
rep(i,0,h)rep(j,0,w)val[i][j]+=m.val[i][j];
return *this;
}
Matrix& operator-=(const Matrix& m){
assert(h==m.h and w==m.w);
rep(i,0,h)rep(j,0,w)val[i][j]-=m.val[i][j];
return *this;
}
Matrix& operator*=(const Matrix& m){
assert(w==m.h);
Matrix<T> res(h,m.w);
rep(i,0,h)rep(j,0,m.w)rep(k,0,w)res.val[i][j]+=val[i][k]*m.val[k][j];
*this=res; return *this;
}
Matrix operator+(const Matrix& m)const{return Matrix(*this)+=m;}
Matrix operator-(const Matrix& m)const{return Matrix(*this)-=m;}
Matrix operator*(const Matrix& m)const{return Matrix(*this)*=m;}
Matrix pow(ll k){
Matrix<T> res(h,h),c=*this; rep(i,0,h)res.val[i][i]=1;
while(k){if(k&1)res*=c; c*=c; k>>=1;} return res;
}
vector<int> gauss(int c=-1){
if(val.empty())return {};
if(c==-1)c=w;
int cur=0; vector<int> res; det=1;
rep(i,0,c){
if(cur==h)break;
rep(j,cur,h)if(val[j][i]!=0){
swap(val[cur],val[j]);
if(cur!=j)det*=-1;
break;
}
det*=val[cur][i];
if(val[cur][i]==0)continue;
rep(j,0,h)if(j!=cur){
T z=val[j][i]/val[cur][i];
rep(k,i,w)val[j][k]-=val[cur][k]*z;
}
res.push_back(i);
cur++;
}
return res;
}
Matrix inv(){
assert(h==w);
Matrix base(h,h*2),res(h,h);
rep(i,0,h)rep(j,0,h)base[i][j]=val[i][j];
rep(i,0,h)base[i][h+i]=1;
base.gauss(h);
det=base.det;
rep(i,0,h)rep(j,0,h)res[i][j]=base[i][h+j]/base[i][i];
return res;
}
bool operator==(const Matrix& m){
assert(h==m.h and w==m.w);
rep(i,0,h)rep(j,0,w)if(val[i][j]!=m.val[i][j])return false;
return true;
}
bool operator!=(const Matrix& m){
assert(h==m.h and w==m.w);
rep(i,0,h)rep(j,0,w)if(val[i][j]==m.val[i][j])return false;
return true;
}
friend istream& operator>>(istream& is,Matrix& m){
rep(i,0,m.h)rep(j,0,m.w)is>>m[i][j];
return is;
}
friend ostream& operator<<(ostream& os,Matrix& m){
rep(i,0,m.h){
rep(j,0,m.w)os<<m[i][j]<<(j==m.w-1 and i!=m.h-1?'\n':' ');
}
return os;
}
};
/**
* @brief Matrix
*/
#line 3 "Math/linearequation.hpp"
template<typename T>pair<vector<T>,Matrix<T>> LinearEquation(Matrix<T> a,vector<T> b){
int h=a.h,w=a.w;
rep(i,0,h)a[i].push_back(b[i]);
a.w++;
vector<int> idx=a.gauss(w);
rep(i,idx.size(),h)if(a[i][w]!=0)return {{},{}};
vector<T> res(w);
rep(i,0,idx.size())res[idx[i]]=a[i][w]/a[i][idx[i]];
Matrix<T> d(w,h+w);
rep(i,0,h)rep(j,0,w)d[j][i]=a[i][j];
rep(i,0,w)d[i][h+i]=1;
int r=d.gauss(h).size();
Matrix<T> basis(w-r,w);
rep(i,r,w)basis[i-r]={d[i].begin()+h,d[i].end()};
return {res,basis};
}
/**
* @brief Linear Equation
*/
#line 5 "FPS/p-recursive.hpp"
template<typename T>Matrix<T> PrefixProdOfPolyMatrix(Matrix<Poly<T>>& m,ll K){
using Mat=Matrix<T>;
int n=m.val.size();
int deg=1;
rep(i,0,n)rep(j,0,n)chmax(deg,(int)m[i][j].size()-1);
ll SQ=1;
while(SQ*SQ*deg<K)SQ<<=1;
T iSQ=T(SQ).inv();
vector<Mat> G(deg+1);
rep(k,0,deg+1){
G[k]=Mat(n,n);
rep(i,0,n)rep(j,0,n)G[k][i][j]=m[i][j].eval(SQ*k);
}
auto process=[&](vector<Mat>& base,T x)->vector<Mat>{
int D=base.size();
vector ret(D,Mat(n,n));
rep(i,0,n)rep(j,0,n){
vector<T> val(D);
rep(k,0,D)val[k]=base[k][i][j];
auto add=SamplePointsShift<T>(val,x);
rep(k,0,D)ret[k][i][j]=add[k];
}
return ret;
};
for(ll w=1;w<SQ;w<<=1){
auto G1=process(G,iSQ*w);
auto G2=process(G,w*deg+1);
auto G3=process(G,iSQ*w+w*deg+1);
rep(i,0,w*deg+1)G1[i]*=G[i],G3[i]*=G2[i];
G1.insert(G1.end(),ALL(G3));
G1.pop_back();
swap(G,G1);
}
Mat ret(n,n);
rep(i,0,n)ret[i][i]=1;
ll k=0;
while(k*SQ+SQ<=K)ret=G[k++]*ret;
k*=SQ;
while(k<K){
Mat mul(n,n);
rep(i,0,n)rep(j,0,n)mul[i][j]=m[i][j].eval(k);
ret=mul*ret;
k++;
}
return ret;
}
// a_{n+i}*f_n(i)+...+a_i*f_0(i)=0
// {f_r}:dec order!!!
template<typename T>vector<Poly<T>> FindPRecursive(vector<T>& a,int d){
int n=a.size();
int k=(n+2)/(d+2)-1;
if(k<=0)return {};
int m=(d+1)*(k+1);
Matrix<T> mat(m-1,m);
rep(i,0,m-1)rep(j,0,k+1){
T base=1;
rep(deg,0,d+1){
mat[i][(d+1)*j+deg]=a[i+j]*base;
base*=(i+j);
}
}
auto basis=LinearEquation(mat,vector<T>(m-1)).second;
if(basis.val.empty())return {};
auto c=basis[0];
vector<Poly<T>> ret;
for(int i=0;i*(d+1)<(int)c.size();i++){
Poly<T> add,base({T(i),T(1)});
for(int j=d;j>=0;j--){
add*=base;
if(c[i*(d+1)+j]!=0)add+=c[i*(d+1)+j];
}
ret.push_back(add);
}
while(ret.back().empty())ret.pop_back();
reverse(ALL(ret));
return ret;
}
template<typename T>T KthtermOfPRecursive(vector<T>& a,vector<Poly<T>>& fs,ll k){
int n=fs.size()-1;
assert(int(a.size())>=n);
if(k<int(a.size()))return a[k];
Matrix<Poly<T>> m(n),den(1);
Matrix<T> base(n);
rep(i,0,n)m[0][i]=-fs[i+1];
rep(i,1,n)m[i][i-1]=fs[0];
den[0][0]=fs[0];
rep(i,0,n)base[i][0]=a[n-1-i];
T ret=(PrefixProdOfPolyMatrix(m,k-n+1)*base)[0][0];
ret/=PrefixProdOfPolyMatrix(den,k-n+1)[0][0];
return ret;
}
template<typename T>T KthtermEsper(vector<T>& a,ll k){
if(k<(int)a.size())return a[k];
int n=a.size()-1;
vector<Fp> b=a;
b.pop_back();
for(int d=0;;d++){
if((n+2)/(d+2)<=1)break;
auto fs=FindPRecursive(b,d);
if(KthtermOfPRecursive(b,fs,n)==a.back()){
return KthtermOfPRecursive(a,fs,k);
}
}
cerr<<"esper Failed"<<'\n';
assert(0);
}
/**
* @brief P-recursive
*/