UOJ156 【清华集训2015】恐怖的奴隶主题解

题目链接：#156. 【清华集训2015】恐怖的奴隶主

题意：

"A fight? Count me in!" 要打架了，算我一个。

"Everyone, get in here!" 所有人，都过来！

冷酷的酒客（Grim Patron）对战歌指挥官（Warsong Commander）的削弱感到很伤心，他打算用一道数学题来纪念战歌指挥官。

设 \(\left \{ u_i \right \}\) 为一数列

\[ \begin{equation} u_i = a + \frac{b}{u_{i-1}} + \frac{c}{u_{i-1}u_{i-2}} \quad (i \geq 2) \end{equation} \] 已知 \(u_0, u_1, a, b, c, n\)。保证 \(x^3 = a x^2 + b x + c\) 有 \(3\) 个不同的正整数解。

求 \(u_n\)。

输入格式：

输入仅一行，包含 \(6\) 个整数 \(u_0, u_1, a, b, c, n\)。

输出格式：

输出仅一行，包含 \(1\) 个整数，保留至少 \(8\) 位至多 \(15\) 位小数。如果你的答案和我们的答案差别不超过 \(10^{-6}\)，则认为正确。考虑到浮点运算本身的误差，当你的答案与真实答案差别不超过 \(10^{-6} - 10^{-8}\) 时，才能保证正确。

数据范围：

对于 \(100\%\) 的数据，\(c \leq 100000\)， \(\lvert u_0 \rvert \leq 1000, \lvert u_1 \rvert \leq 1000\) ，\(0 \leq n \leq 10^9\) 。

很巧妙的一道题，不知道为什么差评那么多（~~谁想得到这种东西啊喂~~

首先从 \(x^3 = a x^2 + b x + c\) 有 \(3\) 个不同的正整数解出发

设 \(x\) 是这个方程的某个解

一般这种题目的数列都是收敛数列，即有极限的数列。

注：数列的极限定义如下：（来自《数学分析》B.A.卓里奇）

数 \(A \in \mathbb{R}\) 称为数列 \(\{x_n\}\) 的极限，如果对于任何 \(\epsilon > 0\) ，都存在序号 \(N\) ，使得对于一切 \(n>N\) ，都有 \(|x_n - A| < \epsilon\)

形式化地， \[ \left(\lim\limits_{n \to \infty} x_n = A\right):= \forall \epsilon > 0 ,\exists N \in \mathbb{N},\forall n > N,(|x_n - A| < \epsilon) \]

于是，我们尝试代入 \(u_i = u_{i-1}=x\)

会发现 \(u_{i+1} = x\) ，~~嗨害嘿~~，原数列落入了不动点

因此 \(\lim\limits_{n \to \infty}u_n = x\) ，现在关键是怎么找到这个 \(x\) （因为有三个根嘛）

对于「 \(n \le 60\) 左右」的情况，可以暴力计算¹

当 \(n\) 特别大的时候，我们只要将得到的数舍入

然后判断那个数是不是根即可，注意要判断其收敛趋势

时间复杂度 \(O(1)\)

代码：

#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cstdarg>
#include <cmath>
#include <iomanip>
#include <random>
using namespace std;
#define int long long
#define double long double
#define INF 0x3f3f3f3f3f3f3f3f
#define N 66

const double eps=1e-10;
int n,a,b,c;
double u[N],rd;
bool ok(int x){return ((x - a) * x - b) * x == c;}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    cout << fixed << setprecision(15);
    cin >> u[0] >> u[1] >> a >> b >> c >> n;
    if(n > N)
    {
        for(int i=2; i<N; i++)
        {
            u[i] = a + (b + c / u[i-2]) / u[i-1];
            rd=roundl(u[i]); // 四舍五入
            if(fabs(u[i] - rd) < 0.01 && fabs(u[i-1] - rd) < 0.01)
                if(ok(rd)) return cout << rd << '\n',0;
        }
    }else
    {
        for(int i=2; i<=n; i++) u[i] = a + (b + c / u[i-2]) / u[i-1];
        cout << u[n] << '\n';
    }
    return 0;
}