Basic Mathematical Algorithm


1 最大公约数

  两个数的最大公约数就是其相同的因数的乘积,其常见的求法有质因数分解法、短除法、辗转相除法、更相减损法。

1.1 辗转相除法

  辗转相除法也就是欧几里算法,其原理为(a,b) = (a, a%b),其中(a,b)指a和b的最大公约数。

证明:
①设$d = gcd(a, b)$,则$a = k_1d$,$b = k_2d$
 由于$a\%b = a - b[a/b]$, 即$a\%b=a - kb$
 所以$a\%b = k_1d - k k_2 d = (k_1 - k k_2)d$
因而$d$也是$a\%b$的公约数。
②$d = gcd(a, b)$,说明$k_1$和$k_2$互质;
 若要证明gcd(a,b)等价于gcd(b,a%b),只需证明$k_2$和$(k_1 - k k_2)$互质;
 假设$k_2$和$(k_1 - k k_2)$不互质,设$gcd(k_2, k_1 - k k_2) = d_1$;
 此时有$k_2 = xd_1$和$k_1 - k k_2) = yd_1$;
 得$k_1 = (kx+y)d_1$,则说明$k_1$和$k_2$不互质,与前边矛盾,故假设错误,$k_2$和$(k_1-k k_2)$互质.
因而$d$为$a,b$的最大公约数,且为$b,a\%b$的最大公约数,即gcd(a,b)=gcd(b,a%b).

/*
*   @param:  所要求最大公因数的两个常数a、b
*   @output: a和b的最大公因数
*/
int gcd(int a, int b) {
    if(!b)  return a;
    return gcd(b, a%b);
}

1.2 更相减损法

/*
*   @param:  所要求最大公因数的两个常数a、b
*   @output: a和b的最大公因数
*/
int gcd(int a, int b) {
    if(a == b)  return a;
    else if(a > b)  return gcd(a-b, b);
    return gcd(a, b-a);
}

1.3 Stein算法

/*
*   @param:  所要求最大公因数的两个常数a、b
*   @output: a和b的最大公因数
*   @note:   结合了辗转相除法和更相减损法的优点
*            a和b同为奇数,则进行一次更相减损;
*            一个为奇数则将偶数方/2;
*            均为偶数则等于二者均/2的公因数的2倍。
*/
int gcd(int a, int b) {
    if(a < b)       return gcd(b, a);
    if(a%b == 0)    return b;
    if(a%2 && b%2)  return gcd(a-b, b);     // 均为奇数,进行更相减损
    else if(a%2)    return gcd(a, b>>1);    
    else if(b%2)    return gcd(a>>1, b);
    return 2*gcd(a>>1, b>>1);
}

2 最小公倍数

  a和b的最小公倍数和最大公因数存在下列关系:

  不难写出:

/*
*   @param: 所要求最小公倍数的两个常数a、b
*   @output: a和b的最小公倍数
*/
int lcm(int a, int b) {
    return a * b / gcd(a, b);
}

3 快速幂与快速乘

3.1 快速幂

  • 递归实现
/*
*   @param: num为底,power为幂
*   @output: num的power次幂
*/
int myPow(int num, int power) {
	if(power == 0)	return 1;
	if(power % 2 == 1)	return num * myPow(num, power-1);
	return myPow(num*num, power/2);
}
  • 迭代实现
/*
*   @param: num为底,power为幂
*   @output: num的power次幂
*/
int myPow(int num, int power) {
	int ans = 1, multi = num;
	while(power) {
		if(power % 2 == 1) {
			ans *= multi;
		}
		power >>= 1;
		multi *= multi;
	}
	return ans;
}

3.2 快速乘

  • 递归实现
/*
*   @param: num1和num2对应两个乘数
*   @output: 输出num1和num2的乘积
*/
int myMulti(int num1, int num2) {
	if(num2 == 0)		return 0;
	else if(num2 == 1)	return num1;
	if(num2 % 2)	return num1 + myMulti(num1, num2-1);
	return myMulti(num1+num1, num2>>1);
}
  • 迭代实现
/*
*   @param: num1和num2对应两个乘数
*   @output: 输出num1和num2的乘积
*/
int myMulti1(int num1, int num2) {
	int sum = 0, add = num1;
	while(num2){
		if(num2 % 2) {
			sum += add;
		}
		add += add;
		num2 >>= 1;
	}
	return sum;
}

4 求根

  • 二分法求解
