Geomerty

preprocess

typedef double db;
const db pi = acos(-1.0);
const db inf = 1e100;
const db eps = 1e-8;
int sign(db a) { return a < -eps ? -1 : a > eps; }
int db_cmp(db a, db b){ return sign(a-b); }
int inmid(db k1,db k2,db k3){return sign(k1-k3)*sign(k2-k3)<=0;}// k3 在 [k1,k2] 内

Point&Vector

struct Point{
    db x,y;
    Point(){}
    Point(db _x,db _y) : x(_x),y(_y){}
    Point operator + (const Point &k1) const{return Point(k1.x+x,k1.y+y);}
    Point operator - (const Point &k1) const{return Point(x-k1.x,y-k1.y);}
    Point operator * (const db k1) const{return Point(x*k1,y*k1);}
    Point operator / (const db k1) const{return Point(x/k1,y/k1);}
    db operator * (const Point b) const{return x * b.x + y * b.y;}//点积
    db operator ^ (const Point b) const{return x * b.y - y * b.x;}//叉积,顺时针为负
    bool operator == (const Point &k1) const{return db_cmp(x,k1.x)==0&&db_cmp(y,k1.y)==0;}
    Point rotate(db k1){return Point(x*cos(k1)-y*sin(k1),x*sin(k1)+y*cos(k1));}// 逆时针旋转 
    Point rotate90(){return Point(-y,x);}
    db abs(){return sqrt(x*x+y*y);}
    db abs2(){return x*x+y*y;}
    Point unit(){return *this/abs();}
    db angle(){return atan2(y,x);}
    bool operator < (const Point k1) const{//极角排序
        int p1=sign(y)==1||(sign(y)==0&&sign(x)==-1),p2=sign(k1.y)==1||(sign(k1.y)==0&&sign(k1.x)==-1);
        return p1<p2||(p1==p2&&sign(*this^k1)>0);
    }
    void out(){printf("%.10f %.10f\n",x,y);}
};
typedef Point Vector;
int inmid(Point k1,Point k2,Point k3){return inmid(k1.x,k2.x,k3.x)&&inmid(k1.y,k2.y,k3.y);}
db angleOfTwoVector(Vector a, Vector b) {return atan2(a ^ b, a * b);}
Point proj(Point q,Point k1,Point k2){ // q 到直线 k1,k2 的投影 
    Point k=k2-k1;
    return k1+k*((q-k1)*k/k.abs2());
}
Point reflect(Point q,Point k1,Point k2){// q 关于直线 k1,k2 的对称点
    return proj(q,k1,k2)*2-q;
}
Point getLL(Point k1,Point k2,Point k3,Point k4){//直线交点
    db w1=(k1-k3)^(k4-k3),w2=(k4-k3)^(k2-k3);
    return (k1*w2+k2*w1)/(w1+w2);
}
int intersect(db l1,db r1,db l2,db r2){//是否有交集
    if (l1>r1) swap(l1,r1); if (l2>r2) swap(l2,r2); return db_cmp(r1,l2)!=-1&&db_cmp(r2,l1)!=-1;
}
int checkSL(Point k1,Point k2,Point k3,Point k4){// 求线段 (S) k1,k2 和直线 (L) k3,k4 的交点 
    return sign((k3-k1)^(k4-k1))*sign((k3-k2)^(k4-k2))<=0;
}
int checkSS(Point k1,Point k2,Point k3,Point k4){//两线段是否相交
    return intersect(k1.x,k2.x,k3.x,k4.x)&&intersect(k1.y,k2.y,k3.y,k4.y)&&checkSL(k1,k2,k3,k4)&&checkSL(k3,k4,k1,k2);
}
db disPP(Point k1,Point k2){return (k2-k1).abs();}
db disPP2(Point k1,Point k2){return (k2-k1).abs2();}
db disPS(Point q,Point k1,Point k2){
    Point k3=proj(q,k1,k2);
    if (inmid(k1,k2,k3)) return disPP(q,k3); 
    return min(disPP(q,k1),disPP(q,k2));
}
db disPL(Point q,Point k1,Point k2){
    if(k1==k2)return disPP(q,k1);
    return fabs((q - k1) ^ (k2 - k1)) / (k2-k1).abs();
}
db disSS(Point k1,Point k2,Point k3,Point k4){
    if (checkSS(k1,k2,k3,k4)) return 0;
    return min(min(disPS(k1,k3,k4),disPS(k2,k3,k4)),min(disPS(k3,k1,k2),disPS(k4,k1,k2)));
}
db disSL(Point k1,Point k2,Point k3,Point k4){
    if (checkSL(k1,k2,k3,k4)) return 0;
    return min(disPL(k1,k3,k4),disPL(k2,k3,k4));
}
int onS(Point q,Point k1,Point k2){return inmid(k1,k2,q)&&sign((k1-q)^(k2-k1))==0;}

PolarPoint

