关于数的一些算法

最新推荐文章于 2021-05-16 21:38:04 发布

kalilili

最新推荐文章于 2021-05-16 21:38:04 发布

阅读量1k

点赞数 1

分类专栏：数论

本文链接：https://blog.csdn.net/kalilili/article/details/42420823

版权

数论专栏收录该内容

24 篇文章 1 订阅

订阅专栏

1.寻找一个数的所有质因数（整数分解）

不要一个个找因数，然后再判断是不是素数，这样太慢了，这是比较高效的算法思路是把某一num的素因数一个个剥离出来：（且从小到大存在vector里）

算法复杂度：最差情况：O（根号n）比如素数100000007，最优情况：O（logN）比如8，16，32，64，3，9，27，81...（因为num/i这步是根号n的复杂度向logN的转变）

vector<int> Prime_Factors;
for(i=2;i*i<=num;i++)              //i*i<num 和判断素数的循环条件 一样作用也相同，当num是素数或者1的时候结束循环。 
{
    if(num%i==0)			//只要num%i==0那么i一定是素数。若i不是素数，他的素因数比自身小，已经从前几次循环中从num中剥离掉了。 
    {
        Prime_Factors.push_back(i);
        while(num%i==0)
            num/=i;           //核心步骤 
    ｝
}
if(num>1) Prime_Factors.push_back(num);		//最后剩下的num如果不是1的话那么也是它的素因数。

也可以打一个根号n以内的素数表复杂度O（根号N）这样就不需要一个一个的i++了，直接用素数表的数当作i即可,一般题目是多组数据所以用素数表来做复杂度会更低。

精练一点的代码有：

map<int ,int> prime_factor(int n)
{
	map<int,int> res;
	for(int i=2;i*i<=n;i++){
		while(n%i==0){
			++res[i];
			n/=i;
		}
	}
	if(n!=1) res[n]=1;
	return res;
}

比如hduoj 的Alexandra and Prime Numbers题

2.寻找一个数n的所有约数

按照约数的性质，一个数n有一半的约数存在于根号n以内，另一半多的约数存在于n-根号n范围内，而且由一个已知的约数m，可以推出另一个约数n/m。所以可以遍历根号n以内的所有约数并且记录下来，再根据现有的推出其他的约数（代码是从小到大储存）复杂度O(根号n)

    cin>>n;
    int Divisor[10000];
    for(int i=1;i*i<=n;i++)
        if(n%i==0) Divisor[p++]= i;
    for(int j=p-1;j>=0;j--)
        if( Divisor[j]*Divisor[j]==n ) continue;
        else Divisor[p++]= n/Divisor[j];

精炼点的代码：

vector<int> divisor(n)
{
	vector<int >res;
	for(int i=1;i*i<=n;i++){
		if(n%i==0) {
			res.push_back(i);
			if(i!=n/i) res.push_back(n/i);
		}
	}
	return res;
}

3.寻找a，b两个整数的所有公约数

由数论知识可知：a和b的所有公约数都是他们的最大公约数的约数，所以可以先求出gcd 复杂度O（logN）再根据标题2求出gcd的所有约数复杂度O（根号（gcd(a,b)））

    cin>>a>>b;
    int n,p=1;
    n=gcd(a,b);
    int Divisor[10000];
    for(int i=1;i*i<=n;i++)
        if(n%i==0) Divisor[p++]= i;
    for(int j=p-1;j>=0;j--)
        if( Divisor[j]*Divisor[j]==n ) continue;
        else Divisor[p++]= n/Divisor[j];

比如 http://blog.csdn.net/kalilili/article/details/42640687（求a和b的第k大公约数）

4.素数表模版

素数代码：（简单易懂）

ll prime[N],primenum;//有primenum个素数 include<math.h>
void PRIME(ll Max_Prime){
	primenum=0;
	prime[primenum++]=2;
	for(ll i=3;i<=Max_Prime;i+=2)
	for(ll j=0;j<primenum;j++)
		if(i%prime[j]==0)break;
		else if(prime[j]>sqrt((double)i) || j==primenum-1)//||前面那个条件用于剪枝,||后面是必须的
		{
			prime[primenum++]=i;
			break;
		}
}

线性筛代码：（较难，效率高）

const int maxn = 100000003;
int p[6666666], tot;
bool vis[maxn];
void get_prime(){
    int i, j;
    for (i = 2;i < maxn; i++) {
        if (!vis[i]) p[++tot] = i;
        for (j = 1; j <= tot; j++){
            if (p[j] * i >= maxn)break;
            vis[p[j]*i] = 1;      //根据线性筛的理论知识知,应该有p[j]*(1...i)均标记,但这里只标记了p[j]*i来优化,原因自己思考 
            if (i % p[j] == 0)break; //如果i不是素数,标记完p[j]*i后就不需再标记p[j...]*i了,又是一个优化. 
        }
    }
}

上面的代码虽然优化得较好但是不易码出来，可以用易理解的模版：

int prime[MAX_N];
bool is_prime[MAX_N+1];
int seive(int n)
{
    int p=0;
    for(int i=0;i<=n;i++) is_prime[i]=true;
    is_prime[0]=is_prime[1]=false;
    for(int i=2;i<=n;i++)
        if(is_prime[i])
        {
            prime[p++]=i;
            for(int j=2*i(或i*i 注意防止爆int);j<=n;j+=i) is_prime[j] = false;
        }
    return p;
}

复杂度：n+n/2+n/3+n/4+...+n/n=nlnn 所以复杂度上界是O(nlnn) 上面那个线性筛的模版经过优化更多，复杂度会更优

打表代码：

int main(){  
    PRIME(10000);  
    FILE *fp=fopen("素数表","w");  
    for(ll i=0;i<primenum;i++)  
    fprintf(fp,"%lld,",prime[i]);  
    fclose(fp);  
    return 0;  
}

5.区间筛法

给定整数a,b 求区间[a,b)的素数，当a，b很大的时候不能直接用筛法，所以要空间和时间优化，先筛去[2，根号b)的合数，再向上取整筛去[a,b)的合数。

typedef long long ll;
bool is_prime[MAX_L];
bool is_prime_small[MAX_SQRT_B];
// 对区间[a,b)内的整数执行筛法,is_prime[i-a]=true <==> i是素数
void segment_sieve(ll a,ll b){
    for(int i=0;(ll)i*i<b;i++) is_prime_small[i]=true;
    for(int i=0;i<b-a;i++) is_prime[i]=true;
    if(1-a>=0) is_prime[1-a]=false;   //易错因为1不是素数也不是合数，这也是区间筛的一个易错bug
     
    for(int i=2;(ll)i*i<b;i++){
        if(is_prime_small[i]){
            for(int j=2*i;(ll)j*j<b;j+=i) is_prime_small[j] = false; //筛[2,根号b)
            for(int j=max((ll)2,(a+i-1)/i)*i;j<b;j+=i) is_prime[j-a]=false;
            // ((a-1)/i+1)*i，向上取整找到a以上的第一个i的倍数，这里易错，有可能a<=i,
            //这个时候a以上的第一个i的倍数就是i。这个时候不要把它筛去了，所以用到了max
        }
    }
}

6.扩展欧几里德算法

a 和 b 的最大公约数是 gcd ，那么，我们一定能够找到这样的 x 和 y ，使得: a*x + b*y = gcd 这是一个不定方程，有多解是一定的，但是只要我们找到一组特殊的解 x0 和 y0 那么，我们就可以用 x0 和 y0 表示出整个不定方程的通解：

x = x0 + (b/gcd)*t

y = y0 –(a/gcd)*t

假设当前我们要处理的是求出 a 和 b的最大公约数，并求出 x 和 y 使得 a*x + b*y= gcd ，而我们已经求出了下一个状态：b 和 a%b 的最大公约数，并且求出了一组x1 和y1 使得： b*x1 + (a%b)*y1 = gcd ，

由a%b = a - (a/b)*b（这里的 “/” 指的是整除）

得：gcd = b*x1 + (a-(a/b)*b)*y1 = b*x1 + a*y1 – (a/b)*b*y1 = a*y1 + b*(x1 – a/b*y1)

解：x = y1，y = x1 – a/b*y1

int e_gcd(int a,int b,int&x,int&y)
{
	if(!a&&!b) return -1;//无最大公约数
	int d=a;
	if(b!=0){
		ans=e_gcd(b,a%b,y,x); //把x,y颠倒，由于引用传参，所以在下一步就省了swap(x,y) 
		y-=(a/b)*x;    //x，y相当已经颠倒了，所以别把x，y写反了 
	}else{
		x=1;y=0;
	}
	return d;
}

一般写成这样即可：

ll e_gcd(ll a,ll b,ll&x,ll&y)  
{  
    ll ans;  
    if(b==0) ans=a,x=1,y=0;  
    else ans=e_gcd(b,a%b,y,x),y-=(a/b)*x;  
    return ans;  
}

相关题目 POJ 1061 http://poj.org/problem?id=1061 HDU 2669 http://acm.hdu.edu.cn/showproblem.php?pid=2669

7.求乘法逆元

对a∈Z*n，存在x∈Z*n，使得a×x ≡1 (mod n)，则称x为a的乘法逆元。可以看出a，n必须互质才有解。

可以等价于求这样的方程的解x： a*x + n*y = 1，所以乘法逆元有多组解，一般都是求最小的非负整数解，所以可以用扩展欧几里德解出一个特解x0，所以通解是 x = x0 + (n/gcd(a,n))*t gcd(a,n)=1，所以通解为x = x0 + n*t

所以求最小的非负整数解x即把特解x0模n即可，这里指数学上的模，在计算机上当x0是负数时%n还是负数，所以还要处理一步，当x0负数那么解为：x0模上 abs（m）再加上abs(m) ，于是，我们不难写出下面的代码求解:

//ax= 1(mod n)  --->  ax+ny=1
long long mod_reverse(long long a,long long n)  
{  
    long long x,y;  
    long long d=extend_gcd(a,n,x,y);  
    if(d==1) return (x%n+n)%n;  // 小技巧 
    else return -1;  
}

到此就结束了，下面稍稍浪一下，由求逆元的解法可以拓展到用扩展欧几里德求不定方程ax≡c(mod n)即ax+ny=c的解，同理如果c不能整除gcd(a,n)那么方程无解，先用扩展欧几里德求 a*x + n*y = gcd(a,n)的特解x0，y0,然后原方程的特解就为x1= x0*c/gcd(a,n)，y1=y0*c/gcd(a,n)。所以原不定方程的通解为 x = x1 + (n/gcd(a,n))*t， y = y1 –(a/gcd(a,n))*t