关闭

解决自然数幂和的各种方法

826人阅读 评论(4) 收藏 举报
分类:

前言

自然数幂和,

暴力

咳咳这个就不用说了。当然不要太暴力了,像求幂可以用快速幂。

高斯消元

我们仔细观察发现,对于求自然数k次方和的公式应该是k+1次的。那么我们可以想到,列个方程组让高斯消元解决,然后得到各项系数后直接代入即可。
复杂度O(k3)

矩阵乘法

发现方法有点问题,先留坑

倍增

我们设f(n,k)=ni=1ik
那么现在我们要计算f(n,k)怎么办呢?
我们采用分治思想。
如果n是奇数那么
f(n,k)=f(n1,k)+nk
否则,我们可以先求出n/2的f,然后对于n/2+1~n中的每个数都可以表示为n/2+i
那么
f(n,k)=f(n/2,k)+n/2i=1kj=0Cjkij(n/2)kj
调换一下
f(n,k)=f(n/2,k)+kj=0Cjkf(n/2,j)(n/2)kj
于是对于每一个n的f我们需要O(k2),一共有log n层,算上快速幂的复杂度
总共是O(k2lognlogk)

插值法

分拉格朗日插值法和牛顿插值法,我不会……

伯努利数

同不会

第一类斯特林数

该方法来自GEOTCBRL小学生数学题题解,%%%。
我们设

Sk(n)=i=0nik

那么这个东东要怎么搞呢?
先引入第一类斯特林数的定义:
我们记f(n,k)=n(n1)(n2)...(nk+1)
f(x,n)=k=0nss(n,k)xk

那么我们很容易推出其递推式
ss(n,m)=ss(n1,m1)(n1)ss(n1,m)
提一个有趣的一点:如果令su(n,m)=|ss(n,m)|
那么su(n,m)=su(n1,m1)+(n1)su(n1,m)
而且有一个性质su(n,m)=(1)n+mss(n,m)
而且su还有另外的含义,su(n,m)表示将 n 个不同元素构成m个圆排列的数目。
圆排列是什么?可以理解为两个排列经过旋转可以相同那么视为等价。
由于标程使用了su,接下来全程都介绍如何用su求自然数幂和。
基本性质:
su(n,n)=1
su(n,0)=0
有一个结论:
k!Ckj=f(j,k)
怎么证明呢?
左边是k!j!k!(jk)!=j!(jk)!=(jk+1)...j=f(j,k)
那这个结论就显然了。
上面的东西有啥用?等等就会用到。
我们知道ss(k,k)=1
那么jk=ss(k,k)jk
于是呢
jk=k!Ckjf(j,k)+ss(k,k)jk

jk=k!Ckj(f(j,k)ss(k,k)jk)

观察括号内的东西
f(j,k)ss(k,k)jk=i=0kss(k,i)jiss(k,k)jk

咦!
f(j,k)ss(k,k)jk=i=0k1ss(k,i)ji

那么
jk=k!Ckji=0k1ss(k,i)ji

于是
Sk(n)=j=0n(k!Ckji=0k1(1)i+ksu(k,i)ji)

Sk(n)=k!j=0nCkji=0k1(1)i+ksu(k,i)j=0nji

运用上面的结论
Sk(n)=k!Ck+1n+1i=0k1(1)i+ksu(k,i)j=0nji

Sk(n)=f(n+1,k+1)k+1i=0k1(1)i+ksu(k,i)j=0nji

我们预处理su需要O(k2),然后可以用O(k2)计算Sk。对于减号前面的东西我们直接O(k)做就好了。
如果要求模,该怎么办?
相对于某些方法需要对模数分解质因数,对每一个p^k做了后用中国剩余定理,该方法可以直接对减号前部分暴力扫一遍相乘并模。分母是k+1,分子是k+1个连续自然数相乘,那么一定有一个自然数是k+1的倍数,因此消除分母的影响,让该方法可以轻松解决任意模数的情况。而某些方法要求模数是质数,还有些方法要求模数拆分后最大的质数幂不能太大。而该方法则比较优秀。代码实现方面也比较简单,时间复杂度为O(k2),因此个人认为十分好用。
下面是该题做WYF的盒子的代码

#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 su[2000+10][2000+10],f[2000+10];
ll i,j,k,l,t,n,m,p;
ll qsc(ll x,ll y){
    ll a1=x/1000000,a2=x%1000000,b1=y/1000000,b2=y%1000000;
    ll t=a1*b1%p*1000000%p*1000000%p;
    t=(t+a2*b2%p);
    t=(t+a1*b2%p*1000000%p);
    t=(t+a2*b1%p*1000000%p);
    return t%p;
}
ll qsm(ll x,ll y){
    if (!y) return 1;
    ll t=qsm(x,y/2);
    t=qsc(t,t);
    if (y%2) t=qsc(t,x);
    return t;
}
ll get(ll n){
    if (!n) return 0;
    if (n<=k){
        t=0;
        fo(i,1,n){
            l=1;
            fo(j,1,k)
                l=qsc(l,i);
            t=(t+l)%p;
        }
        return t;
    }
    fo(i,0,k) f[i]=0;
    if (n%2==0) f[1]=qsc(n/2,n+1);
    else f[1]=qsc(n,(n+1)/2);
    f[1]=(f[1]+1)%p;
    f[0]=n%p;
    fo(i,2,k){
        t=1;
        fo(j,n+1-i,n+1)
            if (j%(i+1)==0) t=qsc(t,j/(i+1));else t=qsc(t,j);
        f[i]=t;
        fo(j,0,i-1){
            if ((j+i)%2==0) l=1;else l=-1;
            f[i]=((f[i]-qsc(l*su[i][j],f[j]))%p+p);
        }
        f[i]=f[i]%p;
    }
    return ((f[k]-1)%p+p)%p;
}
int main(){
    scanf("%lld%lld%lld%lld",&k,&n,&m,&p);
    if (n-m<=5000){
        t=0;
        fo(i,m,n) t=(t+qsm(i,k));
        printf("%lld\n",t%p);
        return 0;
    }
    su[0][0]=1;
    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]+qsc(i-1,su[i-1][j]));
    printf("%lld\n",((get(n)-get(m-1))%p+p)%p);
} 
1
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:251938次
    • 积分:10104
    • 等级:
    • 排名:第1813名
    • 原创:743篇
    • 转载:4篇
    • 译文:0篇
    • 评论:189条
    最新评论
    文章分类