类欧几里得算法小结

版权声明:本文是蒟蒻写出来的,神犇转载也要说一声哦! https://blog.csdn.net/WerKeyTom_FTD/article/details/53812718

基本定义

f(a,b,c,n)=ni=0ai+bc
g(a,b,c,n)=ni=0iai+bc
h(a,b,c,n)=ni=0ai+bc2
m=an+bc
为了方便,接下来下取整都不写下取整,出现除法就是整除

推f

当a>=c或b>=c的时候, 容易得到
f(a,b,c,n)=(a/c)n(n+1)/2+(b/c)(n+1)+f(a%c,b%c,c,n)
如果a和b均小于c呢?
f(a,b,c,n)=ni=0mj=1[(ai+b)/c>=j]
这条式子可以这样理解,通过几何意义。
f的式子的意义其实是一条倾斜直线与x=0和x=n两条铅垂线组成的直角梯形内整点数量(不包括纵坐标为0的点,包括在线上的点)。
接下来我们这条式子就是在枚举一条纵线,统计上面整点的个数,我们现在改为枚举一条横线,统计横线上整点的个数,就是这条式子了。
f(a,b,c,n)=ni=0m1j=0[(ai+b)/c>=j+1]
f(a,b,c,n)=ni=0m1j=0[ai>=cj+cb]
f(a,b,c,n)=ni=0m1j=0[ai>cj+cb1]
f(a,b,c,n)=ni=0m1j=0[i>(cj+cb1)/a]
f(a,b,c,n)=mj=0(n(cj+cb1)/a)
f(a,b,c,n)=nmf(c,cb1,a,m1)
继续递归处理。
可以发现只看第一和第三维表示状态的话,每次由(a,c)变成(c,a%c),因此复杂度与欧几里得算法一致。

推g

a>=c或b>=c的时候,容易得到
g(a,b,c,n)=(a/c)n(n+1)(2n+1)/6+(b/c)n(n+1)/2+g(a%c,b%c,c,n)
注意这里0到n的二次方和公式就是n(n+1)(2n+1)/6
接下来我们来推a和b均小于c的情况。
首先先和f一样把下取整那玩意变换一下
g(a,b,c,n)=ni=0imj=1[(ai+b)/c>=j]
g(a,b,c,n)=ni=0im1j=0[i>(cj+cb1)/a]
g(a,b,c,n)=1/2m1j=0(n+1+(cj+cb1)/a)(n(cj+cb1)/a)
这里看的有点难受,f的时候转化的其实是0次方和,而g是1次方和,所以这里用了个求和公式。
g(a,b,c,n)=1/2m1j=0n(n+1)(cj+cb1)/a[(cj+cb1)/a]2
拆出来了
g(a,b,c,n)=1/2[n(n+1)mf(c,cb1,a,m1)h(c,cb1,a,m1)]
看来g比较麻烦,我们还要会推h。

推h

a>=c或b>=c的时候,h也是很麻烦的……
三项式的平方,拆出来大概是这样
h(a,b,c,n)=(a/c)2n(n+1)(2n+1)/6+(b/c)2(n+1)+(a/c)(b/c)n(n+1)+h(a%c,b%c,c,n)+2(a/c)g(a%c,b%c,c,n)+2(b/c)f(a%c,b%c,c,n)
接下来a和b均小于c的情况,我们要转化一下思路了。
如何获得一个n^2?
n2=2n(n+1)2n=2ni=0in
有了思路我们来推h
h(a,b,c,n)=ni=0(2(ai+b)/cj=1j(ai+b)/c)
可以想到交换主体。
h(a,b,c,n)=2m1j=0(j+1)ni=0[(ai+b)/c>=j+1]f(a,b,c,n)
h(a,b,c,n)=2m1j=0(j+1)ni=0[i>(cj+cb1)/a]f(a,b,c,n)
h(a,b,c,n)=2m1j=0(j+1)(n(cj+cb1)/a)f(a,b,c,n)
h(a,b,c,n)=nm(m+1)2g(c,cb1,a,m1)2f(c,cb1,a,m1)f(a,b,c,n)

实现

我的实现是直接返回f、g和h分别算出来的值,这样是比较简单的。
边界是a为0,返回三个0。
upd:
好像边界我写错了,大家随手推推吧。
写得挺丑,加了的强制转换其实没用……

dong likegcd(ll a,ll b,ll c,ll n){
    if (!a){
        gjx.f=gjx.g=gjx.h=0;
        return gjx;
    }
    if (a>=c||b>=c){
        zlt=likegcd(a%c,b%c,c,n);
        gjx.f=((ll)(a/c)*n%mo*(n+1)%mo*ni2%mo+(ll)(b/c)*(n+1)%mo)%mo;
        (gjx.f+=zlt.f)%=mo;
        gjx.g=((ll)(a/c)*n%mo*(n+1)%mo*(2*n+1)%mo*ni6%mo+(ll)(b/c)*(n+1)%mo*n%mo*ni2%mo)%mo;
        (gjx.g+=zlt.g)%=mo;
        gjx.h=(ll)(a/c)*(a/c)%mo*n%mo*(n+1)%mo*(2*n+1)%mo*ni6%mo;
        (gjx.h+=(ll)(b/c)*(b/c)%mo*(n+1)%mo)%=mo;
        (gjx.h+=(ll)(a/c)*(b/c)%mo*n%mo*(n+1)%mo)%=mo;
        (gjx.h+=(ll)2*(a/c)%mo*zlt.g%mo)%=mo;
        (gjx.h+=(ll)2*(b/c)%mo*zlt.f%mo)%=mo;
        (gjx.h+=zlt.h)%=mo;
        return gjx;
    }
    ll m=((ll)a*n+(ll)b)/(ll)c;
    zlt=likegcd(c,c-b-1,a,m-1);
    gjx.f=((ll)n*m%mo-(ll)zlt.f)%mo;
    gjx.g=((ll)(n+1)*n%mo*m%mo-(ll)zlt.f-(ll)zlt.h)%mo;
    gjx.g=(ll)gjx.g*ni2%mo;
    gjx.h=((ll)n*m%mo*(m+1)%mo-(ll)2*zlt.g%mo-(ll)2*zlt.f%mo-(ll)gjx.f)%mo;
    return gjx;
}

没有更多推荐了,返回首页