类欧几里得算法学习笔记

Larunatrecy

2022-07-06 21:17:07

Algo. & Theory

引入

类欧几里得算法可以在O(\log n)形如如下问题的：

\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor

设 f(n,a,b,c)=\sum\limits_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor

首先，先把 a\geq c 或 b\geq c 的情况通过对 c 取模化为a<c,b<c 的情况

f(n,a,b,c)=\sum_{i=0}^n\lfloor\frac{(c\lfloor\frac{a}{c}\rfloor+(a\bmod c))i+c\lfloor\frac{b}{c}\rfloor+(b\bmod c)}{c}\rfloor =\sum_{i=0}^n\lfloor\frac{a}{c}\rfloor i+\lfloor\frac{b}{c}\rfloor+\lfloor\frac{(a\bmod c)i+(b\bmod c)}{c}\rfloor =\lfloor\frac{a}{c}\rfloor \frac{n(n+1)}{2}+(n+1)\lfloor\frac{b}{c}\rfloor+f(n,a\bmod c,b\bmod c,c)

设 m=\lfloor\frac{an+b}{c}\rfloor

接下来考虑求后半部分的值：

f(n,a,b,c)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor =\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}1

交换求和号：

=\sum_{j=0}^{m-1}\sum_{i=0}^n[j<\lfloor\frac{ai+b}{c}\rfloor]

接下来考虑拆掉里边的取整：

=\sum_{j=0}^{m-1}\sum_{i=0}^n[j+1\leq \frac{ai+b}{c}] =\sum_{j=0}^{m-1}\sum_{i=0}^n[cj+c-b-1< ai] =\sum_{j=0}^{m-1}\sum_{i=0}^n[\lfloor \frac{cj+c-b-1}{a}\rfloor < i]

这样做的好处是里边的i独立了，符合的i的个数可以直接算出来

=\sum_{j=0}^{m-1}n-\lfloor \frac{cj+c-b-1}{a}\rfloor

后边是一个递归的形式

=nm-f(m-1,c,c-b-1,a)

可以发现，我们把a和c的位置交换了，也就是一个类似于辗转相除的过程，复杂度和欧几里得算法的复杂度都是O(logn)，这也是类欧名字的由来。

代码大致长这个样子：

LL f(LL n,LL a,LL b,LL c)
{
    LL ac=a/c,bc=b/c;
    LL m=(a*n+b)/c;
    if(a==0)
    {
        return (n+1)*bc;
    }
    if(a>=c||b>=c)
    {
        return f(n,a%c,b%c,c)+n*(n+1)/2*ac+(n+1)*bc;
    }
    return n*m-f(m-1,c,c-b-1,a);
}

拓展

求

g(n,a,b,c)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor i h(n,a,b,c)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor^2

求g

首先类似的，通过对c取模使得a,b小于c

g(n,a,b,c)=g(n,a\bmod c,b\bmod c,c)+\lfloor \frac{a}{c}\rfloor\frac{n(n+1)(2n+1)}{6}+\lfloor \frac{b}{c}\rfloor \frac{n(n+1)}{2}

然后仿照 f 的求法

g(n,a,b,c)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor i =\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}i

设m=\lfloor\frac{an+b}{c}\rfloor

=\sum_{j=0}^{m-1}\sum_{i=0}^n[j<\lfloor\frac{ai+b}{c}\rfloor]i =\sum_{j=0}^{m-1}\sum_{i=0}^n[\lfloor \frac{cj+c-b-1}{a}\rfloor < i]i

设p=\lfloor \frac{cj+c-b-1}{a}\rfloor

=\sum_{j=0}^{m-1}\frac{(p+1+n)(n-p)}{2} =\frac{1}{2}\sum_{j=0}^{m-1}n+n^2-p^2-p =\frac{1}{2}(mn(n+1)-\sum_{j=0}^{m-1}p-\sum_{j=0}^{m-1}p^2) =\frac{1}{2}(mn(n+1)-f(m-1,c,c-b-1,a)-h(m-1,c,c-b-1,a))

求h

首先还是取模,可以得到:

