HDU.5628.Clarke and math(狄利克雷卷积 快速幂)

(Description)

  $$g(i)=sum_{i_1|i}sum_{i_2|i_1}sum_{i_3|i_2}cdotssum_{i_k|i_{k-1}}f(i_k) mod 1000000007$$
  给出(n,k,f[1sim n]),求(g[1sim n]).

(Solution)

  首先狄利克雷卷积(Dirichlet Product):设(f(n),g(n))是两个数论函数,它们的Dirichlet乘积也是一个数论函数,

[h(n)=sum_{d|n}f(d)g(frac{n}{d}) ]

  简记为(h(n)=f(n)*g(n))

狄利克雷卷积有几个性质:
  1. 满足交换律 (f*g=g*f)
  2. 满足结合律 ((f*g)*h=f*(g*h))
  3. 满足分配率 (f*(g+h)=f*g+f*h)
  4. 存在单位元(e),使得(e*f=f*e=f)

  回到本题。设(I(x)=1).
  将式子依次展开
  $$f'(i_{k-1})=sum_{i_k|i_{k-1}}f(i_k)=sum_{i_k|i_{k-1}}f(i_k)I(frac{i_{k-1}}{i_k}) , 即f'=fI$$
  $$f''(i_{k-2})=sum_{i_{k-1}|i_{k-2}}f'(i_k-1)=sum_{i_{k-1}|i_{k-2}}f'(i_k-1)I(frac{i_{k-2}}{i_{k-1}}) , 即f''=f'*I$$.
  (ldots)
  这样下去可以得到(g=I*I*I*cdots*I*f(k个I))。由于狄利克雷卷积满足结合律,所以(k个I)的狄利克雷卷积可以用快速幂(logk)计算。
  计算狄利克雷卷积时,如果对每个(g(i),1leq ileq n)都按照定义枚举其约数计算,时间肯定爆炸。所以可以枚举约数,再枚举这些约数可以对哪些值给出贡献,那么计算一次狄利克雷卷积的复杂度就是(O(nlogn)),总复杂度(O(nlognlogk))

/*
刚开始要将ans初始化为单位元,即ans[2~n]=0,ans[1]=1,这样最初乘一个函数还是这个函数本身,即1 
初始化x为I,I(n)=1 
注:1.两个函数狄利克雷卷积是个函数 
2.加两个数取模不能直接用-=mod 
*/
#include<cstdio>
#include<cctype>
#include<cstring>
#define gc() getchar()
typedef long long LL;
const int N=1e5+5,mod=1e9+7;

int n,k;
LL f[N],ans[N],tmp[N],x[N];

inline int read()
{
	int now=0,f=1;register char c=gc();
	for(;!isdigit(c);c=gc()) if(c=='-') f=-1;
	for(;isdigit(c);now=now*10+c-'0',c=gc());
	return now*f;
}
#define Mod(x) x>=mod?x-=mod:0
void Dirichlet(LL *a,LL *b)//a*b
{
	memset(tmp,0,sizeof tmp);
	for(int i=1;i*i<=n;++i)
	{
		tmp[i*i]+=a[i]*b[i]%mod, Mod(tmp[i*i]);
		for(int j=i+1;i*j<=n;++j)//下边加上a[i]*b[j]和a[j]*b[i],所以j从i+1开始即可 
			(tmp[i*j]+=a[i]*b[j]%mod+a[j]*b[i]%mod)%=mod;//注意这加两个数不能一步用Mod取模。。
	}
	memcpy(a,tmp,sizeof tmp);
}
void Solve()
{
	for(int i=1;i<=n;++i) x[i]=1,ans[i]=0;//x:I^0
	ans[1]=1;//ans:e
	for(;k;k>>=1,Dirichlet(x,x))
		if(k&1) Dirichlet(ans,x);
	Dirichlet(ans,f);
	for(int i=1;i<=n;++i) printf("%lld%c",ans[i],i==n?'
':' ');//空格及换行符有要求 
}

int main()
{
	for(int t=read();t--;)
	{
		n=read(),k=read();
		for(int i=1;i<=n;++i) f[i]=read();
		Solve();
	}
	return 0;
}
原文地址:https://www.cnblogs.com/SovietPower/p/8351671.html