Keyboard Free （计算几何+自适应辛普森积分/等分圆）2020牛客多校第二场

原题题面

Given three concentric circles whose radiuses are $r_1$,$r_2$,$r_3$respectively, and $A,B,C$ are the moving points on the given three circles respectively. Determine the expected area of $△ABC$ .

输入描述

The first line contains one integer $T(1\le T\le 1000)$, denoting the number of test cases.
For each test case:
One line containing three integers $r_1,r_2,r_3(1\le r1,r2,r3\le 100)$, denoting the radiuses of three given concentric circles.

输出描述

Print $T$ lines each containing one real number with one decimal places after the decimal point, denoting the answer to curresponding test case.
It’s guaranteed that the second decimal place after the decimal point is neither 4 nor 5.

输入数据

2
1 1 1
2 3 5

输出数据

0.5
5.5

题面解析

一道好难好难的几何题。
不妨设$r_1\le r_2\le r_3$。
首先容易考虑到的是让$A$变为定点，然后设点$B,C$在列方程求解，然后就会得到如下的二重积分：
里面有绝对值，而且无法分类讨论（r1,r2,r3是未知的常数），故放弃该方法，采用其他策略。
可以采用如下办法。
我们让$A$作为定点，$B,C$为动点，算出C到AB的期望距离后乘以AB的长度除以2。

对于该图，我们可以旋转坐标系（你可以试着旋转你的头来治好颈椎病 ），让$P_{1}A$垂直于$x$轴。

我们设$\angle O{P}_{1}G=\alpha$，设$C(r_{3}cos\beta ,r_{3}sin\beta )$,那$C$到$AB$的直线距离的期望是
$\frac{1}{2\pi }\ast {\int }_0^{2\pi }|sin\alpha -cos\beta |\mathrm{d}\beta$.

=$\frac{1}{2\pi }\ast 2{\int }_0^{\pi }|sin\alpha -cos\beta |\mathrm{d}\beta$.（这样就可以很好地处理绝对值的问题了）

=$\frac{1}{2\pi }\ast 2({\int }_{\frac{\pi }{2}-\alpha }^{\pi }(sin\alpha -cos\beta )\mathrm{d}\beta +{\int }_0^{\frac{\pi }{2}-\alpha }(cos\alpha -sin\beta )\mathrm{d}\beta )$.

=$\frac{1}{2\pi }\ast 4r_3(\alpha sin\alpha +cos\alpha )$

即$C$到$AB$的直线距离的期望只和$\alpha$有关。
有了这个结论之后，我们可以设$A(r_1,0),B(r_{2}cos\beta ,r_{2}sin\beta )$，把$\alpha$用$\beta$表示。

如图所示，我们做$BD\perp x$轴，则$△OAG\backsim △BAD$(AA,即两个角相等),因为$\angle ABD$可以计算，所以$OG$就可以计算，进而推出$\alpha$，可以得到
$$sin\alpha =\frac{r_1}{r_3}\ast \sqrt{\frac{(r_{2}sin\beta {)}^2}{(r_{2}sin\beta {)}^2+(r{s}_{c}os\beta -r_1{)}^2}}$$
所以最后的答案是
$\frac{1}{2\ast 2\pi \ast 2\pi }{\int }_0^{2\pi }(4r_3(\alpha sin\alpha +cos\alpha )\sqrt{(r_{2}sin\beta {)}^2+(r_{2}cos\beta -r1{)}^2}\mathrm{d}\beta$
当然，这个东西的精确值显然超出笔者范围了，于是笔者采用了自适应辛普森积分来做（其实就是把圆正等分成上千个离散点，分别代入求结果）
需要注意的是，如果采用自适应辛普森积分的话，请注意精确度的大小控制，否则会超时。

AC代码（487ms）

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-4;//精确度要注意不能太小，否则超时
const double PI=acos(-1.0);
long long r1, r2, r3;
double f(double a)//积分内的式子
{
    if(a==0)
        return 0;//特殊处理，否则会RE
//    printf("a=PI/%f\n",PI/a);
    double ans1=sqrt(r2*sin(a)*r2*sin(a)+(r2*cos(a)-r1)*(r2*cos(a)-r1));
//    printf("ans1=%f\n",ans1);
    double angle1=atan((r2*sin(a))/(r2*cos(a)-r1));
    double b=asin(r1*sin(angle1)/r3);
//    printf("b=%f\n");
    double ans2=4*r3*(cos(b)+b*sin(b));
//    printf("ans2=%f\n", ans2);
    return ans1*ans2;
}
double SimpleSimpson(double a, double b)//辛普森积分
{
    return (b-a)/6.0*(f(a)+f(b)+4*f((a+b)/2.0));
}
double Simpson(double l, double r, double ans)//自适应辛普森积分
{
    double mid=(l+r)/2;
//    printf("mid:%f\n",mid);
    double Simpson_l=SimpleSimpson(l, mid);
    double Simpson_r=SimpleSimpson(mid, r);
//    printf("l:%f r:%f ans:%f\n",Simpson_l,Simpson_r,ans);
    if (fabs(Simpson_l+Simpson_r-ans)<=eps)
        return ans;
    return Simpson(l, mid, Simpson_l)+Simpson(mid, r, Simpson_r);
}
void solve1()
{
    int t;
    scanf("%d", &t);
    while(t--)
    {
        scanf("%lld%lld%lld", &r1, &r2, &r3);
        if (r1>r2)
            swap(r1, r2);
        if (r2>r3)
            swap(r2, r3);
        if (r1>r3)
            swap(r1, r3);
//        printf("r1:%lld r2:%lld r3:%lld\n",r1,r2,r3);
//        printf("%f %f %f\n",f(0),f(PI),f(PI/2.0));
//        printf("%f\n",SimpleSimpson(0, PI));
        printf("%.1f\n", Simpson(0, 2*PI, SimpleSimpson(0, 2*PI))/(2*PI*2*PI*2));
    }
}
int main()
{
//    ios_base::sync_with_stdio(false);
//    cin.tie(0);
//    cout.tie(0);
#ifdef ACM_LOCAL
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
    long long test_index_for_debug=1;
    char acm_local_for_debug;
    while(cin>>acm_local_for_debug)
    {
        cin.putback(acm_local_for_debug);
        if (test_index_for_debug>100)
        {
            throw runtime_error("Check the stdin!!!");
        }
        auto start_clock_for_debug=clock();
        solve1();
        auto end_clock_for_debug=clock();
        cout<<"\nTest "<<test_index_for_debug<<" successful"<<endl;
        cerr<<"Test "<<test_index_for_debug++<<" Run Time: "
            <<double(end_clock_for_debug-start_clock_for_debug)/CLOCKS_PER_SEC<<"s"<<endl;
        cout<<"--------------------------------------------------"<<endl;
    }
#else
    solve1();
#endif
    return 0;
}

后记

感谢群聚XLor的指导，耽误了您太多时间真不好意思…
DrGilbert 2020.7.14