h(n,a,b,c)=h(n,a\mod c,b\mod c,c) +2\lfloor \frac{b}{c}\rfloor f(n,a\mod c,b\mod c,c)+2\lfloor \frac{a}{c}\rfloor g(n,a\mod c,b\mod c,c) +\lfloor \frac{a}{c}\rfloor^2\frac{n(n+1)(2n+1)}{6}+\lfloor \frac{b}{c}\rfloor^2(n+1)+\lfloor \frac{a}{c}\rfloor \lfloor \frac{b}{c}\rfloor n(n+1)

接着求h(n,a,b,c)

设 m=\lfloor\frac{an+b}{c}\rfloor,p=\lfloor \frac{cj+c-b-1}{a}\rfloor

为了避免出现求和号相乘的情况，我们换一个形式

n^2=2\frac{n(n+1)}{2}-n=2\sum_{i=1}^ni-n

所以

h(n,a,b,c)=\sum_{i=0}^n(2\sum_{j=1}^{\lfloor\frac{ai+b}{c}\rfloor}j-\lfloor\frac{ai+b}{c}\rfloor) =2\sum_{i=0}^n\sum_{j=1}^{\lfloor\frac{ai+b}{c}\rfloor}j-f(n,a,b,c) =2\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}(j+1)-f(n,a,b,c)

只看前面

=\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[j<\lfloor\frac{ai+b}{c}\rfloor] =\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[\lfloor \frac{cj+c-b-1}{a}\rfloor < i] =\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[p < i] =\sum_{j=0}^{m-1}(j+1)(n-p) =\sum_{j=0}^{m-1}jn+n-jp-p =\frac{1}{2}nm(m+1)-\sum_{i=0}^{m-1}i p-\sum_{i=0}p =\frac{1}{2}nm(m+1)-g(m-1,c,c-b-1,a)-f(m-1,c,c-b-1,a)

总结

发现三个函数互相调用，但是函数参数的变化是相同的，因此可以把三个函数值存到一个结构体里一起递归，总时间复杂度O(logn)

struct node
{
    LL f,g,h;
    node()
    {
        f=g=h=0;
    }
};
node Euclid(LL n,LL a,LL b,LL c)
{
    LL ac=a/c,bc=b/c;
    LL m=(a*n+b)/c;
    LL n1=n+1;
    LL n2=n*2+1;
    node ans;
    if(!a)
    {
        ans.f=bc*n1%mod;
        ans.g=n*n1%mod*inv2%mod*bc%mod;
        ans.h=bc*bc%mod*n1%mod;
        return ans; 
    }
    if(a>=c||b>=c)
    {
        ans.f=n*n1%mod*ac%mod*inv2%mod+n1*bc%mod;ans.f%=mod;
        ans.g=n*n1%mod*n2%mod*inv6%mod*ac%mod+bc*n%mod*n1%mod*inv2%mod;ans.g%=mod;
        ans.h=ac*ac%mod*n%mod*n1%mod*n2%mod*inv6%mod+bc*bc%mod*n1%mod+ac*bc%mod*n%mod*n1%mod;ans.h%=mod;
        node x=Euclid(n,a%c,b%c,c);
        ans.h+=2*bc%mod*x.f%mod+2*ac*x.g%mod+x.h;ans.h%=mod;
        ans.g+=x.g;ans.f+=x.f;ans.g%=mod;ans.f%=mod;
        return ans;     
    }
    node x=Euclid(m-1,c,c-b-1,a);
    ans.f=(n*m%mod-x.f+mod)%mod;
    ans.g=(m*n%mod*n1%mod-x.h-x.f+mod)%mod;
    ans.g=(ans.g%mod+mod)%mod;
    ans.g=ans.g*inv2%mod;
    ans.h=m*n%mod*(m+1)%mod-2*x.g%mod-2*x.f%mod-ans.f;
    ans.h=(ans.h%mod+mod)%mod;
    return ans;
}

例题

P5171 Earthquake

因为

ax+by\leq c

所以

y\leq \lfloor \frac{c-ax}{b}\rfloor

那么对于一个固定的x，合法的y的个数为 \lfloor \frac{c-ax}{b}\rfloor+1,+1 是因为y可以取0

同时，y=0 时 x 最大，因此 x 的上限为\lfloor \frac{c}{a}\rfloor

