[GDKOI2016]小学生数学题

题目大意

求

$$\sum_{i=1}^ni^{-1}$$
对 $p^k$取模
p<=100000,n*p^k<=10^18

分成两组

我们设$f(n,k)表示\sum_{i=1}^ni^{-1}对p^k取模的结果$
我们假设n对p取模为0，不是的话把多余部分暴力然后把n变为模p为0。
然后n以内的可以分成两类：是p的倍数和不是p的倍数。
设第一类的结果为S，那么

$$S=\sum_{i=1}^{\frac{n}{p}}(i*p)^{-1} \mod p^k$$
$$S=\frac{1}{p}\sum_{i=1}^{\frac{n}{p}}i^{-1} \mod p^k$$
有一个数论技巧，如果 $a \mod b=c$，那么 $ka \mod kb=kc$
具体证明很简单，这里不证了。
所以 $S=\frac{f(\frac{n}{p},k+1)}{p}$
递归处理即可！而且这样做n*p^k并没有变化。
设第二类的结果为T，显然我们可以把第二类数弄成以下形式
$a+i*p其中1<=a<=p-1，0<=i<\frac{n}{p}$
那么显然我们现在需要解决的问题是如何求形如a+i*p的逆元
首先显然有 $\frac{1}{1-x}=x^0+x^1+x^2+...$
注意这是模意义下，x含有一个因子p，用等比数列求和就得证了。
没必要做到无穷，后面在模意义下全是0
$\frac{1}{1-x}=x^0+x^1+x^2+...x^k$
好办喽， $(a+i*p)^{-1}=\frac{1}{a+i*p}=\frac{a^{-1}}{1+i*p*a^{-1}}$
于是 $T=\sum_{a=1}^{p-1}a^{-1}*\sum_{i=0}^{\frac{n}{p}-1}\sum_{b=0}^k(-i*p*a^{-1})^b$
交换一下
$T=\sum_{a=1}^{p-1}a^{-1}*\sum_{b=0}^k(-a^{-1}*p)^b\sum_{i=0}^{\frac{n}{p}-1}i^b$
枚举a,b，后面的自然数幂和可以用各种方法求，具体参见我blog。
我的code由于没有考虑爆long long问题有些数据会出错。

#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)%m)%m;
    return t;
}
ll get(ll n,ll k){
    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;
    if (!n) return 0;
    if (n<=k){
        t=0;
        fo(i,1,n){
            l=1;
            fo(j,1,k)
                l=l*i%m;
            t=(t+l)%m;
        }
        return t;
    }
    fo(i,0,k) f[i]=0;
    if (n%2==0) f[1]=n/2*(n+1)%m;
    else f[1]=n*(n+1)/2%m;
    f[0]=n%m;
    f[1]=(f[1]+1)%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]-1)%m+m)%m;
}
void gcd(ll a,ll b,ll &x,ll &y){
    if (!b){
        x=1;
        y=0;
        return;
    }
    else{
        gcd(b,a%b,x,y);
        swap(x,y);
        y-=x*(a/b);
    }
}
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-(n/p)*p;
    n-=r;
    fo(a,1,p-1){
        t=getny(a,m);
        fo(b,0,k)
            ans=(ans+t*qsm(-t*p,b)%m*get(n/p-1,b)%m)%m;
    }
    fo(i,n+1,n+r) 
        if (i%p!=0) ans=(ans+getny(i,m))%m;
    n+=r;
    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));
}