前言

为什么我去学了模拟退火呢？

因为 CSP2021T1 和 NOIP2021T3 都能用模拟退火得到很高的分数，但是我却只打了最简单的暴力。失去了很多分数。

本着补全科技树的想法，我稍微学了一下模拟退火。

下文参考了 OI Wiki。

注：未完工。

前置知识——随机数

众所周知计算机的随机数是伪随机数。

那么我们首先需要让他变得随机一些。

下面依次介绍几个 C++ 库中生成随机数的方法。

（因为不是非常常用，我怕我忘记，所以有这部分内容。）

rand()

需要 <cstdlib> 库。

用法：直接调用 rand() 返回一个 $[0,\text{RAND\_MAX}]$​​ 的随机数，Linux 下这个上界是 $2^{31}-1$​​​​。

随机种子：用 srand() 修改，否则使用默认种子。

优点：好记，常见。

缺点：慢，不是非常随机。

std::mt19937

需要 <random> 库。

用法：这是个类，可以用 std::mt19937 myrand() 来声明一个。重载了括号运算符，实际使用时和 rand() 一样，随机数范围同 unsigned int 范围。

优点：快，更强一些。

缺点：不太好记。

还有个 std::mt19937_64，范围同 unsigned long long。

简介

假设有一个很莫名其妙的函数，要求它的极值，模拟退火算法可以帮助我们找到答案。（但不一定是函数）

首先，假设我们有一个局部最优解。

假设我们随机到了另一个状态，它更优，我们肯定更新最优解。

如果它更劣呢？我们不能把它忽略，因为这样就丧失了找到整体最优解的机会；但也不能直接就跳过去，这样局部最优解就找不到了。

那么我们就以一定概率去接受这个新状态。

定义当前温度 $T$，有一个已知状态，通过已知状态随机得到一个新状态，两者差为 $\Delta E$（$\Delta E \geq0$），那么发生状态转移（更新答案）的概率为： $$ P(\Delta E)= \begin{cases} 1&\text{新状态更优}\\ e^\frac{-\Delta E}{T}&\text{新状态更劣} \end{cases} $$ 解释一下：显然我们需要让差越大转移的概率越低，而且温度越低转移到劣解的概率越低。

实现

有三个参数 $T_0,d,T_{end}$。分别为初始温度（一个较大的数），降温系数（接近 $1$ 的数），终止温度（接近 $0$ 的正数）。

首先让当前温度 $T=T_0$，然后尝试转移一次，接着让 $T\leftarrow T\times d$，直到 $T<T_{end}$。得到最优解。

在找到最优解之后可以再用初始很低的温度再稍微跑一跑模拟退火，可能会找到更优解。

关于这三个参数，需要自己根据样例调整。（或者和暴力对拍）

小技巧——卡时

clock() 函数返回程序运行时间，于是可以这么写：

while ((double)clock() / CLOCKS_PER_SEC < 0.95)
    SA();

后面的 0.95 是自定义的一个数，要小于时限。

这样可以多执行几遍模拟退火，增加准确性。

NOIP2021 方差

题面不说了。

这个题正解不是模拟退火，不过可以用模拟退火得到 $60+$​ 分，而且非常简单，超值。

需要注意的是不能每次对差分数组随机一个排列，这样没意义，不满足模拟退火的下一个状态要从这一个状态转移来的要求。

所以直接随机选一个执行操作即可。

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <iostream>
#include <random>
using std::cout;
using std::endl;

const int maxn = 1e4 + 10;
typedef long long ll;

int a[maxn], n, b[maxn];
ll sumpow, sum;

ll calc()
{
    sumpow = 0, sum = 0;
    for (int i = 1; i <= n; i++)
        sumpow += b[i] * b[i], sum += b[i];
    return sumpow * n - sum * sum;
}

inline double rnd01() { return (double)rand() / RAND_MAX; }

inline void modify(int pos) { b[pos] = b[pos - 1] + b[pos + 1] - b[pos]; }

ll ans = 1e18;
std::mt19937 rng;

void SA()
{
    double T0 = 100000, d = 0.99997, Tend = 0.00001;
    for (int i = 1; i <= n; i++)
        b[i] = a[i];
    ll nowans = 1e18;
    while (T0 > Tend)
    {
        int pos = rng() % (n - 2) + 2;
        modify(pos);
        ll res = calc();
        ans = std::min(ans, res);
        if (res <= nowans)
            nowans = res;
        else if (exp((double)(nowans - res) / T0) > rnd01())
            nowans = res;
        else
            modify(pos);
        T0 *= d;
    }
}

int main()
{
    scanf("%d", &n);
    rng = std::mt19937(time(0));
    srand(time(0));
    for (int i = 1; i <= n; i++)
        scanf("%d", &a[i]);
    while ((double)clock() / CLOCKS_PER_SEC < 0.85)
        SA();
    cout << ans << endl;
}