自适应辛普森

最新推荐文章于 2022-10-16 16:44:23 发布

啊宸

最新推荐文章于 2022-10-16 16:44:23 发布

阅读量349

点赞数 1

分类专栏：计算几何

本文链接：https://blog.csdn.net/qq_36836807/article/details/87380487

版权

计算几何专栏收录该内容

6 篇文章 0 订阅

订阅专栏

文章目录

辛普森公式
自适应辛普森
题目

辛普森公式

$\int_a^bf(x)dx \approx \frac{b-a}{6}*(f(a)+4×f(\frac{a+b}{2})+f(b))$

自适应辛普森

比较左右区间辛普森公式算出的值之和和整个区间辛普森公式算出的值，若差较大则递归左右区间。
这里有一个详细讲解的博客
 模板
没有递归eps,luogu模板2里递归eps似乎A不了啊？？做题的时候递归eps也经常出问题。。

//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<vector>
#include<cmath>
#define For(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 Formylove return 0
#define eps 1e-12
typedef long long LL;
typedef double db;
using namespace std;
db a,b,c,d,L,R;

template<typename T> void read(T &x) {
    char ch=getchar(); T f=1; x=0;
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=-1,ch=getchar();
    for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}

db F(db x) { return (c*x+d)/(a*x+b); }
db sim(db l,db r) { return (r-l)/6.0*(F(l)+4.0*F((l+r)/2.0)+F(r)); }

db sim(db l,db r,db rs) {
    db mid=(l+r)/2.0;
    db lv=sim(l,mid),rv=sim(mid,r);
    if(fabs(lv+rv-rs)<=15.0*eps) return lv+rv+(lv+rv-rs)/15.0;
    return sim(l,mid,lv)+sim(mid,r,rv);
}
    
int main() {
    //freopen("1.in","r",stdin);
    //freopen("1.out","w",stdout);
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R);
    printf("%.6lf\n",sim(L,R,sim(L,R)));
    Formylove;
}

题目

bzoj1502: [NOI2005]月下柠檬树

传送门
圆的任意投影仍是圆。
平行光线下投影为本身。
故本题每个圆面的投影是一个和自身全等的圆。
顶点看成r=eps的圆。把每个圆面投影出来，上下两个圆面之间有无数圆面，也就是相邻两个投影圆作公切线，最终投影就是所有圆面和这些公切线组成的图形。
这个图形的面积很不好算，但是它是连续的，可以用辛普森积分。沿x轴把图形切成两半积分，每个横坐标找它在圆上的最高点和所有切线中的最高点中最高的作为函数值计算。

//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<vector>
#include<cmath>
#define For(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 Formylove return 0
const int N=507;
typedef long long LL;
typedef double db;
using namespace std;
int n,cnt;
db alpha,h[N],r[N];

template<typename T> void read(T &x) {
    char ch=getchar(); T f=1; x=0;
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=-1,ch=getchar();
    for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}

#define eps 1e-10
struct pt {
	db x,y;
	pt(){}
	pt(db x,db y):x(x),y(y){}
	friend bool operator <(const pt&A,const pt&B) {
		return A.x<B.x||(A.x==B.x&&A.y<B.y);
	}
}p[N];
pt operator +(const pt&A,const pt&B) { return pt(A.x+B.x,A.y+B.y); }
pt operator -(const pt&A,const pt&B) { return pt(A.x-B.x,A.y-B.y); }
pt operator *(const pt&A,const db&B) { return pt(A.x*B,A.y*B); }
pt operator /(const pt&A,const db&B) { return pt(A.x/B,A.y/B); }
db cross(pt A,pt B) { return A.x*B.y-A.y*B.x; }
db dot(pt A,pt B) { return A.x*B.x+A.y*B.y; }
db lenth(pt A) { return sqrt(dot(A,A)); }

struct Line{
	pt a,b;
	Line(){}
	Line(pt a,pt b):a(a),b(b){}
}L[N];

