辗转相除法、埃拉托色尼筛选法、牛顿迭代法证明与C++实现

时间：2021-01-23 04:21:46

1.辗转相除法

辗转相除法是用来计算两个数最大公约数的。

对于m，n求最大公约数，公式为: gcd(m,n) = gcd(n, m mod n)

证明：   

//最大公约数用e表示。

对于m=n时，显然e = m = n:

       m>n时，有m = kn + m mod n

令r = m mod n，则m = kn + r，若存在正整数d|n (|表示整除)，且d|r，则必有d|m，所以n，r的公约数也是m，n的公约数。

那么重复该公式：m1 = kn1 + r1

     m2 = kn2 + r2

                             ...

   因为n一直递减，那么该算法终止于n为0的时候，我们知道，任何正整数与0的最大公约数就是改整数。因为0可以整除任何数。

所以，当存在(m, 0)，m是n，r的最大公约数，同时也是m，n的最大公约数，得证。

代码如下，有递归和非递归版本：

#include <iostream>

int gcd_detail(int v1, int v2) 
{
    while(v2 != 0){ 
        int mod = v1 % v2; 
        v1 = v2; 
        v2 = mod;
    }   
    return v1; 
}

#ifdef _RECURSION_
int gcd_detail(int v1, int v2) 
{
    return v2 == 0 ? v1 : gcd_detail(v2, v1 % v2);
}
#endif

int gcd(int v1, int v2) 
{
    if(v1 < 0 || v2 < 0 || (v1 == 0 && v2 == 0)) 
        return -1; 
    
    return v1 > v2 ? gcd_detail(v1, v2) : gcd_detail(v2, v1);
}

int main()
{
    int res = gcd(60, 24);
   std::cout<<res<<std::endl;
    return 0;
}


//update 更相减损术

九章算术说：可半者半之，不可半者，副置分母、子之数，以少减多，更相减损，求其等也。以等数约之。'

举例：(60, 24)

1.60 - 24 = 36, 则转为求（36，24）

2.36 - 24 = 24，(36, 24)

3.36 - 24 = 12, (12, 24)调整为(24,12)继续

4.24 - 12 = 12, (12, 12)

5.12 - 12 = 0,那么最大公约数就是12了。

int gcd_detail(int v1, int v2) 
{
    while(v1 != v2){
        int t = v1 > v2 ? v1 - v2 : v2 - v1; 
        if(t < v2){
            v1 = v2; 
            v2 = t;
        }   
        else
            v1 = t;
    }   
    return v1; 
}


2.埃拉托色尼筛选法

我们前面求最大公约数其实还可以用中学知识求解。

如：(60, 24)

对于这两个数：60 = 2 * 2 * 3 *5

 24 = 2 * 2 * 2 * 3

此时相同的质数因子的乘积就是最大公约数：2*2*3 = 12。

可是，中学教我们的这种方法，又该如何实现呢，埃拉托色尼筛选算法就可以实现目的。

通过消去相应的合数，留下的就是质数。

代码如下：

#include <iostream>
#include <assert.h>
#include <vector>
#include <math.h>

std::vector<int> find_prime_number(int n)
{
    assert(n > 1); 
    
    std::vector<int> A(n+2, 0); 
    for(int p=2; p<=n; ++p)
        A[p] = p;
    
    for(int p=2; p<=sqrt(n); ++p){
        if(A[p] != 0){ 
            int j = p * p;
            while(j <= n){ 
                A[j] = 0;
                j += p;
            }   
        }   
    }   

    std::vector<int> L;
    for(int p=2; p<=n; ++p){
        if(A[p] != 0)
            L.push_back(A[p]);
    }   

    return L;   
}

int main()
{
    auto res = find_prime_number(25);

    for(auto i : res)
        std::cout<<i<<' ';
    std::cout<<std::endl;
    return 0;
}


3.牛顿迭代法

最常见的牛顿迭代法的应用就是来求平方根了。

以下内容转自：http://blog.chinaunix.net/uid-24118190-id-75267.html

求n的平方根，先假设一猜测值X0 = 1，然后根据以下公式求出X1，再将X1代入公式右边，继续求出X2…通过有效次迭代后即可求出n的平方根，Xk+1

辗转相除法、埃拉托色尼筛选法、牛顿迭代法证明与C++实现

先让我们来验证下这个巧妙的方法准确性，来算下2的平方根 (Computed by Mathomatic)

1-> x_new = ( x_old + y/x_old )/2 y (x_old + -----) x_old #1: x_new = --------------- 2 1-> calculate x_old 1 Enter y: 2 Enter initial x_old: 1 x_new = 1.5 1-> calculate x_old 2 Enter y: 2 Enter initial x_old: 1 x_new = 1.4166666666667 1-> calculate x_old 3 Enter y: 2 Enter initial x_old: 1 x_new = 1.4142156862745 1-> calculate x_old 10 Enter y: 2 Enter initial x_old: 1 Convergence reached after 6 iterations. x_new = 1.4142135623731 ...

可见，随着迭代次数的增加，运算值会愈发接近真实值。很神奇的算法，可是怎么来的呢? 查了下wikipedia和wolfram，原来算法的名字叫Newton’s Iteration (牛顿迭代法)。

下面是极其boring的数理介绍，不喜欢数学的言下之意也就是绝大部分人可以略过了。

简单推导

假设f(x)是关于X的函数:

辗转相除法、埃拉托色尼筛选法、牛顿迭代法证明与C++实现

求出f(x)的一阶导，即斜率:

辗转相除法、埃拉托色尼筛选法、牛顿迭代法证明与C++实现

简化等式得到:

辗转相除法、埃拉托色尼筛选法、牛顿迭代法证明与C++实现

然后利用得到的最终式进行迭代运算直至求到一个比较精确的满意值，为什么可以用迭代法呢?理由是中值定理(Intermediate Value Theorem):

如果f函数在闭区间[a,b]内连续，必存在一点x使得f(x) = c，c是函数f在闭区间[a,b]内的一点

我们先猜测一X初始值，例如1，当然地球人都知道除了1本身之外任何数的平方根都不会是1。然后代入初始值，通过迭代运算不断推进，逐步靠近精确值，直到得到我们主观认为比较满意的值为止。例如要求768的平方根，因为252 = 625，而302 = 900，我们可先代入一猜测值26，然后迭代运算，得到较精确值:27.7128。

回到我们最开始的那个求2次方根的公式，令x2 = n，假设一关于X的函数f(x)为:

f(X) = X2 - n

求f(X)的一阶导为:

f'(X) = 2X

代入前面求到的最终式中:

Xk+1 = Xk - (Xk2 - n)/2Xk

化简即得到我们最初提到的那个求平方根的神奇公式了:

辗转相除法、埃拉托色尼筛选法、牛顿迭代法证明与C++实现

 

这里也有一个牛顿迭代法很好的证明： http://www.cnblogs.com/javathread/archive/2011/12/26/2634731.html

代码实现如下：

#include <iostream>

bool cmp(double d1, double d2) 
{
    return d1 - d2 < 0.000001 ? true : false;;
}

int newtons_iteration(int n)
{
    if(n < 0)
        return -1; 
    
    double g = n;    //这个初值是为控制范围，任何数平方根不可能大于它自己，然后一步一步逼近f(0),求的结果 
    while(!cmp(g*g-n, 0.000001)){
        g = (g + n/g) / 2;
    }   
    return g;
}

int main()
{
    int res = newtons_iteration(25);
    std::cout<<res<<std::endl;
    return 0;
}

标签：

相关文章