[硫化铂]高中数列题

高中数列题

题解

首先，我们的$a$的递推式显然是一个前缀和的形式，我们不妨将其全部展开，我们定义$p_i=\frac{i+b}{c},q_i=ci-b$可以得到，
$$a_n=a+\sum _{i=1}^{n}f(i)a_{p_i}=a+\sum _{i=1}^{n}f(i)(\sum _{j=1}^{p_i}f(j)a_{p_j}+a)=\sum _{i=1}^{p_n}{a}_{p_i}f(i)\sum _{j=q_i}^{n}f(j)+(1+\sum _{i=1}^{n}f(i))a$$
显然，我们只需要算一个前缀和，就可以将我们的范围从$n$缩小到$p_n$，但这显然是可以继续缩小的，我们记$sf(n)={\sum }_{i=1}^{n}f(i)$。
那么原式可以被化简，得到
$$a_n=\sum _{i=1}^{p_n}f(i)(sf(n)-sf(q_i-1))a_{p_i}+(1+sf(n))a$$
可以发现，由于$q_i$本身是一个关于$i$的一次式，所以上面的$f(q_i-1)$可以被替换成一个$m+1$次的多项式，将它与$f(i)$相乘，显然我们得到的还是一个关于$i$的多项式。
那么我们之后要计算的${\sum }_{i=1}^{p_n}f(i)(sf(n)-sf(q_i-1))a_{p_i}$显然与前面的${\sum }_{i=1}^{n}f(i)a_{p_i}$是类似的，都相当于是$a_{p_i}$做了一个关于一个多项式做了一个累和，不过两个范围不一样了。
我们可以暴力将新的多项式变成${\sum }_{i=1}^{p_n}g(i)a_{p_i}$的形式，然后继续像刚刚那样化简，继续累和。
关于${\sum }_{i=1}^{n}g(i)a_{p_i}$的化简，我们能得到的是，
$$\sum _{i=1}^{p_n}f(i)(sg(n)-sg(q_i-1))a_{p_i}+sg(n)a$$
这样就可以不断的缩小我们的范围，知道缩小到一个我们可以暴力求的范围得出答案。

事实上，对于一个$m$次项的多项式$g_i$，我们得到$sg(q_i-1)$可以通过拉格朗日插值法的方式，将$0-m+1$的$sg$都求出来，再借以算出$sg(q_i-1)$，时间复杂度是$O\left(k^2\right)$的。
至于$f(i)$乘$sg(n)-sg(q_i-1)$的计算，我们直接暴力多项式乘法就行了，反正也不会比拉格朗日插值更多。
每次化简会使得我们的次数增加$m+1$，总共${\log }_{c}n$次化简。

总时间复杂度$O\left(m^2{\log }_c^{3}n\right)$。
如果用下降幂多项式加FFT去化简的话，可以做到$O\left(m{\log }_c^{2}n{\log }^2(m{\log }_{c}n)\right)$。

源码

