Description
Math is important!! Many students failed in 2+2’s mathematical test, so let's AC this problem to mourn for our lost youth..
Look this sample picture:
A ellipses in the plane and center in point O. the L,R lines will be vertical through the X-axis. The problem is calculating the blue intersection area. But calculating the intersection area is dull, so I have turn to you, a talent of programmer. Your task is tell me the result of calculations.(defined PI=3.14159265 , The area of an ellipse A=PI*a*b )
Input
Input may contain multiple test cases. The first line is a positive integer N, denoting the number of test cases below. One case One line. The line will consist of a pair of integers a and b, denoting the ellipse equation
, A pair of integers l and r, mean the L is (l, 0) and R is (r, 0). (-a <= l <= r <= a).
Output
For each case, output one line containing a float, the area of the intersection, accurate to three decimals after the decimal point.
Sample Input
2
2 1 -2 2
2 1 0 2
Sample Output
6.283
3.142
虽然说是Simpson积分……不过直接用普通的积分也是可行的
由于椭圆上下对称，故我们对上半部分积分即可
将 (frac{x^2}{a^2}+frac{y^2}{b^2}=1) 化简得 (y=bsqrt{1-frac{x^2}{a^2}}) ，将(a)提出来得 (y=frac{b}{a}sqrt{a^2-x^2})
我们考虑求 (intsqrt{a^2-x^2} dx)，令 (x=asin t)，则 (dx=acos t dt)
故原式化为 (intsqrt{a^2-a^2sin^2t} acos t dt=int a^2cos^2t dt=frac{a^2}{4}int (cos 2t+1) d(2t)=frac{a^2}{4}(sin 2t+2t)+C)
再化简得 (frac{asin t}{2}acos t+frac{a^2t}{2}+C)
将 (x=asin t) 代入得 (frac{x}{2}sqrt{a^2-x^2}+frac{a^2}{2}arcsin(frac{x}{a})+C)
最后原式需要计算(S=int_l^ry(x) dx)，我们将(y(x)=frac{b}{a}sqrt{a^2-x^2})代入，再根据上述推导过程推导即可
也可以利用椭圆的左右对称，直接从0积分到(l)或(r)，然后根据其位置加减即可
/*program from Wolfycz*/
#include<map>
#include<cmath>
#include<cstdio>
#include<vector>
#include<cstring>
#include<iostream>
#include<algorithm>
#define Fi first
#define Se second
#define ll_inf 1e18
#define MK make_pair
#define sqr(x) ((x)*(x))
#define pii pair<int,int>
#define int_inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline char gc(){
static char buf[1000000],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
template<typename T>inline T frd(T x){
int f=1; char ch=gc();
for (;ch<'0'||ch>'9';ch=gc()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=gc()) x=(x<<1)+(x<<3)+ch-'0';
return x*f;
}
template<typename T>inline T read(T x){
int f=1; char ch=getchar();
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';
return x*f;
}
inline void print(int x){
if (x<0) putchar('-'),x=-x;
if (x>9) print(x/10);
putchar(x%10+'0');
}
double F(int a,int x){return x/2.0*sqrt(sqr(a)-sqr(x))+sqr(a)/2.0*asin(1.0*x/a);}
int main(){
int T=read(0);
while (T--){
int a=read(0),b=read(0),l=read(0),r=read(0);
double resl=F(a,abs(l)),resr=F(a,abs(r)),Ans=0;
// printf("%lf %lf
",resl,resr);
if (l*r>0) Ans=abs(resl-resr);
else Ans=resl+resr;
printf("%.3lf
",2.0*b/a*Ans);
}
return 0;
}