double nthRoot(int num, int nth) {
    double left = 0, right = num;
    double mid = left + ((right - left) / 2);
    static double epsilon1 = 1e-7;
    while(abs(myPow2(mid, nth) - num) > epsilon1) {
        if(myPow2(mid, nth) == num) 
            return mid;
        else if(myPow2(mid, nth) < num) 
            left = mid;
        else 
            right = mid;
        mid = left + ((right - left) / 2);
    }
    return left;
}
  • 牛顿迭代

  根据牛顿迭代法可知,对于函数$f(x) = x^n - C$,其零点可以通过下面的迭代求得:

当$abs(x_{n+1} - x_{n})<\epsilon$,其中$\epsilon$为很小的常数,则可以认为$x_{n+1}$为所要求的根。

/*
*   @param: num为待解参数,nth为根的次数
*   @output: num的nth次方根,double类型
*/
double nthRoot(int num, int nth) {
    static double epsilon = 1e-9;
    double res = 1.0;
    while(abs(mypow(res, nth) - num) > epsilon) 
        res = ((n-1)*res + num / myPow(res, n-1)) / n;
    return res;
}

5 素数

5.1 判断素数

/*
*   @param: num为待判断常数
*   @output: num为素数则返回true, 否则返回false
*/
bool isPrime(int num) {
    for(int i = 2; i*i < num; ++i) 
        if(num % i == 0)    return false;
    return true;
}

5.2 筛选素数

  • 埃氏筛选法
/*
*   @param: maxn选定0~maxn范围并筛选出该范围内的素数
*   @output: 0~maxn(不包括maxn)范围内的素数, 为vector类型
*/
vector<int> Iprime(int maxn) {
    vector<bool> isPrime(maxn, true);
    vector<int> ret;
    isPrime[0] = isPrime[1] = false;
    for(int i = 2; i < maxn; ++i) {
        if(isPrime[i]) {
            for(int j = 2; j*i < maxn; ++j) {
                isPrime[j*i] = false;
            }
            ret.emplace_back(i);
        }
    }
    return ret;
}
  • 欧拉筛选法
/*
*   @param: maxn选定0~maxn范围并筛选出该范围内的素数
*   @output: 0~maxn(不包括maxn)范围内的素数, 为vector类型
*/
vector<int> Iprime(int maxn) {
    vector<int> ret;
    vector<bool> isPrime(maxn, true);
    isPirme[0] = isPrime[1] = false;
    for(int i = 2; i < maxn; ++i) {
        if(isPrime[i])  ret.emplace_back(i);
        for(int j = 0; j<ret.size() && i*ret[j]<maxn; ++j) {
            isPrime[i*ret[j]] = false;
            if(i%ret[j] == 0)    break;
        }
    }
}

两种筛选法的不同:埃氏筛选法是筛选当前素数的所有倍数,因而可能会出现重复筛出,如24会被2,3所标记;而欧拉筛选法中,合数只会被其最小质因数标记。

6 概率问题

6.1 拒绝采样

  当我们可以通过一种采样方式获取某个区间的值,而我们实际上希望获取的值可能范围是比前边的采样方式所能获得的范围更小,那么如果我们同样希望在小范围中的采样同样也是等概率的,可以直接通过上述的采样方式进行采样,而当采样得到的数据不在我们期望的范围内,那么就拒绝这一次采样;如果落在我们的期望范围内则接收此次采样,这样就可以实现目标范围的等概率采样。

  举一个例子,当我们可以等概率采样的范围为[1, 10],而期望的数值范围是[3, 7]。 那么当采样得到的数据<3>7时拒绝这次采样,采到[3,7]的时候就接收,此时[3,7]中的数字是等概率的。

6.2 从randx()到randy()

  以rand2()为例,我们希望得到rand4()(每个数字采样概率均相等),那么可以通过(rand2()-1)*2+rand2(), 其中:

  • rand2()-1的结果为:0,1
  • (rand2()-1)*2的结果为:0,2
  • rand2()的结果为1,2
  • (rand2()-1)*2+rand2()的结果为1,2,3,4

  推广到从randx(), 可以有:(randx()-1)*x+randx()则可以获得[1,x^2]的值,而如果要获得randy(),那么就加上拒绝采样即可实现。


文章作者: Vyron Su
版权声明: 本博客所有文章除特別声明外,均采用 CC BY 4.0 许可协议。转载请注明来源 Vyron Su !