题解 P5916 【[FJOI2014]病毒防护带】

我们先将要求的东西形式化。 $$ ans = \min_{A,B\in R} \frac{\pi}{A^2+1}\sum_{i=1}^n(Ax_i+B-y_i)^2 $$ 我们发现式子的形式很像最小二乘法最小化的式子，即 $$ \min_{A,B\in R} \sum_{i=1}^n(Ax_i+B-y_i)^2 $$ 如果考虑通过旋转坐标系的方式将 $A$ 转成 $0$ ，那么这两个式子就等价了。 我们就可以知道此时，我们要求的直线就是最小二乘法求出来的直线。 我们可以知道最小二乘法求的直线肯定经过 $(\overline x, \overline y)$ 。又因为旋转坐标系是线性操作，所以我们要求的直线 **必定经过 $(\overline x,\overline y)$** 。 那么就有 $A\overline x + B = \overline y$ ，可以得到 $B=\overline y - A\overline x$ 。这样 $ans$ 就可以表示成关于 $A$ 的函数。我们只需要对这个函数求最值就可以了。 那么转化一下式子: $\displaystyle ans = \min_{A,B\in R} \frac{\pi}{A^2+1}\sum_{i=1}^n(Ax_i+B-y_i)^2$ $\displaystyle = \min_{A,B\in R} \frac{\pi}{A^2+1}\sum_{i=1}^n[Ax_i+(\overline y - A\overline x)-y_i]^2$ $\displaystyle = \min_{A,B\in R} \frac{\pi}{A^2+1}\sum_{i=1}^n[A(x_i-\overline x)+(\overline y-y_i)]^2$ $\displaystyle = \min_{A,B\in R} \frac{\pi}{A^2+1}\sum_{i=1}^nA^2(x_i-\overline x)^2+2A(x_i-\overline x)(\overline y - y_i)+(\overline y-y_i)^2$ 为了方便，我们设 $u =\sum_{i=1}^n(x_i-\overline x)^2$ ， $v = \sum_{i=1}^n2(x_i-\overline x)(\overline y - y_i)$ ， $z =\sum_{i=1}^n(\overline y - y_i)^2$ 。 这样我们需要的就是求函数 $ans(A)=\frac{uA^2+vA+z}{A^2+1}$ 的最小值。可以求导求极值。 $ans'(A)=\frac{(2uA+v)(A^2+1)-2A(uA^2+vA+z)}{(A^2+1)^2}$ $=ans'(A)=\frac{-vA^2+2(u-z)A+v}{(A^2+1)^2}$ 可得到两个极值点 $A=\frac{(u-z)\pm\sqrt{(u-z)^2+v^2}}{v}$ 。 我们可以对两个 $A$ 分别检验答案，然后取较小值。至于 $A=\pm \infty$ 为什么不行，这个把函数图像画一画就出来了。 注意当 $v=0$ 的时候极值点确实是 $A=+\infty$ 。不过数据随机保证了这种情况基本不存在。 顺便一提，就结果来看算出来的所有 $v$ 全都是负数，极值点取 $A=\frac{(u-z)+\sqrt{(u-z)^2+v^2}}{v}$ 。可能是因为数据随机吧。。。 顺便再一提，BZOJ上面的这道题不知道是出了什么鬼问题，连 hzw 的题解交上去都过不了，真的很迷。。。 是不是很复杂? 没事，得到 $B=\overline y -\overline x A$ 之后你就可以对 $A$ 三分了。 或者你就算没得到也可以直接对 $A$ 和 $B$ 同时三分。复杂度并不会改变很多，只是由 $O(\sum n)$ 变成了 $O(\sum n + \log ^2 n)$ 。 代码: ```cpp #include<bits/stdc++.h> #define Rep(i,a,b) for(register int i=(a);i<=(b);++i) #define Repe(i,a,b) for(register int i=(a);i>=(b);--i) #define rep(i,a,b) for(register int i=(a);i<(b);++i) #define pb push_back #define mp make_pair #define mx(a,b) (a>b?a:b) #define mn(a,b) (a<b?a:b) typedef unsigned long long uint64; typedef unsigned int uint32; typedef long long ll; using namespace std; namespace IO { const uint32 Buffsize=1<<15,Output=1<<24; static char Ch[Buffsize],*S=Ch,*T=Ch; inline char getc() { return((S==T)&&(T=(S=Ch)+fread(Ch,1,Buffsize,stdin),S==T)?0:*S++); } static char Out[Output],*nowps=Out; inline void flush(){fwrite(Out,1,nowps-Out,stdout);nowps=Out;} template<typename T>inline void read(T&x) { x=0;static char ch;T f=1; for(ch=getc();!isdigit(ch);ch=getc())if(ch=='-')f=-1; for(;isdigit(ch);ch=getc())x=x*10+(ch^48); x*=f; } template<typename T>inline void write(T x,char ch='\n') { if(!x)*nowps++='0'; if(x<0)*nowps++='-',x=-x; static uint32 sta[111],tp; for(tp=0;x;x/=10)sta[++tp]=x%10; for(;tp;*nowps++=sta[tp--]^48); *nowps++=ch; if(nowps-Out>=1<<23)flush(); } inline void getstr(char*q) { register char ch; for(ch=getc();!isgraph(ch);ch=getc()); for(;isgraph(ch);ch=getc())*q++=ch; *q='\0'; } inline void getwd(char&x){for(x=getc();!isupper(x);x=getc());} } using namespace IO; void file() { #ifndef ONLINE_JUDGE FILE*DSD=freopen("a.in","r",stdin); FILE*CSC=freopen("a.out","w",stdout); #endif } inline void Chkmin(int&u,int v){u>v?u=v:0;} inline void Chkmax(int&u,int v){u<v?u=v:0;} inline void Chkmax(double&u,double v){u<v?u=v:0;} inline void Chkmax(ll&u,ll v){u<v?u=v:0;} inline void Chkmin(ll&u,ll v){u>v?u=v:0;} const int MAXN = 1e7 + 7; static int n; static int X[MAXN], Y[MAXN]; static long double xx, yy; static int _; inline void init() { read(n), read(X[1]), read(Y[1]); assert((ll)_ * n <=2e7); xx = X[1], yy = Y[1]; static int a, b, c; read(a), read(b), read(c); Rep(i, 2, n) { X[i] = (((a * X[i - 1]) + b) * X[i - 1] + c) % 107; Y[i] = (((a * Y[i - 1]) + b) * Y[i - 1] + c) % 107; xx += X[i], yy += Y[i]; } xx /= n, yy /= n; } static long double ans, ans1; const long double pi = acos(-1); inline long double sq(long double x) {return x * x;} inline void solve(int zz) { long double u = 0, v = 0, z = 0, A, B; Rep(i, 1, n) { u += sq(X[i] - xx), v += 2 * (X[i] - xx) * (yy - Y[i]); z += sq(yy - Y[i]); } if(fabs(v) < 1e-10) { Rep(i, 1, n) ans += sq(X[i] - xx); printf("%.8Lf", ans * pi / n); exit(0); } long double delta = sqrt(4 * sq(u - z) + 4 * sq(v)); A = (-2 * (z - u) - delta) / 2 / v; B = yy - A * xx; ans = 0; long double sqq = A * A + 1; Rep(i, 1, n) ans += sq(A * X[i] + B - Y[i]); ans = ans * pi / sqq / n; A = (-2 * (z - u) + delta) / 2 / v; B = yy - A * xx; ans1 = 0; sqq = A * A + 1; Rep(i, 1, n) ans1 += sq(A * X[i] + B - Y[i]); ans1 = ans1 * pi / sqq / n; if(ans1 < ans) ans = ans1; printf("Case %d: %.5Lf\n", zz, ans); } int main() { file(); read(_); Rep(i, 1, _) init(), solve(i); return 0; } ```