发现当把 \(\rm scale×cos\theta，scale×sin\theta,dx,dy\) 当作变量时只有四个，两个方程就行。

当 \(\rm n\le 500\) 时，可以选取两组进行高斯消元，解出答案后回带。

但当 \(n\) 极大时，采用随机化的做法，每次随机选取两个，这样每次选取不正确的概率为 \(\frac{3}{4}\)，选取 50 次后基本就会出答案了。

记得判断 \(\rm sin\) 的正负

Code

#include<bits/stdc++.h>
#define ri register signed
#define p(i) ++i
namespace IO{
    char buf[1<<21],*p1=buf,*p2=buf;
    #define gc() p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?(-1):*p1++
    struct nanfeng_stream{
        template<typename T>inline nanfeng_stream &operator>>(T &x) {
            ri f=0;x=0;register char ch=gc();
            while(!isdigit(ch)) {f|=ch=='-';ch=gc();}
            while(isdigit(ch)) {x=(x<<1)+(x<<3)+(ch^48);ch=gc();}
            return x=f?-x:x,*this;
        }
    }cin;
}
using IO::cin;
namespace nanfeng{
    #define FI FILE *IN
    #define FO FILE *OUT
    template<typename T>inline T cmax(T x,T y) {return x>y?x:y;}
    template<typename T>inline T cmin(T x,T y) {return x>y?y:x;}
    typedef double db;
    static const int N=1e5+7;
    static const db eps=1e-6;
    struct node{db x1,y1,x2,y2;}pnt[N];
    int p[N],n,hl;
    db mt[51][51],x1,x2,x3,x4;
    inline void Guass() {
        for (ri i(1);i<=4;p(i)) {
            ri k=i;
            for (ri j(i+1);j<=4;p(j)) if (fabs(mt[j][i])>fabs(mt[k][i])) k=j;
            if (k!=i) std::swap(mt[k],mt[i]);
            for (ri j(1);j<=4;p(j)) {
                if (i==j) continue;
                db tmp=mt[j][i]/mt[i][i];
                for (ri l(1);l<=5;p(l)) mt[j][l]-=tmp*mt[i][l];
            }
        }
    }
    inline bool judge() {
        ri cnt(0);
        x1=mt[1][5]/mt[1][1];
        x2=mt[2][5]/mt[2][2];
        x3=mt[3][5]/mt[3][3];
        x4=mt[4][5]/mt[4][4];
        for (ri i(1);i<=n;p(i))
            if (fabs(pnt[i].x1*x1-pnt[i].y1*x2+x3-pnt[i].x2)<=eps&&fabs(pnt[i].x1*x2+pnt[i].y1*x1+x4-pnt[i].y2)<=eps) p(cnt);
        return cnt>=hl;
    }
    inline int solve() {
        register db sc=sqrt(x1*x1+x2*x2);
        if (x2/sc<=eps) printf("%.10lf\n",acos(-1)*2.0-acos(x1/sc));
        else printf("%.10lf\n",acos(x1/sc));
        printf("%.10lf\n%.10lf %.10lf\n",sc,x3,x4);
        return 0;
    }
    inline int main() {
        //FI=freopen("nanfeng.in","r",stdin);
        //FO=freopen("nanfeng.out","w",stdout);
        srand(unsigned(time(0)));
        std::cin >> n;
        hl=(n>>1)+(n&1);
        for (ri i(1);i<=n;p(i)) {
            scanf("%lf%lf%lf%lf",&pnt[i].x1,&pnt[i].y1,&pnt[i].x2,&pnt[i].y2);
            p[i]=i;
        }
        if (n>500) {
            std::random_shuffle(p+1,p+n+1);
            for (ri i(1);i<=n;p(i)) {
                ri cur1(p[i]),cur2(p[i+1]);
                mt[1][1]=pnt[cur1].x1,mt[1][2]=-pnt[cur1].y1,mt[1][3]=1.0,mt[1][4]=0.0,mt[1][5]=pnt[cur1].x2;
                mt[2][1]=pnt[cur1].y1,mt[2][2]=pnt[cur1].x1,mt[2][3]=0.0,mt[2][4]=1.0,mt[2][5]=pnt[cur1].y2;
                mt[3][1]=pnt[cur2].x1,mt[3][2]=-pnt[cur2].y1,mt[3][3]=1.0,mt[3][4]=0.0,mt[3][5]=pnt[cur2].x2;
                mt[4][1]=pnt[cur2].y1,mt[4][2]=pnt[cur2].x1,mt[4][3]=0.0,mt[4][4]=1.0,mt[4][5]=pnt[cur2].y2;
                Guass();
                if (judge()) return solve();
            }
        } else {
            for (ri i(1);i<=n;p(i))
                for (ri j(i+1);j<=n;p(j)) {
                    ri cur1(i),cur2(j);
                    mt[1][1]=pnt[cur1].x1,mt[1][2]=-pnt[cur1].y1,mt[1][3]=1.0,mt[1][4]=0.0,mt[1][5]=pnt[cur1].x2;
                    mt[2][1]=pnt[cur1].y1,mt[2][2]=pnt[cur1].x1,mt[2][3]=0.0,mt[2][4]=1.0,mt[2][5]=pnt[cur1].y2;
                    mt[3][1]=pnt[cur2].x1,mt[3][2]=-pnt[cur2].y1,mt[3][3]=1.0,mt[3][4]=0.0,mt[3][5]=pnt[cur2].x2;
                    mt[4][1]=pnt[cur2].y1,mt[4][2]=pnt[cur2].x1,mt[4][3]=0.0,mt[4][4]=1.0,mt[4][5]=pnt[cur2].y2;
                    Guass();
                    if (judge()) return solve();
                }
        }
        return 0;
    }
}
int main() {return nanfeng::main();}