struct circle{
	 pt o; db r;
	 circle(){}
	 circle(pt o,db r):o(o),r(r){}
	 pt get(db th) { return pt(o.x+r*cos(th),o.y+r*sin(th)); }
}C[N];

void inter(circle A,circle B) {
	if(A.r<B.r) swap(A,B);
	db d=lenth(A.o-B.o);
	if(d<=(A.r-B.r)) return;
	db bs=atan2(B.o.y-A.o.y,B.o.x-A.o.x);
	db th=acos((A.r-B.r)/d);
	L[++cnt]=Line(A.get(bs+th),B.get(bs+th));
	if(L[cnt].b<L[cnt].a) swap(L[cnt].a,L[cnt].b);
}

db F(db x) {
	db rs=0;
	For(i,1,n) if((x>=C[i].o.x-C[i].r)&&(x<=C[i].o.x+C[i].r)) 
		rs=max(rs,sqrt(C[i].r*C[i].r-(x-C[i].o.x)*(x-C[i].o.x)));
	For(i,1,cnt) if((x>=L[i].a.x)&&(x<=L[i].b.x))
		rs=max(rs,fabs(L[i].a.y+(L[i].b.y-L[i].a.y)*(x-L[i].a.x)/(L[i].b.x-L[i].a.x)));
	return rs;
}

db sim(db l,db r) { return (r-l)/6.0*(F(l)+4.0*F((l+r)/2.0)+F(r)); }

db solve(db l,db r,db rs) {
	db mid=(l+r)/2.0;
	db lv=sim(l,mid),rv=sim(mid,r);	
	if(fabs(lv+rv-rs)<=15.0*eps) return lv+rv+(lv+rv-rs)/15.0;
	return solve(l,mid,lv)+solve(mid,r,rv);
}
 
int main() {
    //freopen("1.in","r",stdin);
    //freopen("1.out","w",stdout);
 	scanf("%d%lf",&n,&alpha);
 	For(i,0,n) scanf("%lf",&h[i]);
 	db L=1e9,R=-1e9;
 	For(i,1,n+1) {
 		h[i]+=h[i-1];
 		if(i<=n) scanf("%lf",&r[i]);
 		else r[i]=eps;
 		C[i].o.x=h[i-1]/tan(alpha);
 		C[i].o.y=0;
 		C[i].r=r[i];
 		L=min(L,C[i].o.x-C[i].r);
 		R=max(R,C[i].o.x+C[i].r);
 	}
 	For(i,1,n) 
 		inter(C[i],C[i+1]);
 	printf("%.2lf\n",2.0*solve(L,R,sim(L,R)));
    Formylove;
}

bzoj2178:圆的面积并

传送门
dcoi权限号过期啦！Please contact lydsy2012@163.com!
给定n(<=1000)个圆，求面积并。
一开始脑子是木的，不对称不会做啊。看了一下别人博客，有人说上下轮廓线积分相减就可以了，有道理呀！
然后我就WA了
一个卡掉我的数据：

在这里插入图片描述
？？？什么神仙图
然后参考了一下黄学长的代码，f(a)是x=a的直线切得的圆的线段并，把这个函数积分求面积。
感性理解一下，这个函数也是连续的。
因为中间会有空段，我是用并查集把圆分成若干联通块每个联通块内积分。

//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<vector>
#include<cmath>
#define For(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 Formylove return 0
const int N=10007;
typedef long long LL;
typedef double db;
using namespace std;
int n,cnt,no[N];

template<typename T> void read(T &x) {
    char ch=getchar(); T f=1; x=0;
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=-1,ch=getchar();
    for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}

#define eps 1e-13
int dcmp(db x) { return fabs(x)<eps?0:(x>0?1:-1); }

