HDU 4851 Wow! Such Precision!

题目：http://acm.hdu.edu.cn/showproblem.php?pid=4851

题意：
给定一个非负整数$n$，请你构造一个$n$次多项式$f(x)$，使得$\int_{0}^{1}{(f(x)-e^x)^2dx}$取得最小值。给出$f(x)$的$n+1$个项的系数，要求精度误差不超过$10^{-50}$。
$0\le n\le 10$。

题解：
这个$n$实在是太小了，给人一种可以打表的感觉，于是先来打打表试试！
对于$n=0$，设$f(x)=c$，则

$$\int_{0}^{1}{(f(x)-e^x)^2dx}\\=\int_{0}^{1}{(e^{2x}-2ce^x+c^2)dx}\\=(\frac{e^{2x}}{2}-2ce^x+c^2x)|_0^1\\=\frac{e^2-1}{2}-2c(e-1)+c^2$$
这个函数的极小值点满足导数为 $0$，也就是$-2(e-1)+2c=0$，即 $c=e-1$。可以注意到常数项是不影响答案的。
对于 $n=1$，一通强算之后也能得到 $f(x)=(18-6e)x+(4e-10)$，然而对于更高的情况就很难强算啦，所以需要一个更形式化的做法。
令 $f(x)=\sum_{i=0}^{n}{a_i\cdot x^i}$，那么有 $$\int_{0}^{1}{(f(x)-e^x)^2dx}\\=\int_{0}^{1}{(f^2(x)-2f(x)e^x+e^{2x})dx}\\=\sum_{i=0}^{n}\sum_{j=0}^{n}{a_ia_j\int_{0}^{1}{x^{i+j}dx}}-\sum_{i=0}^{n}{2a_i\int_{0}^{1}{x^ie^xdx}}+\int_{0}^{1}{e^{2x}dx}\\=\sum_{i=0}^{n}\sum_{j=0}^{n}{\frac{a_ia_j}{i+j+1}}-\sum_{i=0}^{n}{2a_i\int_{0}^{1}{x^ie^xdx}}+\frac{e^2-1}{2}$$
注意到这里有一个 $\int_{0}^{1}{x^ie^xdx}$没有化简，考虑分部积分的形式 $$\int_{0}^{1}{x^ie^xdx}=(x^ie^x)|_0^1-i\int_{0}^{1}{x^{i-1}e^xdx}=e-i\int_{0}^{1}{x^{i-1}e^xdx}$$
而 $\int_{0}^{1}{e^xdx}=e-1$，所以可以归纳证明 $\int_{0}^{1}{x^ie^xdx}$一定可以被表示成 $b_ie+c_i$的形式（因为之间的关系是数乘），则有 $$\left\{\begin{matrix}b_0=1\\b_i=1-i\cdot b_{i-1}\end{matrix}\right.\left\{\begin{matrix}c_0=-1\\c_i=-i\cdot c_{i-1}\end{matrix}\right.$$
因此积分后的函数就是 $$\sum_{i=0}^{n}\sum_{j=0}^{n}{\frac{a_ia_j}{i+j+1}}-\sum_{i=0}^{n}{2a_i(b_ie+c_i)}+\frac{e^2-1}{2}$$
现在要求这个 $(n+1)$元函数的最小值，可以使用拉格朗日乘数法。
令 $L(a_0,a_1,\cdots,a_n)=\sum_{i=0}^{n}\sum_{j=0}^{n}{\frac{a_ia_j}{i+j+1}}-\sum_{i=0}^{n}{2a_i(b_ie+c_i)}+\frac{e^2-1}{2}$，则有偏导 $\frac{\partial L}{\partial a_i}=2\sum_{j=0}^{n}{\frac{a_j}{i+j+1}}-(b_ie+c_i)$。
极小值的点一定是稳定点，且二阶偏导是正定的，一阶偏导等于 $0$，于是就有了$(n+1)$个形如 $\sum_{j=0}^{n}{\frac{a_j}{i+j+1}}=b_ie+c_i$的方程。
由于方程未知数的系数是有理数，所以 $a_j$也一定能写成 $k_je+t_j$的形式，于是分别对于 $b_i$和 $c_i$进行一次高斯消元，就可以得到对应的 $k_j$和 $t_j$。
然而此题要求的精度十分离谱，所以需要更进一步的分析。可以归纳证明 $k_j$和 $t_j$都是整数，但是数字可能大于 $10^9$，但不会太大。
于是找两个大素数做模意义下的高斯消元，再用中国剩余定理计算出其实际值即可（需要实际值小于模数）。
不过我找的两个大素数是 $1000000007$和 $1000000009$，而在 $n=10$的时候系数已经超过了模数（因为算出的 $f(x)$系数居然出现了负数），所以我选择了第三个素数 $998244353$，得到对应的同余方程组之后，用python做中国剩余定理算出了实际值，大约有 $4\cdot 10^{20}$左右，刚好超出 $64$位整型数的表示范围。
粗略估计为了保证精度，需要 $e$在小数点后$70$位都是精确的，所以泰勒展开求出一个保留小数点后 $70$位的值之后，将所有数字乘以 $10^{70}$转化为整数计算，得出每个系数的值，打入表中输出即可。
分析复杂度的话，最大的复杂度部分是 $O(n^3logn)$的高斯消元，最后的高精度计算是常数的主要部分。
我是打成字符串后交的C++代码，但是可以使用Java的BigDecimal进行高精度浮点数的运算，会方便许多。

代码：
打表的代码感觉还是不要给了，等研究出Java的写法再补。