zoukankan      html  css  js  c++  java
  • SP14175 题解

    UPD:增加了利用斯特林数性质的正宗推式子做法。

    本题首A。

    首先有个显然的式子:

    [egin{aligned}sum_{i}sum_{j|i}j^k&=sum_{j}sum_{j|i}j^k\&=sum_{j}j^kleft[frac{n}{j} ight]end{aligned} ]

    然后你兴高采烈地准备拿拉格朗日插值去水题的时候,你发现一个悲惨的事实:这屑模数不是质数……

    出题人插头地爬!

    于是我们考虑找到

    [sum_{i=1}^{n}i^k ]

    的通项公式。

    首先我们有个引理:

    (d_{(k)}f(x)) 表示 (f(x))(k) 阶差分,有:

    [sum_{i=1}^{n}f(i)=sum_{i=1}^{n}inom{n}{i}d_{(i-1)}f(1) ]

    写个简单的证明:

    显然有

    [d_{(i-1)}f(1)=sum_{k=1}^{i}(-1)^{i-k}f(k)inom{i-1}{k-1} ]

    所以

    [egin{aligned}sum_{i=1}^{n}inom{n}{i}d_{(i-1)}f(1)&=sum_{i=1}^{n}inom{n}{i}sum_{k=1}^{i}(-1)^{i-k}f(k)inom{i-1}{k-1}\&=sum_{k=1}^{n}f(k)sum_{i=k}^{n}inom{i-1}{k-1}(-1)^{i-k}\&=sum_{k=1}^{n}f(k)end{aligned} ]

    得证!

    所以我们就可以计算啦!

    [F(k)=sum_{i=1}^{n}i^k ]

    考虑把每个(F(k))都写成形如

    [sum_{i=1}^{k}inom{n}{i}a_{k,i} ]

    的形式。

    先讲一下怎么计算(inom{n}{i})

    考虑(O(kpi(k)))预处理出(i!)的质因子,然后暴力除,复杂度(O(klog k))

    下面考虑怎么计算(a_{n,k})

    第一种方法:由定义可知:

    [a_{n,k}=d_{(i-1)}f(1)=sum_{k=1}^{i}(-1)^{i-k}f(k)inom{i-1}{k-1} ]

    暴力计算。复杂度(O(k^3+Tk^2log ksqrt n))

    当然如果你不把这玩意展开你可以做到预处理(O(k^2))

    第二种方法:考虑手算出几组小数据,然后试图扔到oeis里面。

    [F(0)=inom{n}{1} ]

    [F(1)=inom{n}{1}+inom{n}{2} ]

    [F(2)=inom{n}{1}+3inom{n}{2}+2inom{n}{3} ]

    [F(3)=inom{n}{1}+7inom{n}{2}+12inom{n}{3}+6inom{n}{4} ]

    取出系数1 1 1 1 3 2 1 7 12 6,扔进 oeis,发现它就是 A028246

    然而他提供的东西太少了……所以我们还是要写个小程序(雾

    翻到 formula 板块,我们惊奇的发现:

    [a_{n,k}=egin{Bmatrix}n\k end{Bmatrix}(k-1)! ]

    于是我们考虑在(O(k^2))的时间内预处理出(a_{n,k}),然后就可以做了。

    总复杂度:(O(k^2+Tk^2log ksqrt n))

    上面那个式子的推导过程太复杂了……我们需要一种新的推式子的做法。

    我们关注一下第二类斯特林数的定义,可以得到一个显然但是很重要的式子:

    [m^n=sum_{k=0}^{m}inom{n}{k}egin{Bmatrix}m\k end{Bmatrix}k! ]

    (m) 求和:

    [sum_{i=1}^{n}i^k=sum_{i=1}^{n}sum_{j=0}^{i}inom{k}{j}egin{Bmatrix}i\j end{Bmatrix}j! ]

    代码:

    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    const ll N=11;
    ll T,n,m,k;
    ll a[N+5][N+5],fac[N+5],pr[5]={2,3,5,7,11},mi[N+5][5];
    char ss[1<<17],*A=ss,*B=ss;
    inline char gc(){return A==B&&(B=(A=ss)+fread(ss,1,1<<17,stdin),A==B)?-1:*A++;}
    inline void read(ll&x){
        x=0;char s=gc();
        while(!isdigit(s))s=gc();
        while(isdigit(s))x=x*10+(s^'0'),s=gc();
    }
    void print(ll x){
        if(x<0)putchar('-'),x=-x;
        if(x>9)print(x/10);
        putchar(x%10+'0');
    }
    namespace Solution{
    	void init(){
    		a[1][1]=1;
    		for(ll i=2;i<=N;i++)
    			for(ll j=1;j<=i;j++)
    				a[i][j]=j*a[i-1][j]+a[i-1][j-1];
    		fac[0]=1;
    		for(ll i=1;i<=N;i++)fac[i]=fac[i-1]*i;
    		for(ll i=1;i<=N;i++){
    			ll t=i;
    			for(ll j=0;j<5;j++){
    				mi[i][j]=mi[i-1][j];
    				while(t%pr[j]==0)mi[i][j]++,t/=pr[j];
    			}
    		}
    		for(ll i=1;i<=N;i++)
    			for(ll j=1;j<=i;j++)
    				a[i][j]*=fac[j-1];
    	}
    	ll t[10],ans;
    	ll C(ll n,ll k){
    		ll ret=1;
    		for(ll i=0;i<k;i++)t[i]=n-i;
    		for(ll i=0;i<5;i++){
    			ll cnt=mi[k][i];
    			for(ll j=n-n/pr[i]*pr[i];j<k&&cnt;j+=pr[i]){
    				while(cnt&&t[j]%pr[i]==0)t[j]/=pr[i],cnt--;
    			}
    		}
    		for(ll i=0;i<k;i++)ret=(__int128)ret*t[i]%m;
    		return ret;
    	}
    	ll calc(ll n,ll k){
    		ll ret=0;
    		for(ll i=1;i<=k+1;i++)ret=(ret+(__int128)C(n,i)*a[k+1][i]%m)%m;
    		return ret;
    	}
    	void solve(){
    		read(n),read(k),read(m);ans=0;
    		for(ll l=1,r;l<=n;l=r+1){
    			r=n/(n/l);
    			ans=(ans+(__int128)n/l*((calc(r,k)-calc(l-1,k)+m)%m)%m)%m;
    		}
    		print(ans),puts("");
    	}
    } 
    int main(){
    	Solution::init();
    	for(read(T);T--;){
    		Solution::solve();
    	}
    	return 0;
    }
    
  • 相关阅读:
    【LeetCode】336.回文对(前缀树,散列,暴力三种方法,java实现)
    Python爬虫入门教程 72-100 分布式爬虫初步解析-配好环境肝完一半
    net.ipv4.tcp_max_tw_buckets=10
    tcp_fin_timeout
    程序停掉close_wait立马回收
    net.ipv4.tcp_tw_reuse = 1 重用time_wait
    做了3年数据报表却毫无进步?看过这3种方法的人,都被领导重视了
    【LeetCode】212.单词搜索 II (前缀树两种方法实现)
    Linux C/C++编程之(十三)系统IO函数
    《java入门第一季》之面向对象匿名内部类面试题
  • 原文地址:https://www.cnblogs.com/happydef/p/13492629.html
Copyright © 2011-2022 走看看