struct pt {
	db x,y;
	pt(){}
	pt(db x,db y):x(x),y(y){}
	friend bool operator <(const pt&A,const pt&B) {
		return A.x<B.x||(A.x==B.x&&A.y<B.y);
	}
}p[N],fp[N];
pt operator +(const pt&A,const pt&B) { return pt(A.x+B.x,A.y+B.y); }
pt operator -(const pt&A,const pt&B) { return pt(A.x-B.x,A.y-B.y); }
pt operator *(const pt&A,const db&B) { return pt(A.x*B,A.y*B); }
pt operator /(const pt&A,const db&B) { return pt(A.x/B,A.y/B); }
db cross(pt A,pt B) { return A.x*B.y-A.y*B.x; }
db dot(pt A,pt B) { return A.x*B.x+A.y*B.y; }
db lenth(pt A) { return sqrt(dot(A,A)); }

struct circle{
	 pt o; db r;
	 circle(){}
	 circle(pt o,db r):o(o),r(r){}
	 pt get(db th) { return pt(o.x+r*cos(th),o.y+r*sin(th)); }
	 void read() { scanf("%lf%lf%lf",&o.x,&o.y,&r); }
}C[N],tpC[N];

int tot,tag;
db F(db x) {
	db rs=0;
	int cc=0;
	For(i,1,tot) if(dcmp(x-(C[i].o.x-C[i].r))>=0&&dcmp(x-(C[i].o.x+C[i].r))<=0) {
		db tp=sqrt(C[i].r*C[i].r-(x-C[i].o.x)*(x-C[i].o.x));
		fp[++cc]=pt(C[i].o.y-tp,C[i].o.y+tp);
	}
	sort(fp+1,fp+cc+1);
	db ll=fp[1].x,rr=fp[1].y;
	For(i,2,cc) {
		if(fp[i].x>rr) {
			rs+=rr-ll;
			ll=fp[i].x; rr=fp[i].y;
		}
		else rr=max(rr,fp[i].y);
	}
	return rs+(rr-ll);
}

db sim(db l,db r) { return (r-l)/6.0*(F(l)+4.0*F((l+r)/2.0)+F(r)); }

db sim(db l,db r,db rs) {
	db mid=(l+r)/2.0;
	db lv=sim(l,mid),rv=sim(mid,r);	
	if(fabs(lv+rv-rs)<=15.0*eps) return lv+rv+(lv+rv-rs)/15.0;
	return sim(l,mid,lv)+sim(mid,r,rv);
}

int in(circle A,circle B) {
	if(A.r>B.r) return 0;
	return lenth(A.o-B.o)<=B.r-A.r;
}

int xj(circle A,circle B) {
	return lenth(A.o-B.o)<=B.r+A.r;
}

int fa[N];
int find(int x) { return x==fa[x]?x:fa[x]=find(fa[x]); }

void solve() {
	read(n); db ans=0;
	For(i,1,n) C[i].read();
	For(i,1,n) {
		int fl=1;
		For(j,1,n) if(!no[j]&&i!=j&&in(C[i],C[j])) {
			fl=0;  break;
		}
		if(fl) tpC[++cnt]=C[i];
		else no[i]=1;
	}
	For(i,1,cnt) fa[i]=i;
	For(i,1,cnt) For(j,i+1,cnt) if(xj(tpC[i],tpC[j])) {
		int u=find(i),v=find(j);
		if(u!=v) fa[u]=v;
	}
	For(i,1,cnt) if(find(i)==i) {
		tot=0;
		For(j,1,cnt) if(find(j)==i)
			C[++tot]=tpC[j];	
		db l=1e9,r=-1e9;
		For(j,1,tot) l=min(l,C[j].o.x-C[j].r),r=max(r,C[j].o.x+C[j].r);
		ans+=sim(l,r,F((l+r)/2.0));	
	}
	printf("%.3lf\n",ans);
}
 
int main() {
    freopen("2178.in","r",stdin);
    freopen("2178.out","w",stdout);
 	solve();
    Formylove;
}

NOIp2018集训test-9-4:Problem B. path

orz llj

