欧拉快速幂/欧几里得扩展算法/欧拉-费马降幂

最新推荐文章于 2022-10-06 01:38:45 发布

无情の学习机器

最新推荐文章于 2022-10-06 01:38:45 发布

阅读量743

点赞数 4

分类专栏：数据结构及算法笔记文章标签：算法 c++ 数学蓝桥杯

本文链接：https://blog.csdn.net/weixin_42034081/article/details/121724993

版权

数据结构及算法笔记专栏收录该内容

19 篇文章 0 订阅

订阅专栏

以下内容都是初等数学的知识…

欧拉快速幂

快速幂

如果大家稍微刷了点力扣应该都很容易能写出来：

using ll = long long;
ll& mypow(ll& a, ll& b) 
{
	ll res = 1;
	while (b) 
	{
		if (b & 1) res *= a;
		a *= a;
		b=b >> 1;
	}
	return res;
}

主要思想就是利用指数的二进制表达方式，比如7次方写成二进制就是1011，那么比如5的7次方就可以写作：

$5\times 5^{2}\times 5^{4}$

因此总体思路就是，每往前推一位，那么底数就相应的平方，如果当前位置上的二进制位为1，那么就乘上当前的底数。

复杂度是log2b，因为是位运算。

欧拉快速幂

当我们求a^b对c取模，其中abc都是整数。

当ab都不大的时候，上面的快速幂能够轻易得出结果，但是问题是，如果ab非常大，比如a=999999999，b=99999999，那么这个b次方数就根本不可能用long long表示出来。

根据以下性质：

$(a\times b)%c=[a%c\times b%c]%c$

我们可以将上面的快速幂改写成以下形式：

using ll = long long;
ll& mypow(ll& a, ll& b,ll&c) 
{
	ll res = 1;
	a %= c;
	while (b) 
	{
		if (b & 1) res = (a * res)%c;
		a = (a*a)%c;
		b=b >> 1;
	}
	return res;
}

时间复杂度一样是log2b

欧几里得(扩展)算法

欧几里得算法

先讲讲欧几里得算法，也就是我们常见的辗转相除法，大概意思就是求俩数的最大公约数，我们可以用以下代码简单实现：

using ll = long long;
ll& Mygcd(ll& a, ll& b) 
{
	if (a == 0) return b;
	if (b == 0) return a;
	if (a > b) 
	{
		a -= b;
		return Mygcd(a, b);
	}
	else
	{
		b -= a;
		return Mygcd(a, b);
	}
}

上面的代码是按照正常人逻辑写的，但是其实有更简便的写法，考虑到a b在运算过程中会不断减小，而且如果a比b大很多倍的话，会反复对b进行减法，我们不如合并成一个取余数的操作，同时在传入参数的时候交替传入a和b的值：

using ll = long long;
constexpr ll& Mygcd(ll& a, ll& b) 
{
	if (b == 0) return a;
	a %= b;
	return Mygcd(b, a);
}

裴蜀定理

在讲欧几里得扩展算法之前，我们需要讲一下裴蜀定理（英语：Bézout's identity）。

借用以下维基百科：

说明了对任何整数a、b 和m，关于未知数 x 和 y 的线性丢番图方程（称为裴蜀等式）：

$ax+by=m$

有整数解时当且仅当m 是a 及b 的最大公约数d 的倍数,证明略。

扩展欧几里得算法

该算法旨在求出裴蜀定理的整数解x和y。

联合欧几里得公式，我们容易得到以下公式：

$gcd(a,b)=ax+by$

$gcd(b,a%b)=bX+(a%b)Y$

$gcd(a,b)=gcd(b,a%b)$

结合上面三个公式及以下等式：

$a%b=a-b[\frac{a}{b}]$

注意此处相除代表向下取整

得到以下等式：

$a(x-Y)-b(X-y-[\frac{a}{b}]Y)=0$

因为要对于ab都成立，因此我们得到：

$x=Y$

$y=X-[\frac{a}{b}]Y$

因此我们发现x y与比a b小一点的数的XY之间是有关系的，因此很容易就想到递归，那么递归的终点就是b==0；此时, gcd(a,b)就是a，那么终点的状态就是x=1,y=0;

根据上面的推论，我们可以写出以下代码来求x与y：

constexpr void extend(ll a, ll b, ll& x, ll& y) 
{
	if (b == 0) {
		x = 1LL;
		y = 0LL;
		return;
	}
	extend(b, a%b, y, x);
	y -= (a / b) * x;
}

那么这段代码有什么用呢，其实就是在求如果等式右边是gcd（a，b）的时候的xy值，跟我们传入任何的x y都没有关系，比如我们传入a=24，b=50，那么她们最大的公约数是2，此时算出来的x=-2，y=1；那么对于是2的倍数的问题，比如是4，那么这个问题的一个解就是简单2倍数xy的值就是。

