随机化算法小结（Miller Rabin,Pollard Rho, 模拟退火, 随机化贪心）
大纲:
算法正文
\(Miller\) - \(Rabin\)
(这是随机化算法)
算法用途： 快速判断单个数是否为质数。
算法核心: 费马小定理 以及 二次探测定理。
费马小定理: 设 \(p\) 为质数，假若 \(a\) 不为 \(p\) 的倍数，则 $a^{p-1} ≡ 1 (mod $ \(p)\) (证明在OIWIKI有，此处略)
那么是不是说，假如现在有一个数 \(x\) 使得多个与 \(x\) 不互为倍数关系的数 \(a\) 满足 $a^{x-1} ≡ 1 (mod $ \(x)\) 就是质数了呢？ 小范围内，很难举出反例，但是后来，人们找出了反例: 卡米切尔数。
但是要完全否定这个办法吗？不是的。
于是后来有了 二次探测定理 优化了上面的判定方法，使得判定的出错率降到几乎为 \(0\)。
二次探测定理:
对于素数 \(p\)，那么对于任意的 \(a\) 均有 \(a^2 ≡ 1 (mod\) \(p)\) ，那么 \(a = p - 1\) 或者 \(a = 1\)。
证明:
∵ \(a^2 ≡ 1 (mod\) \(p)\)
∴ \(a^2 - 1 ≡ 0 (mod\) \(p)\)
∴ \((a + 1)(a - 1) ≡ 0 (mod\) \(p)\)
∴ \(a = p - 1\) 或者 \(a = 1\)
这个定理的用处在于，这提供了一个新的思路（也就是 \(Miller-Rabin\) 算法的核心）：
先约定，假设现在待测的数是 \(p\) 并且其不是偶数（不然可以直接判）
那么有 \(p - 1\) = \(s * 2^r\)。
任意随机一个数 \(a\)，然后进行判断，先通过快速幂求出 $q = $ \(a^s\)，然后不断的对 \(q\) 进行平方运算，当 \(q = 1 (mod\) \(p)\) 的时候，判断 \(q\) 在平方之前的那个数是否等于 \(p - 1\) 或者 \(1\)，如果不等于，说明这个数 \(p\) 一定是合数。
是不是这么判过就万无一失了呢？不是的。这样只可以判断出它可能不是合数，但是并不能代表它一定是素数。最后还需要判断一下 $a^{p-1} ≡ 1 (mod $ \(p)\) 是否成立，这样子又可以判断它是否为素数。
上面这个算法貌似会出错，对，确实会出错。但是出错的概率大约为：\(\frac{1}{4}\) （有人貌似在论文里面证明了出错率似乎比这个更小？见文末扩展资料部分）
证明就略了，因为我也不会。因为单次测试的错误概率是 \(\frac{1}{4}\)，那么取不同的 \(a\) 进行 \(10\) 次上面测试的话，错误的概率就到了 \(10^{-6}\)，如果你追求更高的正确率的话，可以进行更多次的测试。
算法复杂度:
空间复杂度上并不需要辅助数组，可以视作 O\((1)\)。
时间复杂度上，不难看出单次测试的复杂度是 O\((logn)\) 的，同时为了避免出错，一般大概需要使用 \(10\) 次的测试。那么完成一个数的素数判定的复杂度大概是 O\((logn * 10)\) (一般只在需要判断的素数达到 \(10^9\)量级的时候或者需要多次判断的时候使用，不然还不如 O(\(\sqrt{n}\)) 的暴力判断来得直接)
算法模板:
下面的代码仅供参考。
#include <bits/stdc++.h>
using namespace std;
#define int long long
int Mod,M[20] = {0,2,3,5,7,11,13,17,19,23,31,61};
int quickmul(int x,int y) { //快速乘，防止爆 long long
int op = x , Ans = 0;
while(y) {
if(y & 1ll) Ans = (Ans + op) % Mod;
op = (op << 1) % Mod;
y >>= 1ll;
}
return Ans;
}
int quickpower(int x , int y) {
int op = x , Ans = 1;
while(y) {
if(y & 1ll) Ans = quickmul(Ans,op) , Ans %= Mod;
op = quickmul(op,op) , op %= Mod;
y >>= 1ll;
}
return Ans;
}
bool Miller_Rabin(int x) {
Mod = x;
int r = 0 , p = x - 1;
for(int i = 1 ; i <= 11 ; i ++) if(M[i] == x) return 1;
if(x % 2 == 0 || x == 1) return 0;
while(p % 2 == 0) p >>= 1ll , r ++;//将 p 化为 "2^r * s" 的形式，最后余下来的 p 即是对应的 "s"
for(int i = 1 ; i <= 11 ; i ++) {
int B = quickpower(M[i], p);
for(int j = 1 ; j <= r ; j ++) {
int K = (quickmul(B, B)) % Mod;
if(K == 1ll && B != 1ll && B != Mod - 1) return 0; //根据二次探测定理发现不是素数
B = K;
}
if(B != 1ll) return 0; //通过费马小定理判断一定不是素数
}
return 1;
}
补充说明：
对于测试底数 \(a\) 的选取可以是随机的，也可以是自己给定的一些素数，下面给出一些素数：
\(2,3,5,7,11,13,17,19,61\) ，用这 \(9\) 个素数测完的话, \(10^{18}\) 内的数测试是没有问题的（测试已过）。
同时需要注意的是，对于达到 \(2 * 10^9\) 以上的数进行验证的话，运算过程中可能会导致 \(long\) \(long\) 甚至 \(unsigned\) \(long\) \(long\) 溢出，所以需要使用 快速乘* 来防止爆 \(long\) \(long\)，这个没有算在时间复杂度里面，实际使用的时候要注意这个对于复杂度的影响。
快速乘的介绍补充在下面的 扩展资料 部分了。
\(Pollard\) \(Rho\) 算法介绍:
这也是随机化算法。
算法用途: 快速分解一个大数的质因子。
算法核心: 欧几里得辗转相除法，\(Pollard\) \(Rho\) 的随机函数。
引入 ：
现在有一个极大的数 \(n\) (\(1 <= n <= 10^{18}\)) ，你需要找出其的一个非平凡因子。
为什么要找出其中的一个非平凡因子（假设为 \(f\)）？因为如果能找出 \(f\)，那么就可以递归的分别处理 \(f\) 因子以及 \(\frac{n}{f}\) 的因子，这样子就能做到将 \(n\) 整个数因数分解。
解决本问题其中一种就是 O\((\sqrt{n})\) 级别的暴力分解法。当然，这个方法貌似在 \(10^{18}\) 的数据范围面前显得弱小又无助。那么这个问题是无解的吗？并不是。
算法正文:
有一个重要的性质：
- 任意两个数的 \(gcd\) 一定是两个数的因子。
现在相当于找到任意的一个正整数 \(k <= \sqrt{n}\) 使得 \(gcd(k,n) > 1\) 即可以解决一开始的问题。
那么随机找到一个正整数 \(k <= \sqrt{n}\) 使得 \(gcd(k,n) > 1\) 的概率是多大呢？假如 \(n\) 由 \(10^9\) 左右大小的一个质数平方得到，那么正确的概率就会变为 \(\frac{1}{10^9}\)。这微乎其微的正确概率，真的可以优化吗？
首先来了解一个有趣的东西（大家都知道吧？）：
假设指定一个日期，比如 \(2\) 月 \(14\) 日，那么随机从人群里面拉出一个人来，出现生日是这一天的人的概率 \(P = \frac{1}{365}\)，无论取多少次，每一次的概率都是 \(P = \frac{1}{365}\) 的。
但是，假设一个房子里面存在的人数为 \(23\) 个，那么这个房子中出现两个人是同一天生日的概率为 \(50\)%
当人数达到 \(60\) 个，那么房子中出现两个人是同一天生日的概率达到了 \(99.99\) %
(这个在 OIWIKI 上有严格的证明，不再赘述。)
这启发了，思考一下，到底是什么让 生日同一天 的概率大大增加？
来感性理解一下，对于第一个例子，每次都要使得那个人的生日在 \(2\) 月 \(14\) 日，这很难。
但是对于第二个例子，每次只要使得当前的人跟前面的人中的某一个的人生日相同即可，随着人数的增加，这个概率在不断的增大。
\(Pollard\) \(Rho\) 的随机数生成函数
于是想到构造一个伪随机数序列。 \(Pollard\) 发明了一个生成伪随机数的函数:
\(f(x) = (x ^ 2 + c)\) \(mod\) \(n\) (\(n\) 就是待找因数的 \(n, c\) 是随机的一个常数)。
它生成序列的方式为：\(x_i = f(x_{i - 1})(x_0\) 是一个随机的数 )
之所以说这个是个伪随机数序列是因为这个函数生成的随机数序列是由循环节的：
例如：\(x = 1\), \(c = 5\), \(n = 11\) 的时候：
\(1, 6, 8, 3, 3, 3, 3 .....\)
然后就全是 \(3\) 了，如果将这些数如下图一样排列起来，会发现这个图像酷似一个\(\rho\)，算法也因此得名 rho。
(摘自OI Wiki)
根据生日悖论来看的话，这个函数生成的随机数数量约为 O\((\sqrt{n})\) 个。那么假定一个 \(M\) 是 \(n\) 的最小非平凡因子（\(1 < M <= \sqrt{n}\)）,现在已经用 Pollard-Rho 随机函数生成了一个序列为 \(x_1\), \(x_2\),\(x_3\)...\(x_k\)。
考虑生成 \(y\) 序列， \(y_i = x_i\) \(mod\) \(M\)，还是根据生日悖论，可以知道大约生成 O\((\sqrt{M})\)(也就是约 O(\(n^{\frac{1}{4}}\))) 个数就会出现 \(y_i = y_j\) 的情况。那么假如存在 \(y_i = y_j\) 并且 \(x_i != x_j\) ，则有 \(M |\) \(|x_i - x_j|\)，则有 \(gcd(|x_i - y_i|, n) > 1\)，这时候就找到了一个 \(n\) 的非平凡因子。
需要注意的是，并不总是期望能在一次判断中直接就找出 \(n\) 的一个非平凡因子，当然会出现失败的情况，所以需要不断的调整 \(x_0\) 以及 \(c\),多次判断才能找出 \(n\) 的非平凡因子，但是因为每次判断所需要的时间复杂度约为 O\((n^{\frac{1}{4}})\) ,多次判断下来期望复杂度是: O\((n^{\frac{1}{4}})\) 的。
Floyd 判环优化
\(example:\)
假设两个人在赛跑，A 的速度快，B 的速度慢，经过一定时间后，A 一定会和 B 相遇，且相遇时 A 跑过的总距离减去 B 跑过的总距离一定是圈长的 n 倍。 -----节选自 OI Wiki
一开始，让 \(a = f(1)\), \(b = f(f(1))\)，接下来的每一次中，都令 \(a = f(a)\), \(b = f(f(b))\)，然后每次判断 \(gcd(|a - b|, n)\) \(> 1\) 是否成立, 在出现 \(a = b\) 的情况的时候，就相当于出现了环。
为什么 \(a = b\) 的情况出现了的话就相当于出现了环了呢？
首先因为 \(f(1) = (1^2 + c)\) \(mod\) \(n\)，那么一开始 \(a\) 相当于 \(x_1\),则 \(b\) 则是 \(x_2\),随着不断的更改，\(a\) 就会变成 \(x_2\),对应的 \(b\) 就会变成 \(x_4\), \(a\) 变成 \(x_3\) 的时候， \(b\) 就会变成 \(x_6\) .......重复不断下去，因为 \(b\) 始终比 \(a\) 快一倍，当 \(a,b\) 相等就相当于 \(b\) 把 \(a\) 套圈了，于是就出现了环。
参考代码:
LL Pollard_Rho(LL x) {
if(x == 4) return 2;
LL c = rand() % (x - 1) + 1, a = f(0, c, x), b = f(f(0, c, x), c, x), d = 1;
while(a != b) {
d = gcd(Abs(a - b), x);
a = f(a, c, x);
b = f(b, c, x); b = f(b, c, x);
if(d > 1) return d;
}
return x;
}
倍增优化判环:
每次调用 \(gcd\) 函数的话时间复杂度是 O\((logn)\) 的，但是有一个性质：
\(gcd(a,b) > 1\) ,那么对于任意正整数 \(c\) ， \(gcd(a*c,b) > 1\) 一定成立。
在 \(Floyd\) 优化判环的部分里面，每次都用 \(gcd(|a - b|, n)\) 来判断是否找到了 \(n\) 的非平凡因子，一个优化就是将这些 \(|a-b|\) 在 \(mod\) \(n\) 下累乘起来，然后每隔 \(2^7 - 1\) 次就判断一次。同时需要注意的是，如果累乘的过程中出现了累乘的结果 \(mod\) \(n\)为 \(0\)，那么就返回分解失败，重新调整 Pollard-Rho 里面的 \(c\) 来重新判断。
参考代码:
LL Pollard_Rho(LL x) {
LL c = rand() % (x - 1) + 1, a = 0, b = 0, d = 1;
for(int lim = 1; ; lim <<= 1, b = a) {
LL S = 1;
for(int i = 1 ; i <= lim ; i ++) {
a = f(a, c, x);
S = S * Abs(a - b) % x; //累乘起来，因为我开了 int128，所以敢这么大胆，正常下还是要写快速乘的
if(i % 127 == 0) {
d = gcd(S, x);
if(d > 1) return d;
}
}
LL d = gcd(S, x);
if(d > 1) return d;
}
return x;
}
为什么参考代码里面没有判断 \(S\) 累乘变成 \(0\)?
LL gcd(LL A, LL B) {
if(B == 0) return A;
return gcd(B, A % B);
}
这是写的 \(gcd\),倘若 \(S\) 变成 \(0\) 了，在做 \(gcd\) 的时候就会返回 \(x\) 也就是待判断的数，从而终止循环。
算法复杂度分析
空间复杂度很容易看出来是 O\((1)\) 的。
时间复杂度上，期望取 \(n^{\frac{1}{4}}\) 个数就可以跑出答案，所以期望时间复杂度为 O\((n^{\frac{1}{4}})\)，但是要说明的是，当要分解一个数的所有质因数的时候，倘若那个数的质因数比较多，那么 Pollard Rho 算法的时间复杂度将会被卡，也就是 O(被卡死)，要注意这一点。
补充说明:
跟 \(Miller-Rabin\) 算法一样，基本上只应用于大整数的因数拆分，比较小的数进行因数拆分是没有必要的，同时要注意的是程序的常数优化，不然可能会跑得很慢。要求至少掌握倍增优化判环和 \(Floyd\) 判环，能熟练的写出来。
例题：
题意简述
给定一个正整数 \(m\) 以及长度为 \(m\) 的正整数序列 \(a\) ,同时给出 \(n = \prod_{i=1}^{m}{a_i}\)。
你需要找出一个最大的 \(k\) 使得存在一个 \(d\) , \(d > 1\) 并且 \(d^k | n\)。求这个最大的 \(k\) 以及在 \(k\) 最大的情况下有多少个 \(d\) 满足条件。
\(m <= 600, 1 <= a_i <= 10^{18}\)
思路
是不是很板子？对没错就是板子。
关于这道题目，实际上就是让我们找出出现次数最多的质因子以及跟它出现次数相同的质因子的个数（假设为 \(w\)），它出现的次数为 \(k\) 也就是答案为 \(k\) ，同时那么存在的 \(d\) 的个数就是 \(2^w - 1\)。这毫无疑问。
然后我们现在要求出出现次数最多的质因子出现了多少次。那么很显然的 \(10^{18}\) 的数据范围无法使用 \(\sqrt{n}\) 级别的算法来做这一道题目。于是我们用 Pollard Rho 和 Miller Rabin 算法来优化这道题目。
但是，因为 Pollard Rho 如果要去分解一个大数的所有质因子的话，倘若这个大数的质因子很多，那么 Pollard Rho 的时间复杂度将无法保证。考虑怎么优化这个过程。
Pollard Rho 算法其实在使用的时候有一个小 \(trick\)，如果一个数含有很小的质因子，那么没有必要去用 Pollard Rho 算法去找这些很小的因子，反而不如直接先用暴力跑出这些小的质因子。于是定义 "很小的质因子"为小于等于 \(10^6\) 的质因子 *，对于给定的 \(a_i\) ，我们先用 \(10^6\) 以内的所有质数去判断一下看这个数是不是含有 \(10^6\) 以内的质因子，如果有，我们就直接将其除掉，这样操作后，很容易知道剩下的质因子不超过 \(2\) 个。
注：并非所有的题目都是定义为 \(10^6\) 以内的质因数为“很小的质因数”，根据题目而定，这里这么定义是因为 \(m\) 只有 \(600\) ,而 \(10^6\) 以内的质数不超过 \(80000\) 个，所以可以这么定义。
然后再使用 Pollard Rho 算法进行分解质因数的话就可以通过本题。
同时 Pollard Rho 用的时候一般都要搭配 Miller Rabin 素性测试算法，具体使用见下面代码的
fac(x) 函数里面。
另外值得注意的是，因为本题的答案中 \(2^w - 1\) 可能会达到 \(2^{600}\) ，所以需要手写高精度。
下面给出代码:
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define ll long long
#define ull unsigned long long
const int MAXN = 1e3;
int m, Mod;
int A[MAXN], M[25] = {0,2,3,5,7,11,13,17,19,23,31,61}; //M 即 Miller Rabin 探测底数
int book[1 << 20], Prime[1 << 20], tot = 0;
int P[MAXN << 10], cnt = 0, len = 0;
char s[114514];
map <int, int> mp;//开一个 map 来统计每个质因子出现的次数
inline ll mul(ll x,ll y,ll p){ //快速乘，防止爆long long
ll z = (long double)x / p * y;
ll ans = (ull)x * y - (ull)z * p;
return (ans + p) % p;
}
int f(int x, int c, int MOD) {return (mul(x, x, Mod) + c) % MOD;} //这个是 Pollard Rho 伪随机函数
int gcd(int A, int B) {return B == 0 ? A : gcd(B, A % B);}
int Abs(int x) {return x < 0 ? -x : x;}
int quickpower(int x, int y) {
int ans = 1, op = x;
while(y) {
if(y & 1) ans = mul(ans, op, Mod);
op = mul(op, op, Mod), op %= Mod;
y >>= 1;
}
return ans;
}
bool Miller_Rabin(int x) { //Miller Rabin 素性测试模板
Mod = x;
for(int i = 1 ; i <= 11 ; i ++) if(M[i] == x) return 1;
if(x % 2 == 0 || x == 1) return 0;
int p = x - 1, r = 0;
while(p % 2 == 0) p >>= 1, r ++;
for(int i = 1 ; i <= 10 ; i ++) {
int B = quickpower(M[i], p);
for(int j = 1 ; j <= r ; j ++) {
int K = mul(B, B, Mod);
if(K == 1 && B != Mod - 1 && B != 1) return 0;
B = K;
}
if(B != 1ll) return 0;
}
return 1;
}
int Pollard_Rho(int x) { // Pollard Rho 模板
if(x == 4) return 2;
int c = rand() % (x - 1) + 1, a = 0, b = 0, d = 1;
for(int lim = 1; ; lim <<= 1, b = a) { //采用了倍增优化
int S = 1;
for(int i = 1 ; i <= lim ; i ++) {
a = f(a, c, x);
S = mul(S, Abs(a - b), x); //累乘
if(i % 127 == 0) { // 隔 127 次判一下
d = gcd(S, x);
if(d > 1) return d;
}
}
int d = gcd(S, x);
if(d > 1) return d;
}
return x;
}
void GetPrime() { //线性筛
for(int i = 2 ; i <= 1e6 ; i ++) {
if(!book[i]) Prime[++ tot] = i;
for(int j = 1 ; j <= tot && Prime[j] * i <= 1e6; j ++) {
book[Prime[j] * i] = 1;
if(i % Prime[j] == 0) break;
}
}
return ;
}
void fac(int x) {
if(x == 1) return ;
if(Miller_Rabin(x)) {
if(mp[x] == 0) P[++cnt] = x; //保证每一个质因数只进队一次，防止重复统计
mp[x] ++;
return ;
}
int p = x;
while(p >= x) p = Pollard_Rho(x);
int flag = Miller_Rabin(p);
while(x % p == 0) {
if(mp[p] == 0 && flag) P[++ cnt] = p; //防止漏掉质因子
if(flag) mp[p] ++;
x /= p;
}
fac(x);
if(!flag) fac(p);//如果 p 是质数就没必要进行下去了，已经被丢进质因数队列中了
return ;
}
void Times() { //高精度乘法
int len = strlen(s + 1);
reverse(s + 1, s + 1 + len);
int add = 0;
for(int i = 1 ; i <= len ; i ++) {
int x = s[i] - '0';
x *= 2;
x += add; add = 0;
if(x >= 10) add = 1, x %= 10;
s[i] = x + '0';
}
if(add != 0) s[++ len] = '1';
reverse(s + 1, s + 1 + len);
return ;
}
signed main() {
cin >> m;
for(int i = 1 ; i <= m ; i ++) cin >> A[i];
sort(A + 1, A + 1 + m);
GetPrime();
for(int i = 1 ; i <= m ; i ++) { //这里是再处理较小的质因数
for(int j = 1 ; j <= tot ; j ++)
while(A[i] % Prime[j] == 0) {
if(mp[Prime[j]] == 0)
P[++ cnt] = Prime[j];
A[i] /= Prime[j];
mp[Prime[j]] ++;
}
}
for(int i = 1 ; i <= m ; i ++) fac(A[i]); //分解质因数
sort(P + 1, P + 1 + cnt);
int Max = 0, total = 0;
for(int i = 1 ; i <= cnt ; i ++) Max = max(Max, mp[P[i]]);
for(int i = 1 ; i <= cnt ; i ++) if(mp[P[i]] == Max) total ++;
cout << Max << endl;
s[1] = '1';
while(total -- ) Times();
int len = strlen(s + 1);
int x = s[len]; x --; s[len] = char(x); // 2 的任意正整数次方末尾只可能是 2,4,8,6 直接减没有问题
for(int i = 1 ; i <= len ; i ++) cout << s[i];
return 0;
}
习题:
洛谷 P4718 【模板】Pollard-Rho算法 (需要O(1) 快速乘 + 倍增优化判环)
P5071 [Ynoi2015] 此时此刻的光辉
随机化贪心
这东西真的就是玄学了......但是这个其实很简单。
算法用途：考试拿部分分
算法核心：随机化，贪心
引入:
思考下面的一个问题:
给定 \(n\) 个正整数，你现在要把 \(n\) 个正整数分成两组，定义某一组的价值为其组内所有数的和，你现在需要在两组选的数的个数差不超过 \(1\) 的前提下，最小化两组的价值差，输出最小的价值差。
\(n <= 100,1 <= num_i <= 450\)
然后一上来，神仙们肯定就开始秒切了这一题，不就是个 \(naive\) 的 \(dp\) 吗？
设置状态 \(dp[i][j][k]\) 表示.....balabala。 然后您就秒杀了这一题。
模拟退火
算法用途:抢劫一般的卷走大量部分分甚至AC
算法核心:随机化
模拟退火是什么？
模拟退火是一种随机化算法。当一个问题的方案数量极大（甚至是无穷的）而且不是一个单峰函数时，常使用模拟退火求解。 -----摘自 OiWiKi
模拟退火适用的问题通常是一些求最优解的问题（不是啥都能退火的）
引入
现在要求解一个最优解问题，可以把答案的函数抽象一下形成一个这样的函数：
然后对应的最优解假设就是函数的最低点。现在假设从 \(A\) 点开始找最优解。
(上面两张图片均摘自: FlashHu的模拟退火总结)
会在 \(B\) 点停下，因为它比两侧都要低，但是它并非是全局最优解，它只是局部最优解。
怎么解决这个只能找局部最优解的问题？（相当于对于爬山算法进行优化）
可以发现，爬山算法在实现的时候总是贪心的取临近的两个状态中的最优解，但是这样子就不会取接受更劣的解，这就导致很难到达全局最优解。所以现在，要想出一种办法来接受这个比较劣的解但是又不能只接受劣的解。
算法原理：
对于状态的选取加入概率的因素：如果新状态的解更优则修改答案，否则以一定概率接受新状态。
为了实现这一过程，科学家们决定模拟物理的退火工艺来实现一套算法。
一个处于很高温度的物体，现在要给它降温，使物体内能降到最低。
常规的思维是，越快越好，让它的温度迅速地降低。
然而，实际上，过快地降温使得物体来不及有序地收缩，难以形成结晶。而结晶态，才是物体真正内能降到最低的形态。
正确的做法，是徐徐降温，也就是退火，才能使得物体的每一个粒子都有足够的时间找到自己的最佳位置并紧密有序地排列。
开始温度高的时候，粒子活跃地运动并逐渐找到一个合适的状态。
在这过程中温度也会越降越低，温度低下来了，那么粒子也渐渐稳定下来，相较于以前不那么活跃了。这时候就可以慢慢形成最终稳定的结晶态了。
引用自 FlashHu 的博客。
因此得名，模拟退火。
算法思路
首先定义概念：
-
\(T\) —— 表示当前温度
-
\(\Delta T\) —— 温度变化率，温度变化率，每次温度等于上一次温度乘上 \(\Delta T\)， 一般 \(\Delta T\) 设置为 0.95 ~0.99 之间。
-
\(x\) 表示当前选择的点（也就是最开始提到的，现在处在的点）
-
\(f(x)\) 表示 \(x\) 对应的函数值，也就是目前的解。
-
\(x'\) 表示当前待选择的解
-
\(\Delta f\) 表示的是 \(f(x') - f(x)\) 也就是移动到待选择点的话答案会怎么变。
(下面假设是要求一个最优解是要尽可能小的)
那么假如 \(f(x') < f(x)\) , 那么说明新的状态比当前状态更加优秀，所以应该转移到 \(x'\),令 \(x = x'\)，准备开始下一次的退火。
那么假如 \(f(x') >= f(x)\), 那么取 \(e^\frac{-\Delta f}{T}\) 的概率随机判断取这个 \(x'\) 还是不取 (注意后面的 \(\frac{-\Delta f}{T}\) 是 \(e\) 的次数而不是将 \(e\) 与 \(\frac{-\Delta f}{T}\)乘起来)。至于为什么是这个概率，你得去问问科学家们是怎么得来的了，本蒟蒻实在不会证！
过程中记录出现的解中的最优解当作是全局最优解，所有同时得到最优解后，在最优解的附近继续随机几个点来判断是否更优。实现完这一套过程后，就认定得到的是全局最优解，然后存下来。
具体实现
因为是退火，所以需要设定一个初始温度 \(T\)(通常很大)，温度变化率 \(\Delta T\),终止的温度: \(T_d\)
其实上面的参数在于控制随机的次数，\(T\) 的作用则是调整每次随机取点的范围。
首先选一个点当作 \(x\) (也就是一开始在的点)，当然具体题目要具体调整，因为可能可以通过人工计算让一开始得到的这个起点就本身是比较优秀的点。
具体流程：
取 \(x'\) = \(x\) + \(T\) * \(Rand()\) (这里人为调整 \(Rand()\) 的范围为 \([1, -1]\) 的实数即可)
然后判断按照前面提到的，如果 \(x'\) 更优，毫无疑问的选择 \(x'\), 否则采用 \(exp\) 函数（包含于 #include < cmath>）求出 \(e^\frac{-\Delta f}{T}\) ，然后再随机的去判断要不要取。
值得注意的是，要在整个退火的过程中记录下出现过的，最优的解，这个才是最后得到的 \(x\), 而非退火完最后的 "\(x\)"
然后 \(T = T * \Delta T\)，直到 \(T < T_d\) 就终止算法。然后在最后得到的最优的 \(x\) 附近继续用上面的方法使用剩下的余温 \(T\) 来取 \(x'\), 看能不能有更优的答案即可。
最后就可以输出答案了。
例题：
POJ2420. 星星还是树
问题描述:
给定 \(n, (1 <= n <= 100)\) 个二维坐标系中的点的坐标，其中第 \(i\) 个点的坐标为 \((x_i,y_i)\) (\(0 <= x_i,y_i <= 10000\))。你现在要在平面中任意选取一个点，使得 \(\displaystyle \sum_{i = 1}^{i = n}{dis_i}\) 最小，问你最小距离和是多少。
solution
经典模拟退火入门题，听了上面的讲课应该不难解决这道题目了。先将答案坐标设置为:
\(nowx = \frac{\sum{x_i}}{n}, nowy = \frac{\sum{y_i}}{n}\)
然后就进行模拟退火，按照上面说的退火流程来搞就不难，每次计算答案的时候更新最优解，看计算结果是否比目前最优解更优，如果更优就更改，放在这里主要是想让同学们看看具体的代码实现。
const int MAXN = 105;
int n;
double ans,nowx,nowy,ansx,ansy;//nowx 和 nowy 表示是当前解，并不是最优解
double X[MAXN],Y[MAXN];
double sqr(double x) {return x * x;}
double Rand() {return (double)rand() / (double) RAND_MAX;}
double calc(int x, int y) {
double dis = 0;
for(int i = 1 ; i <= n ; i ++)
dis += sqrt(sqr(x - X[i]) + sqr(y - Y[i]));
if(dis < ans)ans = dis, ansx = x, ansy = y;//ansx记录的是目前最优的解的横坐标，ansy同理
return dis;
}
void SA() { //退火
double T = 10000, D = 0.988, Td = 1e-5; //设定初温，变化率以及终止温度
while(T >= Td) {
double newx = nowx + T * (Rand() * 2 - 1.0);//产生一个新的点对，判断跟原来点对的谁优谁劣
double newy = nowy + T * (Rand() * 2 - 1.0);
double Delta = calc(newx, newy) - calc(nowx, nowy);// f(x') - f(x)
if(exp(-Delta / T) > Rand()) nowx = newx, nowy = newy;
T *= D;
}
for(int i = 1 ; i <= 6666 ; i ++) { //退火完成后继续随机那么个几千次什么的。
double newx = ansx + T * (Rand() * 2 - 1.0);
double newy = ansy + T * (Rand() * 2 - 1.0);
calc(newx, newy);
}
return ;
}
一点 trick
while ((double)clock()/CLOCKS_PER_SEC<MAX_TIME) SA();
其中 MAX_TIME 设置为 0.7 ~ 0.8
在代码里面加入这一句话，然后就可以不用顾虑了超时，程序就会尽量多的运行退火算法而且不超时。
但是貌似好像只能在 windows 环境下使用？
然后就是分块模拟退火：
这里的分块并非是分成 \(\sqrt{n}\) 什么的，而是分成比较小的段，然后对于每一段跑一次模拟退火，取最优解，然后就作为答案。
习题
HAOI2006 均分数据 （这个玩意貌似可以不用退火，就乱随机就完事了）
其实所有的找最优解的题目都可以用模拟退火试一试。
参考博客：
- 自为风月马前卒 - \(MillerRabin\) 算法详解
- \(Pollard Rho\) 算法简介
- OI Wiki ---质因数分解
- FlashHu's模拟退火总结（模拟退火）
- OI Wiki ---模拟退火
- M-sea的模拟退火总结
扩展资料:
- \(Miller-Rabin\) 素性测试的错误率分析论文
- 快速乘的介绍
- 非平凡因子 : 假若有一个数 \(x\) (\(1 < x < n\)), 并且 \(x | n\) 则称 \(x\) 是 \(n\) 的非平凡因子。
