GDKOI2016 Day2 T4 小学生数学题

T4 小学生数学题

求

$$\sum_{i=1}^{n}{1\over i}\mod p^k$$

首先，我们把原式分成两部分，一部分对于p有逆元，另一部分没有。
即原式=

$${1\over p}\sum_{i=1}^{\lfloor {n\over p}\rfloor}{1\over i}+\sum_{a=1}^{p-1}\sum_{i=0}^{\lfloor {n\over p}\rfloor -1}{1\over {a+i*p}}\mod p^k$$
剩余一部分，由于个数少于p，暴力做就行了。
然后分类讨论

part 1

$${1\over p}\sum_{i=1}^{\lfloor {n\over p}\rfloor}{1\over i}\mod p^k$$
一个明显的结论
若 $$a\mod b=c$$
则 $$ak\mod bk=ck$$
设 $$f(n,k)=\sum_{i=1}^{n}{1\over i}\mod p^k$$
那么，原式= $$f({\lfloor {n\over p}\rfloor},k+1)\mod p^{k+1}\over p$$
递归下去做就好了。

part 2

$$\sum_{a=1}^{p-1}\sum_{i=0}^{\lfloor {n\over p}\rfloor -1}{1\over {a+i*p}}\mod p^k$$
原式= $$\sum_{a=1}^{p-1}\sum_{i=0}^{\lfloor {n\over p}\rfloor -1}a^{-1}*{1\over1+i*a^{-1}*p}\mod p^k$$
链接一个泰勒展开，当 $x^∞=0$时
$${1\over {1-x}}=x^0+x^1+x^2+...+x^∞$$
所以原式可以写成
$$\sum_{a=1}^{p-1}\sum_{i=0}^{\lfloor {n\over p}\rfloor -1}a^{-1}*\sum_{b=0}^{k-1}({-i*a^{-1}*p})^b\mod p^k$$
为什么只到k-1次幂呢，因为这是在模 $p^k$的意义下进行的，大于等于k的话就为0了。
考虑交换枚举顺序，原式=
$$\sum_{a=1}^{p-1}{a^{-1}}*\sum_{b=0}^{k-1}(-a^{-1}*p)^b*\sum_{i=0}^{\lfloor {n\over p}\rfloor -1}i^b\mod p^k$$
这样求一个 自然数幂和就可以解决了。

#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
ll i,j,k,l,t,n,m,a,b,p;
ll su[150][150];
ll f[150];
ll qsm(ll x,ll y){
    if (!y) return 1;
    ll t=qsm(x,y / 2);
    t=t*t%m;
    if (y%2) t=t*(x%m)%m;
    return t;
}
ll get(ll n,ll k){
    if (n<=0) return 0;
    su[0][0]=1;ll i,j;
    fo(i,1,k) su[i][0]=0,su[i][i]=1;
    fo(i,1,k)
        fo(j,1,i-1)
            su[i][j]=(su[i-1][j-1]+(i-1)*su[i-1][j]%m)%m;
    fo(i,0,k) f[i]=0;f[0]=(n+1)%m;
    if (n%2==0) f[1]=n/2*(n+1)%m;
    else f[1]=n*(n+1)/2%m;
    fo(i,2,k){
        t=1;
        fo(j,n+1-i,n+1)
            if (j%(i+1)==0) t=t*j/(i+1)%m;else t=t*j%m;
        f[i]=t;
        fo(j,0,i-1){
            if ((j+i)%2==0) l=1;else l=-1;
            f[i]=((f[i]-l*su[i][j]%m*f[j]%m)%m+m)%m;
        }
        f[i]=f[i]%m;
    }
    return f[k];
}
void gcd(ll a,ll b,ll &x,ll &y){
    if (!b){
        x=1;y=0;
        return;
    }
    ll xx,yy;
    gcd(b,a%b,xx,yy);
    x=yy;y=xx-a/b*yy;   
}
ll getny(ll a,ll b){
    ll x,y;
    gcd(a,b,x,y);
    x=(x%b+b)%b;
    return x;
}
ll solve(ll n,ll k){
    if (!n) return 0;
    ll ans=0,t,mo,r,i;
    r=(n/p)*p;
    fo(a,1,p-1){
        t=getny(a,m);
        fo(b,0,k-1) ans=(ans+t*qsm(-t*p,b)%m*get(n/p-1,b)%m)%m;
    }
    fo(i,r+1,n) ans=(ans+getny(i,m))%m;
    mo=m;
    m=m*p;
    ans=(ans+solve(n/p,k+1)/p)%mo;
    ans=(ans%mo+mo)%mo;
    return ans;
}
int main(){
    scanf("%lld%lld%lld",&p,&k,&n);
    m=1;
    fo(i,1,k) m=m*p;
    printf("%lld\n",solve(n,k));
}