zoukankan      html  css  js  c++  java
  • 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}$ 那就查表,否则那就是一个质数,分类讨论一下就行了。

    伪代码:

    image

    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=10000000;
    const int sn=sqrt(N);
    bool np[N+5];
    int ps[N+5],pn=0;
    int cs[N+5][3];
    void xs()
    {
        np[1]=cs[1][0]=cs[1][1]=cs[1][2]=1;
        for(int i=2;i<=N;i++)
        {
            if(!np[i]) {cs[i][0]=cs[i][1]=1; cs[i][2]=i; ps[++pn]=i;}
            for(int j=1;j<=pn&&i*ps[j]<=N;j++)
            {
                np[i*ps[j]]=1;
                int cm=cs[i][0]*ps[j];
                if(cm<cs[i][1])
                {
                    cs[i*ps[j]][0]=cm;
                    cs[i*ps[j]][1]=cs[i][1];
                    cs[i*ps[j]][2]=cs[i][2];
                }
                else if(cm<cs[i][2])
                {
                    cs[i*ps[j]][0]=cs[i][1];
                    cs[i*ps[j]][1]=cm;
                    cs[i*ps[j]][2]=cs[i][2];
                }
                else
                {
                    cs[i*ps[j]][0]=cs[i][1];
                    cs[i*ps[j]][1]=cs[i][2];
                    cs[i*ps[j]][2]=cm;
                }
                if(i%ps[j]);else break;
            }
        }
    }
    int gcdd[sn+3][sn+3];
    void smgcd()
    {
        for(int i=0;i<=sn;i++) gcdd[i][0]=gcdd[0][i]=i;
        for(int i=1;i<=sn;i++)
        {
            for(int j=1;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
    ");
            return -1;
        }
        int *x=cs[a],g=1;
        for(int i=0;i<3;i++)
        {
            int d;
            if(x[i]<=sn) d=gcdd[x[i]][b%x[i]];
            else if(b%x[i]) d=1;
            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=-1;
    void gc()
    {
        if(tmd==-1) tmd=clock();
        else
        {
            printf("Passed: %dms
    ",clock()-tmd);
            tmd=-1;
        }
    }
    int main()
    {
        int seed=time(0); 
        //1kw个随机数测试 
        int ans;
        printf("Euclid gcd...
    ");
        srand(seed);
        gc();
        ans=0;
        for(int i=1;i<=10000000;i++)
        {
            int a=(rand()*32768+rand())%N+1,b=(rand()*32768+rand())%N+1;
            ans^=euclid_gcd(a,b);
        }
        printf("Ans = %d
    ",ans);
        gc();
        printf("New gcd...
    ");
        srand(seed);
        gc();
        pre_gcd();
        ans=0;
        for(int i=1;i<=10000000;i++)
        {
            int a=(rand()*32768+rand())%N+1,b=(rand()*32768+rand())%N+1;
            ans^=gcd(a,b);
        }
        printf("Ans = %d
    ",ans);
        gc();
    }
  • 相关阅读:
    Shell: 定期存档日志文件
    canva实践小实例 —— 马赛克效果
    canvas API ,通俗的canvas基础知识(五)
    canvas实践小实例二 —— 扇形
    canvas API ,通俗的canvas基础知识(四)
    canvas实践小实例一 —— 画板工具
    canvas API ,通俗的canvas基础知识(三)
    canvas API ,通俗的canvas基础知识(二)
    canvas API ,通俗的canvas基础知识(一)
    JavaScript小实例:拖拽应用(二)
  • 原文地址:https://www.cnblogs.com/zzqsblog/p/5436775.html
Copyright © 2011-2022 走看看