每条边的长度为 $a*x_i+(1 − a)*y_i$ , $a\in[0,1]$ ,求s到t的期望最短路。
设 $a$ 的最短路为 $d i s (a)$ 。

$ans=\int_0^1dis(x) dx$

边长化为 $x_i-y_i)*a+y_i$ 的形式，那么任意路径的长度都可以化为 $k * a + b$ 的形式， $d i s (a)$ 就是若干直线在 $a$ 处的最低点, $d i s (a)$ 的图形即为连续的凸包形状。
于是直接上自适应辛普森。

//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<vector>
#include<cstdio>
#include<queue>
#include<cmath>
#include<set>
#include<map>
#define Formylove return 0
#define For(i,a,b) for(int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(int i=(a);i>=(b);i--)
const int N=807,M=807;
typedef long long LL;
typedef long double db;
using namespace std;
int n,m,s,t,x[M],y[M];

template<typename T>void read(T &x)  {
    char ch=getchar(); x=0; T f=1;
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=-1,ch=getchar();
    for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}

int ecnt,fir[N],nxt[M],to[M]; 
db val[N];
void add(int u,int v) {
    nxt[++ecnt]=fir[u]; fir[u]=ecnt; to[ecnt]=v;
    nxt[++ecnt]=fir[v]; fir[v]=ecnt; to[ecnt]=u;
}

queue<int>que;
db dis[N];
int vis[N];
#define inf 1e18
#define EPS 1e-5
int dcmp(db x) { return fabs(x)<=EPS?0:(x>0?1:-1); }

db spfa() {
    que.push(t);
    For(i,1,n) dis[i]=inf;
    dis[s]=0;
    que.push(s);
    while(!que.empty()) {
        int x=que.front();
        que.pop(); vis[x]=0;
        for(int i=fir[x];i;i=nxt[i]) {
            if(dcmp(dis[to[i]]-dis[x]-val[i])>0) {
                dis[to[i]]=dis[x]+val[i];
                if(!vis[to[i]]) {
                    vis[to[i]]=1;
                    que.push(to[i]);
                }
            }
        }
    } 
    return dis[t];
}


double f(double a) {
    For(i,1,m) {
        val[i*2]=a*x[i]+(1.0-a)*y[i];
        val[i*2-1]=val[i*2];
    }
    return spfa();
}

double sim(double x,double y) {
    double mid=((x+y)/2.0);
    return  (y-x)/6.0*(f(x)+4.0*f(mid)+f(y));
}

double calc(double l,double r) {
    double mid=((l+r)/2.0);
    double tp=sim(l,mid)+sim(mid,r),tpp=sim(l,r);
    if(dcmp(tp-tpp)==0) return tp+(tp-tpp)/15.0;
    else return calc(l,mid)+calc(mid,r);
}
/*
double calc(double l,double r) {
    double mid=((l+r)/2.0);
    db ls=f(l),rs=f(r),ms=f(mid);
    if(fabs(ls+rs-ms*2.0)<=EPS) return ms*(r-l);
    else return calc(l,mid)+calc(mid,r);
}
*/

#define ANS 
int main() {
#ifdef ANS
    freopen("path.in","r",stdin);
    freopen("path.out","w",stdout);
#endif
    read(n); read(m); read(s); read(t);
    For(i,1,m) {
        int u,v;
        read(u); read(v);
        add(u,v); 
        read(x[i]); read(y[i]);
    }
    db ans=calc(0,1);
    printf("%Lf\n",ans);
    Formylove;
}

啊宸

关注

1
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
自适应辛普森

文章目录辛普森公式自适应辛普森题目bzoj1502: [NOI2005]月下柠檬树bzoj2178:圆的面积并NOIp2018集训test-9-4:Problem B. path辛普森公式∫abf(x)dx≈b−a6∗(f(a)+4×f(a+b2)+f(b))\int_a^bf(x)dx \approx \frac{b-a}{6}*(f(a)+4×f(\frac{a+b}{2})+f(b))...
复制链接

扫一扫