struct PolarPoint {
    db Angle, Length;
    PolarPoint(db a,db b):Angle(a),Length(b){}
};
Vector PolarPointToVector(PolarPoint a) {
    return (Vector){a.Length * cos(a.Angle),a.Length * sin(a.Angle)};
}
PolarPoint VectorToPolarPoint (Vector a){
	return PolarPoint(a.angle(),a.abs());
}

Line&Segment

struct Line {
    Point s, t;
    Line(){}
    Line(Point _s,Point _t):s(_s),t(_t){}
    Vector dir(){return (t-s);}
    Vector unitDir(){return (t-s).unit();}
    int place(Point k){return sign((t-s)^(t-k));}
    db len(){return (t-s).abs();}
};
typedef Line Segment;
Point proj(Point q,Line k){return proj(q,k.s,k.t);}
Point reflect(Point q,Line k){return reflect(q,k.s,k.t);}
Point getLL(Line k1,Line k2){return getLL(k1.s,k1.t,k2.s,k2.t);}
bool parallel(Line k1,Line k2){return sign(k1.dir()^k2.dir())==0;}
bool sameDir(Line k1,Line k2){return parallel(k1,k2)&&sign(k1.dir()*k2.dir())==1;}
int checkSS(Segment k1,Segment k2){return checkSS(k1.s,k1.t,k2.s,k2.t);}
int checkSL(Segment k1,Line k2){return checkSL(k1.s,k1.t,k2.s,k2.t);}
db disPS(Point a, Segment b) {return disPS(a,b.s,b.t);}
db disPL(Point a, Line b) {return disPL(a,b.s,b.t);}
db disSS(Segment k1,Segment k2){return disSS(k1.s,k1.t,k2.s,k2.t);}
db disSL(Segment k1,Line k2){return disSL(k1.s,k1.t,k2.s,k2.t);}
int onS(Point q,Segment k){return onS(q,k.s,k.t);}

Ciercle

struct Circle {
    Point O;
    db R;
    Circle(Point _O,db _R):O(_O),R(_R){}
};

Simpson

用二次曲线逼近原函数，在平面直角坐标系中，由三点 $(x_1,y_1),(x_2,y_2),(x_3,y_3) (x_1<x_2<x_3,x_2=\frac{x_1+x_3} 2)$ 确定的抛物线 $y=f(x)$ 的定积分为

$$\large\int_{x_1}^{x_3}f(x)dx=\frac 1 6(x_3-x_1)(y_1+4y_2+y_3)$$

将要积分的区域 $[L, R]$ 划分成 $n$ 份，横坐标为 $x_0\sim x_n$，对应的函数值分别为 $y_0\sim y_n$ 其结果为

$$\large\int_L^Rf(x)dx\approx \frac {R-L}{3n}[(y_0+y_n)+2(y_2+y_4+\cdots+y_{n-2})+4(y_1+y_3+\cdots+y_{n-1})]$$

自适应Simpson

对要积分的区域 $[L, R]$ 分成两半 $[L, mid]$ 和 $[mid, R]$，如果用二次曲线对 $(L, mid, R)$ 求出的积分值，和对 $(L, \frac {L+mid} 2, mid)$ 和 $(mid,\frac {L+mid} 2, R)$ 分别积分再加起来的值，相差不超过某个精确度，那么就可以停止划分。

// 1
db simpson (db a,db b){
    db c = a + (b-a)/2;
    return (F(a)+4*F(c)+F(b))*(b-a)/6;
}
db asr (db a,db b,db eps,db A){
    db c = a + (b-a)/2;
    db L=simpson(a,c),R=simpson(c,b);
    if(fabs(L+R-A)<=15*eps)return L+R+(L+R-A)/15.0;
    return asr(a,c,eps/2,L)+asr(c,b,eps/2,R);
}
db asr(db a,db b,db eps){//a-左端点, b-右端点
    return asr(a,b,eps,simpson(a,b));
}

// 2
db calc(db len,db fL,db fM,db fR){ //求长度为len的[L,R]区间，中点为M的Simpson近似面积
    return (fL+4*fM+fR)*len/6;
}
db Simpson(db L,db R) {//Simpson积分求区间[L,R]的面积并,F(L)=L,F(R)=R,F(M)=M,把[L,R]当成整体来拟合得到的面积是sqr
    db M=(L+R)/2,fL=F(L),fM=F(M),fR=F(R),sqr=calc(R-L,fL,fM,fR);
    db g1=calc(M-L,fL,F((L+M)/2),fM),g2=calc(R-M,fM,F((M+R)/2),fR);
    if(fabs(sqr-g1-g2)<=eps) //把当前区间分成2半再拟合得到的答案差别很小，就不再递归下去了
        return g1+g2;
    return Simpson(L,M)+Simpson(M,R);
}

Others

n 维球体的体积递推公式

$$V_{n+1}=2 \int_{0}^{\frac{\pi}{2}} V_{n} \cos ^{n} \theta\ d({\sin \theta})$$