【JZOJ5395】Count

Description
令f(i)表示与i互质的数的平均数*2
求∑ri=lf(i)k

Input
三个数l,r,k(1≤l≤r≤108,1≤k≤106)

Output
一个数答案，对998244353取模

Sample Input
2 3 1

Sample Output
5

Hint
1与2互质，f(2)=11∗2=2
1,2与3互质，f(3)=1+22∗2=3

---

Solution

易知f(1)=2，f(2)=2
当n>2时，f(n)=n，证:

与n互质的数有φ(n)个，又因为若a<n,gcd(a,n)=1则gcd(n−a,n)=1，且n>2时2|φ(n)，所以互质数和为φ(n)∗n2，f(i)=φ(n)n2φ(n)∗2=n

设S(n)=∑ni=1ik，则Ans=S(n)−1+2k

接下来问题就变成了求自然数幂和。

---

差分法

(a+1)k−ak=∑i=0k−1Cikai

∴(n+1)k−1=∑i=1n(i+1)k−ik=∑i=1n∑j=0k−1Cjkij=∑i=0k−1Cik∑j=1nji=∑i=0k−1CikS(i)

∴(n+1)k+1−1=∑i=0kCik+1S(i)

（把i=k之外的移项,Ckk+1=k+1）

(k+1)S(k)=(n+1)k+1−1−∑i=0k−1Cik+1S(i)

∴S(k)=(n+1)k+1−1−∑k−1i=0Cik+1S(i)k+1

但这种方法放到原问题上会超时

---

拉格朗日(Lagrange)插值法

引入

拉格朗日插值法的目的就是让你求一个多项式(或者是化简的通式)：
已知k+1个x以及对应多项式的值，求该k次k项式的通式

定义自变量值为x的多项式取值为L(x)
对于我们已知的xi对应多项式的值为yi
有L(x)=∑k+1i=0yili(x)
我们要找出l_i(x)的式子使得：
1. 当x=xd(某个已知的x),li(x)=[x==xi]
2. 其余情况多项式成立
那么，我们可以找到一种：li(x)=∏k+1i=1,i≠jx−xjxi−xj

当x=xi时，对于各个xj有x−xj=xi−xj即x−xjxi−xj=1∴li(x)=1
否则因为x一定是x1..k+1中的一个，那么∃x=xj(i≠j)，此时li(x)=0
既然如此，此式子对于xi都可得出yi，又因k+1个点可确定一个k次多项式，所以，其他的(x′,y′)依然满足

于是我们得出一种通式

L(n)=∑i=1k+1yi∏j=1,i≠jk+1n−xjxi−xj

回到原问题

我们知道自然数k次幂和的求和公式是个k+1次多项式

由差分表法的，S(k)=(n+1)k+1−1−∑k−1i=0Cik+1S(i)k+1，此时n是个自变量，(x+1)k+1二项式展开后有k+1次项，而S(i)(i<k)也依赖于前面的S(i)，数学归纳，可证S(i)不会超过k次项，故自然数幂和求和公式是一个k+1次多项式

那么我们要先制定一些x1...k+2求出y1...k+2，为了方便，指定xi=i，用快速幂求yi
利用L(n)=∑k+2i=1yi∏k+2j=1,i≠jn−xjxi−xj求解。

L(n)=∑i=1k+2yi∏j=1,i≠jk+2n−xjxi−xj=∑i=1k+2yi∏k+2j=1,i≠j(n−j)∏k+2j=1,i≠j(i−j)

i−j
只能是1,2,3,4,5,6,7,8,9,...,−1,−1,−2,−3,...,−k−1+i,−k−2+i
这可以用阶乘处理[(−1)k+2−i(i−1)!(k+2−i)!]−1

n−j只有k+2种，预处理∏(n−j)就好了，利用前缀积pre[i−1]，后缀积suf[i+1]排除掉xj=xi的情况。pre[i−1]∗suf[i+1]

S(n)=L(n)=∑i=1k+2yipre[i−1]suf[i+1][(−1)k+2−i(i−1)!(k+2−i)!]−1

预处理O(klogk)插值O(k)

#include<cstring>#include<cstdio>#include<cctype>#define mo 998244353using namespace std;typedef long long ll;ll l,r,k,frac[1001001],nyfrac[1001001],ny[1001001],s[1001001],f[1001001],g[1001001];ll qpow(ll a,ll ind){    ll ans=1;    for(;ind;ind>>=1,a=a*a%mo)if(ind&1)ans=ans*a%mo;return ans;}ll ask(ll n){    if(!n)return 0;    ll ans=0;    if(n<=k+2)ans=s[n];else{        f[0]=g[k+3]=1;        for(ll i=1;i<=k+2;i++)f[i]=f[i-1]*(n-i)%mo;        for(ll i=k+2;i;i--)g[i]=g[i+1]*(n-i)%mo;        for(ll i=1;i<=k+2;i++){            ans=(ans+nyfrac[i-1]*nyfrac[k+2-i]*((k+2-i&1)?-1:1)%mo*s[i]%mo*f[i-1]%mo*g[i+1])%mo;        }    }    ans=(ans-1+qpow(2,k)+mo)%mo;return ans;}int main(){    freopen("count.in","r",stdin);    freopen("count.out","w",stdout);    frac[0]=nyfrac[0]=1;    for(ll i=1;i<=1000010;i++)ny[i]=qpow(i,mo-2);    for(ll i=1;i<=1000010;i++)frac[i]=frac[i-1]*i%mo,nyfrac[i]=nyfrac[i-1]*ny[i]%mo;    scanf("%lld %lld %lld",&l,&r,&k);    for(ll i=1;i<=k+2;i++)s[i]=(s[i-1]+qpow(i,k))%mo;    printf("%lld",(ask(r)-ask(l-1)+mo)%mo);    fclose(stdin);fclose(stdout);    return 0;}