因此等价于

\sum_{x=0}^{\lfloor \frac{c}{a}\rfloor}\lfloor \frac{c-ax}{b}\rfloor+1

不过因为 -a 小于0，所以不能直接上类欧的板子。

考虑加上一个 bx ,再在外面减掉一个 x

\sum_{x=0}^{\lfloor \frac{c}{a}\rfloor}\lfloor \frac{c+(b-a)x}{b}\rfloor-x+1

前面的直接用类欧，后边的可以直接算出来。

因为原题的限制对 a,b 是对称的，因此为了保证 b-a\geq0,如果 b<a,交换一下就好了。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
LL f(LL n,LL a,LL b,LL c)
{
    LL ac=a/c,bc=b/c;
    LL m=(a*n+b)/c;
    if(a==0)
    {
        return (n+1)*bc;
    }
    if(a>=c||b>=c)
    {
        return f(n,a%c,b%c,c)+n*(n+1)/2*ac+(n+1)*bc;
    }
    return n*m-f(m-1,c,c-b-1,a);
}
int main()
{
    LL a,b,c;
    cin>>a>>b>>c;
    if(a>b) swap(a,b);
    LL m=(c/a);
    LL res=f(m,b-a,c,b)-m*(m+1)/2+m+1;
    cout<<res;
    return 0;
}

P5172 Sum

设x=\sqrt r，注意这里没有取整

若x为整数时

若x是偶数，那答案为n

若x是奇数，那当n为奇数时答案为-1，为偶数时答案为0

若x不是整数

首先有一个小引理：

(-1)^n=1-2(n\bmod2)=1-2(n- 2\lfloor \frac{n}{2}\rfloor)=1-2n+4\lfloor \frac{n}{2}\rfloor

所以原式为

\sum_{i=1}^n1-2\lfloor ix\rfloor+4\lfloor \frac{\lfloor ix \rfloor}{2}\rfloor =n-2\sum_{i=1}^n\lfloor ix\rfloor+4\sum_{i=1}^n\lfloor \frac{ix}{2}\rfloor

设

f(n,a,b,c)=\sum_{i=1}^n\lfloor \frac{ax+b}{c} i\rfloor

那么题目所求即为，n-2f(n,1,0,1)+4f(n,1,0,2)

记 p=\frac{ax+b}{c},q=\lfloor p\rfloor

那么

f(n,a,b,c)=\sum_{i=1}^n\lfloor pi\rfloor =\sum_{i=1}^n\lfloor pi-qi+qi\rfloor =\sum_{i=1}^n\lfloor (p-q)i\rfloor+qi =\frac{n(n+1)}{2}q+\sum_{i=1}^n\lfloor (\frac{ax+b}{c}-q)i\rfloor =\frac{n(n+1)}{2}q+f(n,a,b-cq,c)

这样当 q 不为0的时候我们可以把 q 化为等于 0 的情况

f(n,a,b,c)=\sum_{i=1}^n\sum_{j=1}^{\lfloor pn \rfloor}[j<pi] f(n,a,b,c)=\sum_{i=1}^n\sum_{j=1}^{\lfloor pn \rfloor}[\frac{j}{p}<i] f(n,a,b,c)=\sum_{j=1}^{\lfloor pn \rfloor}\sum_{i=1}^n[\frac{j}{p}<i] f(n,a,b,c)=\sum_{j=1}^{\lfloor pn \rfloor}n-\lfloor \frac{j}{p} \rfloor f(n,a,b,c)=\sum_{j=1}^{\lfloor pn \rfloor}n-\lfloor \frac{j}{\frac{ax+b}{c}} \rfloor f(n,a,b,c)=\sum_{j=1}^{\lfloor pn \rfloor}n-\lfloor \frac{c}{ax+b}j \rfloor

考虑把 x 翻上去