#include<bits/stdc++.h>
using namespace std;
#define MAXN 100005
#define lowbit(x) (x&-x)
#define reg register
#define pb push_back
#define mkpr make_pair
#define fir first
#define sec second
#define lson (rt<<1)
#define rson (rt<<1|1)
typedef long long LL;
typedef unsigned long long uLL; 
typedef long double ld;
typedef pair<int,int> pii;
const int INF=0x3f3f3f3f;
const int mo=1004535809;
const int mod=1e5+3;
const int inv2=5e8+4;
const int jzm=2333;
const int zero=2000;
const int n1=1000;
const int M=100000;
const int orG=3,ivG=332748118;
const long double Pi=acos(-1.0);
const double eps=1e-12;
template<typename _T>
_T Fabs(_T x){return x<0?-x:x;}
template<typename _T>
void read(_T &x){
	_T f=1;x=0;char s=getchar();
	while(s>'9'||s<'0'){if(s=='-')f=-1;s=getchar();}
	while('0'<=s&&s<='9'){x=(x<<3)+(x<<1)+(s^48);s=getchar();}
	x*=f;
}
template<typename _T>
void print(_T x){if(x<0){x=(~x)+1;putchar('-');}if(x>9)print(x/10);putchar(x%10+'0');}
int gcd(int a,int b){return !b?a:gcd(b,a%b);}
int add(int x,int y,int p){return x+y<p?x+y:x+y-p;}
void Add(int &x,int y,int p){x=add(x,y,p);}
int qkpow(int a,int s,int p){int t=1;while(s){if(s&1)t=1ll*t*a%p;a=1ll*a*a%p;s>>=1;}return t;}
struct poly{
	int a[1265],len;poly(){len=0;a[0]=0;}
	void clear(){for(int i=0;i<=len+1;i++)a[i]=0;}
	poly operator + (const poly &rhs)const{
		poly res;res.len=max(len,rhs.len);res.clear();int i;
		for(i=0;i<=min(rhs.len,len);i++)
			res.a[i]=add(a[i],rhs.a[i],mo);
		while(i<=len)res.a[i]=a[i],i++;
		while(i<=rhs.len)res.a[i]=rhs.a[i],i++;
		return res;
	}
	poly operator - (const poly &rhs)const{
		poly res;res.len=max(len,rhs.len);res.clear();int i;
		for(i=0;i<=min(rhs.len,len);i++)
			res.a[i]=add(a[i],mo-rhs.a[i],mo);
		while(i<=len)res.a[i]=a[i],i++;
		while(i<=rhs.len)res.a[i]=mo-rhs.a[i],i++;
		return res;
	}
	poly operator * (const poly &rhs)const{
		poly res;res.len=len+rhs.len;res.clear();
		for(int i=0;i<=len;i++)
			for(int j=0;j<=rhs.len;j++)
				Add(res.a[i+j],1ll*a[i]*rhs.a[j]%mo,mo);
		return res;
	}
	poly operator * (const int rhs)const{
		poly res;res.len=len;res.clear();
		for(int i=0;i<=len;i++)res.a[i]=1ll*rhs*a[i]%mo;
		return res;
	}
	int ask(int x){
		int res=0,now=1;
		for(int i=0;i<=len;i++,now=1ll*x*now%mo)
			Add(res,1ll*now*a[i]%mo,mo);
		return res;
	}
}P,G,F;
struct point{
	int x,y;point(){x=y=0;}
	point(int X,int Y){x=X;y=Y;}
}d[MAXN];
poly sakura(point *p,int len){
	poly res;res.len=0;res.clear();res.a[0]=p[1].y;
	poly tmp;tmp.len=1;tmp.clear();tmp.a[1]=1;tmp.a[0]=mo-p[1].x;
	for(int i=2;i<=len;i++){
		int dif=add(p[i].y,mo-res.ask(p[i].x),mo);
		int dmp=tmp.ask(p[i].x);
		int jzg=1ll*dif*qkpow(dmp,mo-2,mo)%mo;
		res=res+tmp*jzg;tmp.a[tmp.len+1]=0;
		for(int j=tmp.len;j>=0;j--)
			Add(tmp.a[j+1],tmp.a[j],mo),
			tmp.a[j]=1ll*(mo-p[i].x)*tmp.a[j]%mo;
		tmp.len++;
	}
	return res;
}
int a,b,c,m,ans,dp[MAXN];LL n;
signed main(){
	//freopen("seq.in","r",stdin);
	//freopen("seq.out","w",stdout);
	read(n);read(a);read(b);read(c);read(m);ans=0;
	P.len=m;for(int i=0;i<=m;i++)read(P.a[i]);F=P;
	dp[0]=a;for(int i=1;i<=500;i++)dp[i]=add(dp[i-1],1ll*P.ask(i)*dp[(i+b)/c]%mo,mo);
	while(n>500){
		for(int i=1;i<=P.len+2;i++)d[i].x=i-1,d[i].y=P.ask(i-1);
		for(int i=2;i<=P.len+2;i++)Add(d[i].y,d[i-1].y,mo);
		poly tmp=sakura(d,P.len+2),tp,t;tmp.a[0]=0;
		int summ=tmp.ask(n%mo);tp.len=0;tp.a[0]=summ;
		poly num;num.clear();num.len=1;num.a[1]=c;num.a[0]=mo-b-1;
		poly sum;sum.clear();sum.len=0;sum.a[0]=1;G.len=G.a[0]=0;
		for(int i=0;i<=tmp.len;i++)G=G+sum*tmp.a[i],sum=sum*num;
		G=tp-G;P=F*G;n=(n+b)/c;Add(ans,1ll*a*summ%mo,mo);
	}
	for(int i=1;i<=n;i++)Add(ans,1ll*P.ask(i)*dp[(i+b)/c]%mo,mo);
	Add(ans,a,mo);printf("%d\n",ans);
	return 0;
}

谢谢！！！