欧拉-费马降幂

欧拉函数

写作phi（x），或者 $\varphi (x)$ ,其中x是正数，表示0到x中与x互质的数的个数‘

我们会有引理，对于素数x：

$\varphi (x)=x-1$

$\varphi (i\times x)=p\times \varphi(i).....(i%x==0)$

$\varphi (i\times x)=(p-1)\times \varphi(i).....(i%x!=0)$

还有人说可以用这个公式来表示：

$\varphi (x)= x\times\prod (1-\frac{1}{ai})$ 其中ai是x 的质因子

我们可以用下面代码简单实现欧拉函数，找到因子了，就要除掉，一直除除到除不动为止，这样就能保证我们找到的i都是质因子

using ll = long long;
constexpr ll Euler(ll x) 
{
	ll res = x;
	for (ll i = 2LL; i < x / i; ++i) 
	{
		if (!(x % i)) 
		{
			res = res / i * (i - 1);
			while (!(x % i)) x /= i;
		}
	}
	if (x > 1) res = res / x * (x - 1);
	return res;
}

最后if语句是为了保证我们已经除完了n的所有的素因子，有可能还会出现一个我们未除的因子，如果结尾出现n>1 ,说明我们还剩一个素因子木有除。比如x=8，那么到if的时候就是x=1；但是如果x=4，那么到if的时候x=4；

欧拉定理

一个公式，当a与n互质的时候：

$a^{\varphi (n)}%n=1$

这个定理最大的应用就是来简化幂的模运算：

比如我们要计算一个数非常大 $a^{b\varphi (n)}$ 对某个数n 的模，如果恰好a与n互质，那么我们可以写成. $(a^{\varphi (n)})^{b}%n$ , 又因为欧拉公式，我们知道可以将括号里面的分别取模再相乘，因此最后结果为1.

对于任意的指数c，我们都可以将其拆成 $\varphi (n)\times b+d$ 的形式，最后只需要求 $a^{c}%n$ 即可。

比如计算 $7^{222}$ 的个位数，实际是求 $7^{222}$ 被10除的余数。7和10互质，且φ(10)=4。由欧拉定理知：

$7^{4}\equiv 1(mod 10)$

所以:

$7^{222}=(7^{4} )^{55}\times 7^{2}\equiv 7^{2}\equiv 9(mod10)$

欧拉扩展定理

所谓扩展就是推广到更一般的情况，如果a与n不互质那么有以下定理：

$a^{b}\equiv a^{b}(mod n)......(b<\varphi (n))$

$a^{b}\equiv a^{b%\varphi (n)+\varphi (n )}(mod n)......(b>=\varphi (n))$

我们直接上题：

力扣https://leetcode-cn.com/problems/super-pow/这个题怎么做呢，评论里面有一个很高赞的回答说是根据欧拉-费马降幂，a^b %c == a^(b%phi(c)) % c（c是素数） phi(c)是欧拉函数，表示小于c的和c互质的数的个数。这道题给定c是1137，所以可以提前求出phi(c) = 1140。你也可以直接求（自己想想怎么实现）

其实这种做法是错误的，因为1337并不是质数，如果底数为7或者是1337的倍数，那么结果就是错误的，本身题的测试案例就设计就有缺陷。

因此我们要不互质的情况单独提出来说：

using ll = long long;
class Solution {
    constexpr int Euler(int x)
    {
        ll res = x;
        for (ll i = 2LL; i < x / i; ++i)
        {
            if (!(x % i))
            {
                res = res / i * (i - 1);
                while (!(x % i)) x /= i;
            }
        }
        if (x > 1) res = res / x * (x - 1);
        return res;
    }
    const int mod = 1337;
    constexpr int mypow(int a, int b, int c)
    {
        ll res = 1;
        a %= c;
        while (b)
        {
            if (b & 1) res = (a * res) % c;
            a = (a * a) % c;
            b = b >> 1;
        }
        return res;
    }
    constexpr int Mygcd(int a, int b)
    {
        if (b == 0) return a;
        a %= b;
        return Mygcd(b, a);
    }

public:
    int superPow(int a, std::vector<int>& b) {
        if (b.empty()) return 1;
        int n = Euler(mod);
        auto f = [a,this](std::vector<int>& b, int n)
        {
            bool flag = false;
            int pownum = 0;
            for (int i = 0; i < b.size(); ++i)
            {
                pownum = pownum * 10 + b[i];
                if (pownum > n)
                {
                   pownum %= n;
                   flag = true;
                }
            }
            if (flag&&Mygcd(a,n)!=1) pownum += n;
            return pownum;
        };
        int pownum = f(b, n);
        return mypow(a, pownum, mod);
    }
};