O(1) 查询gcd

我们来安利一个黑科技。（其实是Claris安利来的

比如我现在有一坨询问，每次询问两个不超过n的数的gcd。

n大概1kw，询问大概300w（怎么输入就不是我的事了，大不了交互库

http://mimuw.edu.pl/~kociumaka/files/stacs2013_slides.pdf

http://drops.dagstuhl.de/opus/volltexte/2013/3938/pdf/26.pdf

我们定义一个数k的一种因式分解k=k1*k2*k3为“迷之分解”当且仅当k1、k2、k3为质数或小于等于$\sqrt{k}$ 。

我们发现线筛的时候对于一个数x，设x最小的质因子为p，x/p=g，那么x的“迷之分解”可以通过g的“迷之分解”中三个数最小的一个乘上p得到。

证明似乎可以用数学归纳法证(然而我证不出来啊

然后对于每两个小于等于$\sqrt{n}$ 的数我们可以打一张gcd表出来。

最后如果我们要询问gcd(x,y)，我们找到x的“迷之分解”，然后如果分解的一部分小于等于$\sqrt{n}$ 那就查表，否则那就是一个质数，分类讨论一下就行了。

伪代码:

UPD:实际测试了一下随机数据跑得并没有沙茶gcd快。可能是我实现的姿势不够优越(雾

大家可以测试一下跑gcd(5702887,9227465)这个算法比沙茶gcd不知道快到哪里去了

//跑得比谁都快的gcd？

#include <iostream>

#include <stdio.h>

#include <stdlib.h>

#include <algorithm>

#include <string.h>

#include <vector>

#include <math.h>

#include <time.h>

#include <limits>

#include <set>

#include <map>

using namespace std;

const int N=;

const int sn=sqrt(N);

bool np[N+];

int ps[N+],pn=;

int cs[N+][];

void xs()

{

    np[]=cs[][]=cs[][]=cs[][]=;

    for(int i=;i<=N;i++)

    {

        if(!np[i]) {cs[i][]=cs[i][]=; cs[i][]=i; ps[++pn]=i;}

        for(int j=;j<=pn&&i*ps[j]<=N;j++)

        {

            np[i*ps[j]]=;

            int cm=cs[i][]*ps[j];

            if(cm<cs[i][])

            {

                cs[i*ps[j]][]=cm;

                cs[i*ps[j]][]=cs[i][];

                cs[i*ps[j]][]=cs[i][];

            }

            else if(cm<cs[i][])

            {

                cs[i*ps[j]][]=cs[i][];

                cs[i*ps[j]][]=cm;

                cs[i*ps[j]][]=cs[i][];

            }

            else

            {

                cs[i*ps[j]][]=cs[i][];

                cs[i*ps[j]][]=cs[i][];

                cs[i*ps[j]][]=cm;

            }

            if(i%ps[j]);else break;

        }

    }

}

int gcdd[sn+][sn+];

void smgcd()

{

    for(int i=;i<=sn;i++) gcdd[i][]=gcdd[][i]=i;

    for(int i=;i<=sn;i++)

    {

        for(int j=;j<=i;j++) gcdd[i][j]=gcdd[j][i]=gcdd[i-j][j];

    }

}

void pre_gcd() {xs(); smgcd();}

int gcd(int a,int b)

{

    if(a>N||b>N)

    {

        puts("Fuck You\n");

        return -;

    }

    int *x=cs[a],g=;

    for(int i=;i<;i++)

    {

        int d;

        if(x[i]<=sn) d=gcdd[x[i]][b%x[i]];

        else if(b%x[i]) d=;

        else d=x[i];

        g*=d; b/=d;

    }

    return g;

}

int euclid_gcd(int x,int y)

{

    while(y)

    {

        int t=x%y; x=y; y=t;

    }

    return x;

}

int tmd=-;

void gc()

{

    if(tmd==-) tmd=clock();

    else

    {

        printf("Passed: %dms\n",clock()-tmd);

        tmd=-;

    }

}

int main()

{

    int seed=time();

    //1kw个随机数测试

    int ans;

    printf("Euclid gcd...\n");

    srand(seed);

    gc();

    ans=;

    for(int i=;i<=;i++)

    {

        int a=(rand()*+rand())%N+,b=(rand()*+rand())%N+;

        ans^=euclid_gcd(a,b);

    }

    printf("Ans = %d\n",ans);

    gc();

    printf("New gcd...\n");

    srand(seed);

    gc();

    pre_gcd();

    ans=;

    for(int i=;i<=;i++)

    {

        int a=(rand()*+rand())%N+,b=(rand()*+rand())%N+;

        ans^=gcd(a,b);

    }

    printf("Ans = %d\n",ans);

    gc();

}

巴特西

O(1) 查询gcd

最新文章

热门文章