GDKOI2016 Day2 T4 國小生數學題

T4 國小生數學題

求

∑i=1n1imodpk

首先，我們把原式分成兩部分，一部分對于p有逆元，另一部分沒有。

即原式=

1p∑i=1⌊np⌋1i+∑a=1p−1∑i=0⌊np⌋−11a+i∗pmodpk

剩餘一部分，由于個數少于p，暴力做就行了。

然後分類讨論

part 1

1p∑i=1⌊np⌋1imodpk

一個明顯的結論

若 amodb=c

則 akmodbk=ck

設 f(n,k)=∑i=1n1imodpk

那麼，原式= f(⌊np⌋,k+1)modpk+1p

遞歸下去做就好了。

part 2

∑a=1p−1∑i=0⌊np⌋−11a+i∗pmodpk

原式= ∑a=1p−1∑i=0⌊np⌋−1a−1∗11+i∗a−1∗pmodpk

連結一個泰勒展開，當 x∞=0 時

11−x=x0+x1+x2+...+x∞

是以原式可以寫成

∑a=1p−1∑i=0⌊np⌋−1a−1∗∑b=0k−1(−i∗a−1∗p)bmodpk

為什麼隻到k-1次幂呢，因為這是在模 pk 的意義下進行的，大于等于k的話就為0了。

考慮交換枚舉順序，原式=

∑a=1p−1a−1∗∑b=0k−1(−a−1∗p)b∗∑i=0⌊np⌋−1ibmodpk

這樣求一個 自然數幂和就可以解決了。

#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[][];
ll f[];
ll qsm(ll x,ll y){
    if (!y) return ;
    ll t=qsm(x,y / );
    t=t*t%m;
    if (y%2) t=t*(x%m)%m;
    return t;
}
ll get(ll n,ll k){
    if (n<=) return ;
    su[][]=;ll i,j;
    fo(i,,k) su[i][]=,su[i][i]=;
    fo(i,,k)
        fo(j,,i-)
            su[i][j]=(su[i-][j-]+(i-)*su[i-][j]%m)%m;
    fo(i,,k) f[i]=;f[]=(n+)%m;
    if (n%2==) f[]=n/*(n+)%m;
    else f[]=n*(n+)/%m;
    fo(i,,k){
        t=;
        fo(j,n+-i,n+)
            if (j%(i+)==) t=t*j/(i+)%m;else t=t*j%m;
        f[i]=t;
        fo(j,,i-){
            if ((j+i)%2==) l=;else l=-;
            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=;y=;
        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 ;
    ll ans=,t,mo,r,i;
    r=(n/p)*p;
    fo(a,,p-){
        t=getny(a,m);
        fo(b,,k-) ans=(ans+t*qsm(-t*p,b)%m*get(n/p-,b)%m)%m;
    }
    fo(i,r+,n) ans=(ans+getny(i,m))%m;
    mo=m;
    m=m*p;
    ans=(ans+solve(n/p,k+)/p)%mo;
    ans=(ans%mo+mo)%mo;
    return ans;
}
int main(){
    scanf("%lld%lld%lld",&p,&k,&n);
    m=;
    fo(i,,k) m=m*p;
    printf("%lld\n",solve(n,k));
}