题目大意
给一个简单多边形和一条抛物线,问抛物线在简单多边形内部的长度和是多少。
Tag
积分·辛普森
计算几何·点,线
解析几何
解法
给一个例子——多边形 ABCDEF和抛物线
我们算出抛物线和多边形的焦点 GHIJ。
那么就能很简单确定我们需要的积分区间 G->H 和 I -> J
于是我们要做的事情就简单了——
1.算出来每条线段与抛物线交点的 X 坐标
2.区分哪些是进入多边形的X哪些是离开多边形的X
3.对于进入-> 离开这段区间,计算抛物线长度,累加即可。
我们把这三步骤分成两步——1,2 算是计算几何,3算是积分。
计算几何部分
分以下几种情况——
1.边不与多边形相交 或者 交点不在线段上(也不在端点上) 或者 边与抛物线相切。
做法——直接忽略。
2. 边平行于Y轴。 直接计算交点
3. 边与抛物线交于 >= 1 个顶点,解方程。
有三种特殊情况——
对于上面两种情况,交点应该被忽略。
对于第三种情况,交点不应该被忽略!!!因为抛物线的确由此进入了多边形。
对于这种判断的方法,我的方法有些蠢——取交点左右线段上各 0.5 左右长度,看它们俩是不是一个在抛物线上方一个在抛物线下方(code中check函数)。如果是那么就属于第三种——跨立 状态。
对于要不要加上这个点,很难判断。。。那我们干脆这样,碰到交点就加入,对于前两种情况再加入一个一样的点使得把一个错误的点去掉。要特殊判断的是如果当前交点和下一个或者上一个点重合,那么我们需要 ignore 它。。
比如前两种 X 序列可能是 1 , 2 , 3 , 3 , 4 , 5.。。这样就3 -> 3 直接可以忽略了
对于第三种 X 序列可能是 1 , 2 , 3 , 3 , 3 , 4.。这样 3-》4 就能被计算上了。。
判断的方法是——按照一个方向,比如对于线段 p[i] - p[i + 1] 我每次都判断 p[i] 看是否属于 “跨立点”。
那么计算几何部分就没有压力了。。
积分部分
首先是 曲线积分长度 = sqrt(1 + (f'(x)) ^ 2) , 即根号下 一 加上 导数的平方。由于我高数几乎属于不及格的状态,我就直接套的 辛普森积分。
要特殊判断的是给定的 [l , r] 区间
我们这么考虑——
枚举每个区间 h[i] - h[i + 1] , 我们积分 L = fmax(l , h[i]) , R = fmin(r , h[i + 1]); L-> R 这段区间,如果 L >= R 那么说明积分结束。。
于是这个题目就这样解决了。。细节很多。
Code
namespace Geo{
#define typec double
const typec eps=1e-8;
int dblcmp(double d){
if (fabs(d)<eps)return 0;
return d>eps?1:-1;
}
int sgn(double a) {return a<-eps?-1:a>eps;}
inline double sqr(double x){return x*x;}
struct Point2D{
typec x,y;
Point2D(){}
Point2D(typec _x,typec _y):x(_x),y(_y){};
void input(){
scanf("%lf%lf",&x,&y);
}
void output(){
printf("%.2f %.2f\n",x,y);
}
bool operator==(Point2D a)const{
return dblcmp(a.x-x)==0&&dblcmp(a.y-y)==0;
}
bool operator<(Point2D a)const{
return dblcmp(a.x-x)==0?dblcmp(y-a.y)<0:x<a.x;
}
typec len(){
return hypot(x,y);
}
typec len2(){
return x*x+y*y;
}
Point2D operator + (const Point2D &A) const{
return Point2D(x + A.x , y + A.y);
}
Point2D operator - (const Point2D &A) const{
return Point2D(x - A.x , y - A.y);
}
Point2D operator * (const typec _x) const{
return Point2D(x * _x , y * _x);
}
typec operator * (const Point2D &A) const{
return x * A.x + y * A.y;
}
typec operator ^ (const Point2D &A) const{
return x * A.y - y * A.x;
}
Point2D operator / (const typec _p) const{
return Point2D(x / _p , y / _p);
}
typec distance(Point2D p){
return hypot(x-p.x,y-p.y);
}
Point2D add(Point2D p){
return Point2D(x+p.x,y+p.y);
}
Point2D sub(Point2D p){
return Point2D(x-p.x,y-p.y);
}
Point2D mul(typec b){
return Point2D(x*b,y*b);
}
Point2D div(typec b){
return Point2D(x/b,y/b);
}
typec dot(Point2D p){
return x*p.x+y*p.y;
}
typec det(Point2D p){
return x*p.y-y*p.x;
}
typec rad(Point2D a,Point2D b){
Point2D p=*this;
return fabs(atan2(fabs(a.sub(p).det(b.sub(p))),a.sub(p).dot(b.sub(p))));
}
Point2D trunc(typec r){
typec l=len();
if (!dblcmp(l))return *this;
r/=l;
return Point2D(x*r,y*r);
}
Point2D rotleft(){
return Point2D(-y,x);
}
Point2D rotright(){
return Point2D(y,-x);
}
Point2D rotate(Point2D p,typec angle)//èÆμãpÄæê±ÕëDy×aangle½Ç¶è
{
Point2D v=this->sub(p);
typec c=cos(angle),s=sin(angle);
return Point2D(p.x+v.x*c-v.y*s,p.y+v.x*s+v.y*c);
}
};
typec cross(Point2D a,Point2D b,Point2D c){
return (b.sub(a)).det(c.sub(a));
}
}using namespace Geo;
/* .................................................................................................................................. */
double a , b , c , l , r;
DB f0(DB x){ // Input the function
return sqrt(1. + sqr(2. * a * x + b));
}
/* helper function of adaptive_simpsons */
DB adaptive_simpsons_aux(DB (*f)(DB), DB a, DB b,
DB eps, DB s, DB fa, DB fb, DB fc, int depth) {
DB c = (a + b) / 2, h = b - a;
DB d = (a + c) / 2, e = (c + b) / 2;
DB fd = f(d), fe = f(e);
DB sl = (fa + 4 * fd + fc) * h / 12;
DB sr = (fc + 4 * fe + fb) * h / 12;
DB s2 = sl + sr;
if (depth <= 0 || fabs(s2 - s) <= 15 * eps)
return s2 + (s2 - s) / 15;
return adaptive_simpsons_aux(f, a, c, eps / 2, sl, fa, fc, fd, depth - 1) +
adaptive_simpsons_aux(f, c, b, eps / 2, sr, fc, fb, fe, depth - 1);
}
/* Adaptive Simpson's Rule, integral of f at [a, b],
* max error of eps, max depth of depth */
DB adaptive_simpsons(DB (*f)(DB), DB a, DB b,
DB eps, int depth) {
DB c = (a + b) / 2, h = b - a;
DB fa = f(a), fb = f(b), fc = f(c);
DB s = (fa + 4 * fc + fb) * h / 6;
return adaptive_simpsons_aux(f, a, b, eps, s, fa, fb, fc, depth);
}
DB f(DB x){
return a * x * x + b * x + c;
}
int n;
const int N = 2e4 + 9;
Point2D p[N];
vector<DB> h;
bool inSeg(Point2D l , Point2D r , Point2D mid){
return dblcmp(mid.x - fmin(l.x , r.x)) >= 0
&& dblcmp(fmax(l.x , r.x) - mid.x) >= 0
&& dblcmp(mid.y - fmin(l.y , r.y)) >= 0
&& dblcmp(fmax(l.y , r.y) - mid.y) >= 0;
}
bool check(Point2D pre , Point2D now , Point2D nex){
pre = pre - now;
nex = nex - now;
pre = pre.trunc(0.1);
nex = nex.trunc(0.1);
pre = pre + now;
nex = nex + now;
if(dblcmp(pre.y - f(pre.x)) * dblcmp(nex.y - f(nex.x)) < 0) return true;
return false;
}
void solve(){
scanf("%lf%lf%lf%lf%lf" , &a , &b , &c , &l , &r);
for (int i = 0 ; i < n ; ++i) p[i].input();
h.resize(0);
if (n <= 2){
puts("0.00");
return;
}
for (int i = 0 ; i < n ; ++i){
Point2D now , nex , pre , Inter , Inter2;
pre = p[(i + n - 1) % n];
now = p[i];
nex = p[(i + 1) % n];
// if (now == nex) continue;
if (dblcmp(now.x - nex.x) == 0){
Inter = Point2D(now.x , f(now.x));
if (inSeg(now , nex , Inter)) {
h.PB(Inter.x);
if (Inter == now) {
if (!check(pre , now , nex)){
h.PB(Inter.x);
}
}
}
}
else{
DB kk = (nex.y - now.y) / (nex.x - now.x);
DB bb = nex.y - nex.x * kk;
DB A = sqr(b - kk) - 4. * (c - bb) * a;
if (dblcmp(A) > 0){
DB x1 = (-(b - kk) + sqrt(A)) / 2. / a;
DB x2 = (-(b - kk) - sqrt(A)) / 2. / a;
Inter = Point2D(x1 , f(x1));
Inter2 = Point2D(x2 , f(x2));
if (inSeg(now , nex , Inter)) {
if (!(Inter == pre) && !(Inter == nex))h.PB(Inter.x);
if (Inter == now && !check(pre , now , nex)) h.PB(Inter.x);
}
if (inSeg(now , nex , Inter2)) {
if (!(Inter2 == pre) && !(Inter2 == nex))h.PB(Inter2.x);
if (Inter2 == now && !check(pre , now , nex))h.PB(Inter2.x);
}
}
}
}
sort(ALL(h));
DB ans = 0;
for (int i = 0 ; i + 1 < SZ(h) ; i += 2){
DB L = fmax(l , h[i]) , R = fmin(h[i + 1] , r);
if (dblcmp(R - L) > 0){
ans += adaptive_simpsons(f0 , L , R , eps , 10);
// ans += simpson(f0 , L , R);
// ans += romberg(f0 , L , R );
}
if (dblcmp(R - L) < 0) break;
}
double out = ans;
printf("%.2lf\n" , out + eps);
}
int main(){
// freopen("in.txt" , "r" , stdin);
// freopen("output.txt" , "w" , stdout);
while(~scanf("%d" , &n)) solve();
}
/**
4 1 0 0 -10 10
4 4
1 1
1 0
0 1
Answer = 1.90
*/