f(n,a,b,c)=\sum_{j=1}^{\lfloor pn \rfloor}n-\lfloor \frac{c(ax-b)}{(ax+b)(ax-b)}j \rfloor f(n,a,b,c)=\lfloor pn \rfloor n-\sum_{j=1}^{\lfloor pn \rfloor}\lfloor \frac{acx-bc)}{a^2x^2-b^2}j \rfloor f(n,a,b,c)=\lfloor pn \rfloor n-\sum_{j=1}^{\lfloor pn \rfloor}\lfloor \frac{acx-bc)}{a^2r-b^2}j \rfloor f(n,a,b,c)=\lfloor pn \rfloor n-f(\lfloor pn \rfloor,ac,-bc,a^2r-b^2)

显然 n 是越来越小的，复杂度O(\log n)

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double db;
db x;
LL r; 
LL gcd(LL a,LL b)
{
    if(!b) return a;
    return gcd(b,a%b);  
}
LL f(LL a,LL b,LL c,LL n)
{
    if(n==0) return 0;
    LL d=gcd(a,gcd(b,c));
    a/=d;b/=d;c/=d;
    LL t=(a*x+b)/c;
    if(t!=0) return n*(n+1)/2*t+f(a,b-c*t,c,n);
    LL k=(a*x+b)*n/c;
    return n*k-f(a*c,-b*c,a*a*r-b*b,k); 
}
int main()
{
    int T;
    cin>>T;
    while(T--)
    {
        LL n;
        scanf("%lld %lld",&n,&r);
        x=sqrt(r);
        int s=(int)x;
        if(1ll*s*s==r)
        {
            if(r%2==0) printf("%lld\n",n);
            else if(n%2==0) printf("0\n");
            else printf("-1\n");
        }
        else printf("%lld\n",n-2ll*f(1,0,1,n)+4ll*f(1,0,2,n));
    }
    return 0;
}

P5179 Fraction

一道有意思的题目

首先有一个小引理：使得 p 最小等价于使得 q 最小

否则的话,设 p 最小时 q 为 x，q 最小时 p为 y 。

x>q,y>p \frac{a}{b}<\frac{p}{x}<\frac{p}{q}<\frac{y}{q}<\frac{c}{d}

矛盾，故命题成立。

接着分类讨论:

\lfloor \frac{a}{b}\rfloor+1\leq \lfloor \frac{c-1}{d}\rfloor

那么这两个分数之间有一个整数，此时 q 取到最小值 1，根据上述引理等价于 p 最小

那么就是

\frac{p}{q}<\frac{c}{d} q>\frac{dp}{c}

取 p=1, q=\lfloor \frac{d}{c} \rfloor+1

即为

\frac{a\mod b}{b}<\frac{p-q\lfloor \frac{a}{b} \rfloor}{q}< \frac{c-d\lfloor \frac{a}{b}\rfloor}{d}

递归即可

复杂度O(\log n)

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
void f(LL a,LL b,LL c,LL d,LL &p,LL &q)
{
    if(a/b+1<=(c-1)/d) p=a/b+1,q=1;
    else if(a==0) p=1,q=(d/c+1);
    else if(a<=b&&c<=d) f(d,c,b,a,q,p);
    else f(a%b,b,c-d*(a/b),d,p,q),p+=q*(a/b); 
}
int main()
{
    LL a,b,c,d;
    while(scanf("%lld %lld %lld %lld",&a,&b,&c,&d)!=EOF)
    {
        LL p=0,q=0;
        f(a,b,c,d,p,q);
        printf("%lld/%lld\n",p,q);
    }
    return 0;
}

[COCI2009-2010#1] ALADIN

首先可以把区间离散化，然后用线段树维护区间覆盖，因此，我们只需要要解决形如下列问题即可

\sum_{i=l}^r((i-L+1)\times A\mod B) \sum_{i=l}^r((i-L+1)\times A-B\lfloor \frac{(i-L+1)\times A}{B}\rfloor)

设

F(n)=\sum_{i=0}^n((i-L+1)\times A-B\lfloor \frac{(i-L+1)\times A}{B}\rfloor)

那么直接用 F(r)-F(l-1) 即可

前边可以 O(1)，看后面实际上等价于 B\times f(r,A,A-LA,B)，不过第二项有可能小于零，可以类似于上面处理的做法做，加上一个值再减掉一个值。

其他的用线段树维护即可

其他题目

CF1098E Fedya the Potter

CF1182F Maximum Sine

CF868G El Toll Caves