[GDKOI2016]小学生数学题

记

$F(n)=sumlimits_{i=1}^{n}i^{-1}$

$G(n)=sumlimits_{i=1,i eq jp}^{n}i^{-1}$

我们要算$F(n)\%p^k$

那么

$F(n)\%p^k=frac{F( left lfloor frac{n}{p}  ight  floor )}{p}\%p^k+G(n)\%p^k$

我们知道$frac{F( left lfloor frac{n}{p}  ight  floor )}{p}\%p^k=frac{F( left lfloor frac{n}{p}  ight  floor )\%p^{k+1}}{p}$，其中$F( left lfloor frac{n}{p}  ight  floor )\%p^{k+1}$可以可以递归算，所以我们重点要考虑的是$G(n)\%p^k$

不妨设$p|n$，那么

$G(n)\%p^k=sumlimits_{a=1}^{p-1}sumlimits_{b=0}^{left lfloor frac{n}{p}  ight  floor -1}(a+bp)^{-1}\%p^k$

用广义二项式定理展开：

$(a+bp)^{-1}=sumlimits_{i=0}^{oo}C_{-1}^{i}a^{-1-i}b^ip^i=sumlimits_{i=0}^{oo}(-1)^ia^{-1-i}b^ip^i$

又因为是在模$p^k$意义下的，所以

$(a+bp)^{-1}\%p^k=sumlimits_{i=0}^{k-1}(-1)^ia^{-1-i}b^ip^i\%p^k$

所以

$G(n)\%p^k=sumlimits_{a=1}^{p-1}sumlimits_{b=0}^{left lfloor frac{n}{p}  ight  floor -1}sumlimits_{i=0}^{k-1}(-1)^ia^{-1-i}b^ip^i\%p^k$

$=sumlimits_{i=0}^{k-1}(-1)^ip^isumlimits_{a=1}^{p-1}a^{-1-i}sumlimits_{b=0}^{left lfloor frac{n}{p}  ight  floor -1}b^i\%p^k$

需要注意的一点是，用二项式定理的时候，规定$0^0=1$

我们枚举$i$和$a$，要$O(kp)$的时间复杂度；剩下的$sumlimits_{b=0}^{left lfloor frac{n}{p}  ight  floor -1}b^i\%p^k$就是自然数幂和，记$S_{i}(n)=sumlimits_{b=0}^{n}b^i\%p^k$，可以用$O(k^2)$的时间内预处理出所有的$S_{i}(left lfloor frac{n}{p}  ight  floor-1)\%p^k(0leq ileq k-1)$

如果$n$不是$p$的倍数，剩下的零头乱搞就行了

所以每递归一次的时间复杂度是$O(kp+k^2)$

所以总的时间复杂度是$O(log_{p}n(kp+k^2))$

#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<fstream>
#include<algorithm>
#include<cstring>
#include<string>
#include<cmath>
#include<queue>
#include<stack>
#include<map>
#include<utility>
#include<set>
#include<bitset>
#include<vector>
#include<functional>
#include<deque>
#include<cctype>
#include<climits>
#include<complex>
#include<cassert> 
//#include<bits/stdc++.h>适用于CF,UOJ，但不适用于poj
 
using namespace std;

typedef long long LL;
typedef double DB;
typedef pair<int,int> PII;
typedef pair<DB,DB> PDD;
typedef complex<DB> CP;
typedef vector<int> VI;

#define mmst(a,v) memset(a,v,sizeof(a))
#define mmcy(a,b) memcpy(a,b,sizeof(a))
#define fill(a,l,r,v) fill(a+l,a+r+1,v)
#define re(i,a,b)  for(i=(a);i<=(b);i++)
#define red(i,a,b) for(i=(a);i>=(b);i--)
#define fi first
#define se second
#define mp(a,b) make_pair(a,b)
#define pb(a) push_back(a)
#define SF scanf
#define PF printf
#define two(k) (1<<(k))
#define SZ(x) (int(x.size()))
#define all(x) (x).begin(),(x).end()
#define ire(i,v,x) for(i=0,v=i<SZ(x)?x[i]:0;i<SZ(x);v=x[++i])


template<class T>inline T sqr(T x){return x*x;}
template<class T>inline void upmin(T &t,T tmp){if(t>tmp)t=tmp;}
template<class T>inline void upmax(T &t,T tmp){if(t<tmp)t=tmp;}

inline int sgn(DB x){if(abs(x)<1e-9)return 0;return(x>0)?1:-1;}
const DB Pi=acos(-1.0);

int gint()
  {
        int res=0;bool neg=0;char z;
        for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar());
        if(z==EOF)return 0;
        if(z=='-'){neg=1;z=getchar();}
        for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar());
        return (neg)?-res:res; 
    }
LL gll()
  {
      LL res=0;bool neg=0;char z;
        for(z=getchar();z!=EOF && z!='-' && !isdigit(z);z=getchar());
        if(z==EOF)return 0;
        if(z=='-'){neg=1;z=getchar();}
        for(;z!=EOF && isdigit(z);res=res*10+z-'0',z=getchar());
        return (neg)?-res:res; 
    }

const LL maxp=100100;
const LL maxk=1010;

LL multiply(LL a,LL b,LL MOD)
  {
      LL t=(LL)((DB)a*b/MOD);
      return a*b-t*MOD;
  }

LL su[maxk][maxk],sb[maxk];
void calc(LL n,LL d,LL m)
  {
      LL i,j;
      su[0][0]=1;
      re(i,1,d)
        {
            su[i][0]=0;
            re(j,1,d-1)su[i][j]=(multiply(i-1,su[i-1][j],m)+su[i-1][j-1])%m;
            su[i][i]=1;
        }
      sb[0]=n%m;
      re(i,1,d)
        {
            sb[i]=1;
            red(j,n+1,n+1-i)if(j%(i+1)==0)sb[i]=multiply(sb[i],j/(i+1),m); else sb[i]=multiply(sb[i],j,m);
            LL l=(i%2==0)?1:-1;
            re(j,0,i-1){(sb[i]-=l*multiply(su[i][j],sb[j],m))%=m;l=-l;}
        }
  }

LL reva[maxp],powreva[maxp];

LL F(LL n,LL p,LL k,LL m)
  {
      LL i,j,res=0;
      if(n<p)
        {
            reva[1]=1;re(j,2,n)reva[j]=multiply(reva[m%j],(m-m/j),m);
            re(j,1,n)(res+=reva[j])%=m;
            return res;
        }
      calc(n/p-1,k-1,m);
        (sb[0]+=1)%=m;
      reva[1]=1;re(j,2,p-1)reva[j]=multiply(reva[m%j],(m-m/j),m);
      re(j,1,p-1)powreva[j]=1;
      LL sp=1,sa;
      re(i,0,k-1)
        {
            sa=0;
            re(j,1,p-1)powreva[j]=multiply(powreva[j],reva[j],m),(sa+=powreva[j])%=m;
            (res+=multiply(multiply(sp,sa,m),sb[i],m))%=m;
            sp=multiply(sp,-p,m);
        }
      if(n%p!=0)
        {
            re(j,1,n%p)powreva[j]=1;
            LL sp=1,sa;
            re(i,0,k-1)
              {
                  sa=0;
                  re(j,1,n%p)powreva[j]=multiply(powreva[j],reva[j],m),(sa+=powreva[j])%=m;
                  (res+=multiply(sp,sa,m))%=m;
                  sp=multiply(sp,(-p)*(n/p),m);
              }
          }
        (res+=F(n/p,p,k+1,m*p)/p)%=m;
        return res;
  }

int main()
  {
      freopen("math.in","r",stdin);
      freopen("math.out","w",stdout);
      LL i,p=gll(),k=gll(),n=gll(),m=1;re(i,1,k)m*=p;
      cout<<(F(n,p,k,m)%m+m)%m<<endl;
      return 0;
  }