普通视图

发现新文章，点击刷新页面。

昨天以前vegetable1024的博客

计算几何模板整理

vegetable1024的博客

作者 vegetable1024

2022年4月5日 20:35

前言

最近快昆明了，所以整理了一下自己做的一些题目，把板子整合了一下。

主要学习了Kuangbin，CF，OI-wiki，还有俊杰的代码。

~~因为自己的板子太杂现在还整合的不太好，有不少冗余的地方不能联系起来。~~

模板

PDF版下载链接：几何模板V2FSign.pdf

#include <bits/stdc++.h>
#define maxn 800005
#define int long long
using namespace std;


// +--------------------------------+
// |                                |
// |  Geometry Template BASIC(Sgn)  |
// |                                |
// +--------------------------------+


//const long long eps = 0;
const long double eps = 1e-8;

//long double情况下使用的sgn函数
int sgn(long double x){
    if(fabs(x) <= eps)    return 0;
    return x > 0 ? 1 : -1;
}
/*
//long long情况下使用的sgn函数
int sgnLL(long long  x){
    if(abs(x) == eps)    return 0;
    return x > 0 ? 1 : -1;
}
*/

// +----------------------------+
// |                            |
// |  Geometry Template Struct  |
// |                            |
// +----------------------------+

/* *** 点（基础模板） *** */
template<typename T> struct TP{
    T x, y; 
    TP(){}
    TP(T _x, T _y){ x = _x; y = _y; }
    TP operator -() const {
        return {-x, -y};
    }
    friend TP operator +(const TP &a, const TP &b){
        return {a.x + b.x, a.y + b.y};
    }
    friend TP operator -(const TP &a, const TP &b){
        return {a.x - b.x, a.y - b.y};
    }
    friend T operator *(const TP &a, const TP &b){
        return a.x * b.x + a.y * b.y;
    }
    friend T operator ^(const TP &a, const TP &b){
        return a.x * b.y - a.y * b.x;
    }
    friend bool operator ==(const TP &a, const TP &b){
        return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;
    }
    friend bool operator <(const TP &a, const TP &b){
        if(sgn(a.x - b.x) == 0)    return sgn(a.y - b.y) < 0;
        return sgn(a.x - b.x) < 0;
    }
    TP operator *(const long double &k) const{
        return {x * k, y * k};
    }
    TP operator /(const long double &k) const{
        return {x / k, y / k};
    }
    //a在逆时针方向: 1，顺时针方向：-1，其他：0
    int toleft(const TP &a) const {
        auto t = (*this) ^ a; 
        return (t > eps) - (t < -eps);
    }
    //返回极角，(-PI, PI]
    long double angle(){
        return atan2(y, x);
    }
    //返回长度
    long double len() const {
        return sqrt(len2());
    }
    //返回两点间距离
    long double dis(const TP &a) const {
        return sqrt(dis2(a));
    }
    //返回长度的平方
    T len2() const {
        return (*this) * (*this);
    }
    //返回两点间距离的平方
    T dis2(const TP &a) const {
        return TP(x - a.x, y - a.y).len2();
    }
    //返回两向量的夹角
    long double ang(const TP &a) const {
        return acos(max(-1.0, min(1.0, ((*this) * a) / (len() * a.len()))));
    }
    //返回逆时针旋转rad后的结果
    TP rot(const long double rad) const {
        return {x * cos(rad) - y * sin(rad), x * sin(rad) + y * cos(rad)};
    }
    //旋转（指定cos和sin）
    TP rot(const long double cosr,const long double sinr) const {
        return {x * cosr - y * sinr, x * sinr + y * cosr};
    }
    //返回长度为R的向量
    TP trunc(long double r){
        long double l = len();
        if(!sgn(l))    return (*this);
        r /= l;
        return {x * r, y * r};
    }
    void input(){
        cin >> x >> y;
    }
    void print(){
        cout << "[Point]\n";
        cout << x << " " << y << '\n';
    }
}; using Point = TP<long double>;

/* *** 线（基础模板） *** */
template<typename T> struct TL{
    TP<T> s, e;
    TL(){}
    TL(TP<T> _s, TP<T> _e){ s = _s; e = _e; }
    friend T operator *(const TL &la, const TL &lb){
        return (la.e - la.s) * (lb.e - lb.s);
    }
    friend T operator ^(const TL &la, const TL &lb){
        return (la.e - la.s) ^ (lb.e - lb.s);
    }
    friend bool operator ==(const TL &la, const TL &lb){
        return la.parallel(lb) && la.isOnSeg(lb.s);
    }
    //点到直线的距离
    long double length(){
        return (e - s).len();
    }
    //点到直线的距离
    long double disLine(const TP<T> &p) const {
        return fabs((p - s) ^ (e - s)) / length();
    }
    //点到线段的距离
    long double disSeg(const TP<T> &p) const{
        if(sgn((p - s) * (e - s)) < 0 || sgn((p - e) * (s - e)) < 0){
            return min(p.dis(s), p.dis(e));
        }
        return disLine(p);
    }
    //点到直线的投影
    TP<T> proj(const TP<T> &p) const {
        return s + (((e - s) * ((e - s) * (p - s))) / ((e - s).len2()));
    }
    //关于直线的对称点
    TP<T> symmetryPoint(TP<T> p){
        Point q = proj(p);
        return Point(2 * q.x - p.x, 2 * q.y - p.y);
    }
    //两直线平行
    bool parallel(TL v){
        return sgn((e - s) ^ (v.e - v.s)) == 0;
    }
    //两直线交点
    TP<T> crosspoint(TL v){
        auto a1 = (v.e - v.s) ^ (s - v.s);
        auto a2 = (v.e - v.s) ^ (e - v.s);
        return {(s.x * a2 - e.x * a1) / (a2 - a1), (s.y * a2 - e.y * a1) / (a2 - a1)};
    }
    //点在线段上，端点：-1，线段内：1，其他：0
    int isOnSeg(const TP<T> &p){
        if(p == s || p == e)    return -1;
        return sgn((p - s) ^ (e - s)) == 0 && sgn((p - s) * (p - e)) <= 0;
    }
    //2 -> 规范相交，1 -> 非规范相交，0 -> 不相交
    int segCrossSeg(TL v){
        int d1 = sgn((e - s) ^ (v.s - s));
        int d2 = sgn((e - s) ^ (v.e - s));
        int d3 = sgn((v.e - v.s) ^ (s - v.s));
        int d4 = sgn((v.e - v.s) ^ (e - v.s));
        if((d1 ^ d2) == -2 && (d3 ^ d4) == -2)    return 2;
        return (d1 == 0 && sgn((v.s - s) * (v.s - e)) <= 0) ||
               (d2 == 0 && sgn((v.e - s) * (v.e - e)) <= 0) ||
               (d3 == 0 && sgn((s - v.s) * (s - v.e)) <= 0) ||
               (d4 == 0 && sgn((e - v.s) * (e - v.e)) <= 0);
    }
    // this -> Line, v -> Seg, 2，规范相交，1，非规范相交，0，不相交
    int lineCrossSeg(TL v){
        int d1 = sgn((e - s) ^ (v.s - s));
        int d2 = sgn((e - s) ^ (v.e - s));
        if((d1 ^ d2) == -2)    return 2;
        return (d1 == 0 || d2 == 0);
    }
    //两直线关系，0，平行，1，重合，2，相交
    int lineRelation(TL v){
        if((*this).parallel(v)){
            return v.toleft(s) == 0;
        }
        return 2;
    }
    //a在直线的，逆时针方向: 1，顺时针方向：-1，其他：0
    int toleft(const TP<T> &p) const {
        int c = sgn((p - s) ^ (e - s));
        if(c < 0)       return 1;
        else if(c > 0)  return -1;
        else            return 0;
    }
    void print(){
        cout << "[Line]\n";
        cout << s.x << " " << s.y << " " << e.x << " " << e.y << '\n';
    }
}; using Line = TL<long double>;

/* *** 圆（基础模板） *** */
const long double PI = acos(-1.0);
template<typename T> struct TC{
    TP<T> c;
    long double r;
    TC(){}
    TC(TP<T> _c, long double _r){ c = _c; r = _r; }
    //根据极角返回圆上一点
    TP<T> point(long double a) {
        return TP<T>(c.x + cos(a) * r, c.y + sin(a) * r);
    }
    long double area(){
        return PI * r * r;
    }
    //5 -> 相离，4 -> 外切， 3 -> 相交，2 -> 内切，1 -> 内含
    int relationCircle(TC v){
        long double d = c.dis(v.c);
        if(sgn((d - r - v.r)) > 0)    return 5;
        if(sgn((d - r - v.r)) == 0)   return 4;
        long double l = fabs(r - v.r);
        if(sgn((d - r - v.r)) < 0 && sgn(d - l) > 0)    return 3;
        if(sgn(d - l) == 0)    return 2;
        if(sgn(d - l) < 0)     return 1;
    }
    //【已测试：ZOJ1597】
    //两圆相交得到的面积
    long double areaCircle(TC v){
        int rel = relationCircle(v);
        if(rel >= 4)    return 0.0;
        if(rel <= 2)    return min(area(), v.area());
        long double d = c.dis(v.c);
        long double hf = (r + v.r + d) / 2.0;
        long double ss = 2 * sqrt(hf * (hf - r) * (hf - v.r) * (hf - d));
        long double a1 = acos((r * r + d * d - v.r * v.r) / (2.0 * r * d));
        a1 = a1 * r * r;
        long double a2 = acos((v.r * v.r + d * d - r * r) / (2.0 * v.r * d));
        a2 = a2 * v.r * v.r;
        return a1 + a2 - ss;
    }
}; using Circle = TC<long double>;

/* *** 多边形（基础模板） *** */
template<typename T> struct TG{
    vector<TP<T>> p;
    size_t nxt(const size_t i) const {return i == p.size() - 1 ? 0 : i + 1;}
    size_t pre(const size_t i) const {return i == 0 ? p.size() - 1 : i - 1;}
    //求面积的二倍（逆时针存点则为正）
    T getArea2(){
        int siz = p.size();
        T sum = 0;
        for(int i = 0; i < siz; i++){
            sum += (p[i] ^ p[(i + 1) % siz]);
        }
        return sum;
    }
    //Winding，判断点与多边形关系，{true, 0} -> 点在边上，{false, cnt} -> 回转数为0 -> 外部，其他 -> 内部
    pair<bool, int> winding(const Point &a) {
        int cnt = 0;
        for(int i = 0; i < p.size(); i++){
            Point u = p[i], v = p[nxt(i)];
            if(sgn((a - u) ^ (a - v)) == 0 && sgn((a - u)*(a - v)) <= 0)    return {true, 0};
            if(sgn(u.y - v.y) == 0)    continue;
            Line uv = {u, v - u};
            if(u.y < v.y - eps && uv.toleft(a) <= 0)    continue;
            if(u.y > v.y + eps && uv.toleft(a) >= 0)    continue;
            if(u.y < a.y - eps && v.y >= a.y - eps)     cnt++;
            if(u.y >= a.y - eps && v.y < a.y - eps)     cnt--;
        }
        return {false, cnt};
    }
    void print(){
        cout << "[Polygon]\n";
        for(int i = 0; i < p.size(); i++){
            cout << i << " " << p[i].x << " " << p[i].y << '\n';
        }
    }
}; using Polygon = TG<long double>;


// +------------------------------+
// |                              |
// |  Geometry Template Function  |
// |                              |
// +------------------------------+


//凸包点集调整 -> 起点变为下凸壳最左侧的点；
void adjustConvexHull(vector<Point> &P, vector<Point> &tmp){
    int n = P.size(); tmp.resize(n);    
    int pos = -1;
    long double minX = 1e60, maxY = -1e60;
    for(int i = 0; i < n; i++){
        if(P[i].x < minX || (fabs(P[i].x - minX) <= eps && P[i].y > maxY)){
            pos = i;
            minX = P[i].x; maxY = P[i].y;
        }
    }
    int cnt = 0;
    for(int i = pos; i < n; i++)    tmp[cnt++] = P[i];
    for(int i = 0; i < pos; i++)    tmp[cnt++] = P[i];
    for(int i = 0; i < n; i++)      P[i] = tmp[i];
}
//求解点集p的凸包（Andrew算法），逆时针存于点集ans中。
#define bk1(x) (x.back())
#define bk2(x) (*(x.rbegin() + 1))
void findConvexHull(vector<Point> p, vector<Point> &ans){
    vector<Point> st;
    sort(p.begin(), p.end(), [&](const Point &A, const Point &B){ return sgn(A.x - B.x) ? A.x < B.x : A.y < B.y; });
    for(Point u : p){
        while(st.size() > 1 && ((bk1(st) - bk2(st)).toleft(u - bk2(st))) <= 0){
            st.pop_back();
        }
        st.push_back(u);
    }
    int k = st.size();
    p.pop_back(); reverse(p.begin(), p.end());
    for(Point u : p){
        while(st.size() > k && ((bk1(st) - bk2(st)).toleft(u - bk2(st))) <= 0){
            st.pop_back();
        }
        st.push_back(u);
    }
    st.pop_back();
    ans.clear();
    for(auto x : st)    ans.push_back(x);
}

//叉积排序函数
bool argcmpC(const Point &a, const Point &b){
    auto Quad = [](const Point &a){
        if(a.y < -eps)    return 1;
        if(a.y > +eps)    return 4;
        if(a.x < -eps)    return 5;
        if(a.x > +eps)    return 3;
        return 2;
    };
    int qa = Quad(a), qb = Quad(b);
    if(qa != qb)    return qa < qb;
    auto cross = (a ^ b);
    return cross > eps;
}

//极角序，自动去重，返回<方向，个数>的集合，需依赖上方的argcmpC函数
vector<pair<Point, int>> polarUniqueTrans(vector<Point> &p){
    map<Point, int, decltype(&argcmpC)> uni{&argcmpC};
    vector<pair<Point, int>> res;
    for(auto x : p)     uni[x]++;
    for(auto x : uni)   res.push_back(x);
    return res;
}
//【已测试：P4525 「模板」自适应辛普森法 1】
//自适应simp积分，近似求int[a, b]
long double functionVal(long double x){
//returns f(x)
    return 0;
}
long double simp(long double l, long double r){
long double mid = (l + r) / 2.0;
return (r - l) * (functionVal(l) + 4 * functionVal(mid) + functionVal(r)) / 6.0;
}
//需注意eps精度问题
long double asr(long double l, long double r, long double ans){
long double mid = (l + r) / 2.0;
long double vL = simp(l, mid), vR = simp(mid, r), tmp = vL + vR - ans;
if(fabs(tmp) <= eps)    return ans;
else                    return asr(l, mid, vL) + asr(mid, r, vR);
}
//求解两圆公切线：返回切线的条数，-1表示无穷多条切线，a -> A上的切点，b -> B上的切点
int getCircleTangents(Circle A, Circle B, vector<Point> &a, vector<Point> &b){
    int cnt = 0;
    if (A.r < B.r) {
        swap(A, B);
        swap(a, b);
    }
    double d2 = (A.c.x - B.c.x) * (A.c.x - B.c.x) + (A.c.y - B.c.y) * (A.c.y - B.c.y);
    double rdiff = A.r - B.r;
    double rsum = A.r + B.r;
    if (sgn(d2 - rdiff * rdiff) < 0) return 0;  // 内含
    double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
    //无限多条切线
    if (sgn(d2) == 0 && sgn(A.r - B.r) == 0) return -1;
    //内切，一条切线
    if (sgn(d2 - rdiff * rdiff) == 0) {  
        a.push_back(A.point(base));
        b.push_back(B.point(base));
        cnt++;
        return cnt;
    }
    //有外公切线
    double ang = acos(rdiff / sqrt(d2));
    a.push_back(A.point(base + ang));
    b.push_back(B.point(base + ang));
    a.push_back(A.point(base - ang));
    b.push_back(B.point(base - ang));
    cnt += 2;
    if(sgn(d2 - rsum * rsum) == 0){  // 一条内公切线
        a.push_back(A.point(base));
        b.push_back(B.point(PI + base));
        cnt++;
    } else if(sgn(d2 - rsum * rsum) > 0){  // 两条内公切线
        double ang = acos(rsum / sqrt(d2));
        a.push_back(A.point(base + ang));
        b.push_back(B.point(PI + base + ang));
        a.push_back(A.point(base - ang));
        b.push_back(B.point(PI + base - ang));
        cnt += 2;
    }
    return cnt;
}

/* 圆的反演：C2C,C2L,L2C */
/*
 * 反演变换
 * 适用于题目中存在多个圆/直线之间的相切关系的情况。
 * 1. 圆O外的点的反演点在圆O内，反之亦然，圆O上的点的反演点为其自身。
 * 2. 不过点O的圆，其反演图形也是不过点O的圆。
 * 3. 过点O的圆，其反演图形是不过点O的直线。
 * 4. 两个图形相切，则他们的反演图形也相切。（*）
 * 5. 两个不经过反演点的外切的圆，反演之后的图形为相交的两条直线。
 *    如果其中一个圆经过反演点，那么反演之后的图形为一个圆和它的一条切线并且反演点和反演后的圆的圆心在切线的【同一侧】。
 *    内切的话反演中心和反演圆的圆心在【异侧】。
 */
//【已测试：HDU4773: Problem of Apollonius】
//点O在圆A外，求圆A的反演圆B，R是反演半径
Circle inversionC2C(Point O, long double R, Circle A){
    long double OA = (A.c - O).len();
    long double RB = 0.5 * ((1 / (OA - A.r)) - (1 / (OA + A.r))) * R * R;
    long double OB = OA * RB / A.r;
    long double Bx = O.x + (A.c.x - O.x) * OB / OA;
    long double By = O.y + (A.c.y - O.y) * OB / OA;
    return Circle(Point(Bx, By), RB);
}
//【已测试：HDU4773: Problem of Apollonius】
//直线反演为过O点的圆B，R是反演半径
Circle inversionL2C(Point O, long double R, Point A, Point B){
    Point P = Line(A, B).proj(O);
    long double d = (O - P).len();
    long double RB = R * R / (2 * d);
    Point VB = (P - O) / d * RB;
    return Circle(O + VB, RB);
}
//【已测试：HDU4773: Problem of Apollonius】
//圆A经过反演中心O，反演得到直线L
Line inversionC2L(Point O, long double R, Circle A){
    long double angle = (O - A.c).angle();
    if(sgn(angle) < 0)    angle += 2 * PI;
    long double angleL = angle + PI / 2;
    long double angleR = angle - PI / 2;
    if(angleL < 0)    angleL += 2 * PI;
    if(angleR < 0)    angleR += 2 * PI;
    Point PL = A.point(angleL), PR = A.point(angleR), dirL = PL - O, dirR = PR - O;
    long double disL = O.dis(PL), disrL = R * R / disL, disR = O.dis(PR), disrR = R * R / disR;
    return Line(O + dirL.trunc(disrL), O + dirR.trunc(disrR));
}
/*
 * 根据两个圆的位置关系来确定情况：
 * (1) 两个圆内含，没有公共点，没有公切线
 * (2) 两圆内切，有一个条公切线
 * (3) 两圆完全重合，有无数条公切线
 * (4) 两圆相交。有2条公切线
 * (5) 两圆外切，有3条公切线
 * (6) 两圆相离，有4条公切线
 */


// +------------------------------+
// |                              |
// |  Geometry Template ExStruct  |
// |                              |
// +------------------------------+


/* 拓展自Polygon：求解点是否在凸包内，以及凸包外一点对该凸包的切线 */
struct Convex : Polygon{ 
    //闵可夫斯基和，对应凸包
    //【已测试：BZOJ2564. 集合的面积】
    Convex operator +(const Convex &c){
        const auto &p = this->p;
        vector<Line> e1(p.size()), e2(c.p.size()), edge(p.size() + c.p.size());
        Convex res; 
        res.p.reserve(p.size() + c.p.size());
        for(int i = 0; i < p.size(); i++){
            e1[i] = {p[i], p[this -> nxt(i)]};
        }
        for(int i = 0; i < c.p.size(); i++){
            e2[i] = {c.p[i], c.p[c.nxt(i)]};
        }
        const auto cmp = [](const Line &u,const Line &v) { return argcmpC(u.e - u.s, v.e - v.s); };
        rotate(e1.begin(), min_element(e1.begin(), e1.end(), cmp), e1.end());
        rotate(e2.begin(), min_element(e2.begin(), e2.end(), cmp), e2.end());
        merge(e1.begin(), e1.end(), e2.begin(), e2.end(), edge.begin(), cmp);
        const auto check = [](const vector<Point> &p, const Point &u){
            const auto back1 = p.back(), back2 = *prev(p.end(), 2);
            return (back1 - back2).toleft(u - back1) == 0 && (back1 - back2) * (u - back1) >= -eps;
        };
        auto u = e1[0].s + e2[0].s;
        for(const auto &v : edge){
            while(res.p.size() > 1 && check(res.p, u)){
                res.p.pop_back();
            }
            res.p.push_back(u);
            u = u + v.e - v.s;
        }
        if(res.p.size() > 1 && check(res.p, res.p[0]))    res.p.pop_back();
        return res;
    }
    //【已测试：Enclosure】
    //O(logN)判断点是否在凸包内，1：在凸包内，0：在凸包外，-1：在凸包上
    int inConvex(const Point &a){
        auto &p = this->p;
        int l = 1, r = (int)(p.size()) - 2;
        while(l <= r){
            auto mid = (l + r) / 2;
            auto t1 = (p[mid] - p[0]).toleft(a - p[0]);
            auto t2 = (p[mid + 1] - p[0]).toleft(a - p[0]);
            if(t1 >= 0 && t2 <= 0){
                if(mid == 1 && Line(p[0], p[mid]).isOnSeg(a))    return -1;
                if(mid + 1 == (int)(p.size()) - 1 && Line(p[0], p[mid + 1]).isOnSeg(a))    return -1;
                if(Line(p[mid], p[mid + 1]).isOnSeg(a))    return -1;
                return (p[mid + 1] - p[mid]).toleft(a - p[mid]) > 0;
            }    
            if(t1 < 0)    r = mid - 1;
            else          l = mid + 1;
        }
        return false;
    }
    //【已测试：USACO03FALL - Beauty Contest G】
    //旋转卡壳，求解内容取决于传入的函数F
    template<typename F> void rotcaliper(const F &func) {
        const auto &p = this->p;
        const auto area = [](const Point &u, const Point &v, const Point &w){ return fabs((w - u) ^ (w - v)); };
        for(int i = 0, j = 1; i < p.size(); i++){
            const auto nxti = this -> nxt(i);
            func(p[i], p[nxti], p[j]);
            while(area(p[this -> nxt(j)], p[i], p[nxti]) >= area(p[j], p[i], p[nxti])){
                j = this -> nxt(j);
                func(p[i], p[nxti], p[j]);
            }
        }
    }
    //【已测试：USACO03FALL - Beauty Contest G】
    //旋转卡壳，求凸包直径（平方），需根据选定类型确定返回值（long long double/long long）
    long double diameter2(){
        const auto &p = this -> p;
        if(p.size() == 1) return 0;
        if(p.size() == 2) return p[0].dis2(p[1]);
        long double ans = 0;
        auto func = [&](const Point &u, const Point &v, const Point &w){
            ans = max(ans, max(w.dis2(u), w.dis2(v)));
        };
        rotcaliper(func);
        return ans;
    }
    //【已测试：Enclosure】
    //O(logN)求解凸包外一点切线（返回其中一个切点），配合tangent使用
    template<typename F> int extreme(const F &dir){
        auto &p = this -> p;
        auto check = [&](const int i){
            return dir(p[i]).toleft(p[this->nxt(i)] - p[i]) >= 0;
        };
        auto dir0 = dir(p[0]);
        auto check0 = check(0);
        if(check0 == 0 && check((int)(p.size()) - 1))    return 0;
        int l = 0, r = p.size() - 1;
        while(l < r){
            auto mid = (l + r) / 2;
            auto checkm = check(mid);
            if(checkm == check0){
                auto t = dir0.toleft(p[mid] - p[0]);
                if((check0 == 0 && t <= 0) || (check0 && t < 0))    checkm ^= 1;
            }
            if(checkm)    l = mid + 1;
            else          r = mid;
        }
        return r;
    }
    //【已测试：Enclosure】
    //凸包外一点切点（返回两个切点下标）
    pair<int, int> tangent(const Point &a){
        int i = extreme([&](const Point &u){ return u - a; });
        int j = extreme([&](const Point &u){ return a - u; });
        return {i, j};
    }
    //【已测试：A highway and the seven dwarfs】
    //求直线与凸包上点的关系
    pair<int, int> tangent(const Line &l, const Point &dir){
        int i = extreme([&](...){ return dir; });
        int j = extreme([&](...){ return -dir; });
        return {i, j};
    }
    //O(logN)求直线是否穿过凸包
    bool isLineCrossConvex(const Line &l, const Point &dir){
        if(p.size() <= 1)    return true;
        if(p.size() == 2)    return l.toleft(p[0]) == l.toleft(p[1]);
        auto t = tangent(l, dir);
        return l.toleft(p[t.first]) == l.toleft(p[t.second]);
    }
};

//【已测试：Enclosure】
/* 拓展自Convex：利用前缀和求解凸包中若干【连续】点构成的小凸包的面积 */
struct sumConvex : Convex{
    vector<long double> sum;
    void init(){
        getSum();
    }
    void getSum(){
        auto &p = this->p;
        vector<long double> a(p.size());
        for(int i = 0; i < p.size(); i++){
            a[i] = p[this->pre(i)] ^ p[i];
        }
        sum.resize(p.size());
        partial_sum(a.begin(), a.end(), sum.begin());
    }
    long double queryTangentSum(const Point &a){
        auto &p = this->p;
        pair<int, int> result = this->tangent(a);
        int l = result.second, r = result.first;
        return querySum(l, r);
    }
    long double querySum(){
        return sum.back();
    }
    long double querySum(int l, int r){
        if(l <= r)    return sum[r] - sum[l] + (p[r] ^ p[l]);
        return sum[p.size() - 1] - sum[l] + sum[r] + (p[r] ^ p[l]);
    }
};
//【已测试：[HNOI2012]射箭】
/* 半平面（单个） */
struct Halfplane : Line{
    Halfplane(){}
    Halfplane(Point _s, Point _e){s = _s; e = _e;}
    //叉积排序，减少精度损失
    bool operator <(const Halfplane &b){
        Point A = e - s, B = b.e - b.s;
        return argcmpC(A, B);   
    }
};

/* 半平面交（集合） */
struct Halfplanes{
    int n, st, ed, que[maxn];
    Point p[maxn]; Halfplane hp[maxn];
    //去重 & 便于判断非法情况
    void unique(){
        int m = 1;
        for(int i = 1; i < n; i++){
            if(!(sgn(hp[i] ^ hp[i - 1]) == 0 && sgn(hp[i] * hp[i - 1]) >= 0)){
                hp[m++] = hp[i];
            }
            else if(sgn((hp[m - 1].e - hp[m - 1].s) ^ (hp[i].s - hp[m - 1].s)) > 0){
                hp[m - 1] = hp[i];
            }
        }
        n = m;
    }
    //True -> 存在半平面交
    bool Halfplaneinsert(){
        sort(hp, hp + n); unique();
        que[st = 0] = 0; que[ed = 1] = 1;
        p[1] = hp[0].crosspoint(hp[1]);
        for(int i = 2; i < n; i++){
            while(st < ed && sgn((hp[i].e - hp[i].s) ^ (p[ed] - hp[i].s)) < 0)    ed--;
            while(st < ed && sgn((hp[i].e - hp[i].s) ^ (p[st + 1] - hp[i].s)) < 0)    st++;
            que[++ed] = i;
            if(hp[i].parallel(hp[que[ed - 1]]))    return false;
            p[ed] = hp[i].crosspoint(hp[que[ed - 1]]);
        }
        while(st < ed && sgn((hp[que[st]].e - hp[que[st]].s) ^ (p[ed] - hp[que[st]].s)) < 0)    ed--;
        while(st < ed && sgn((hp[que[ed]].e - hp[que[ed]].s) ^ (p[st + 1] - hp[que[ed]].s)) < 0)    st++;
        if(st + 1 >= ed)    return false;
        return true;
    }
    //Halfplaneinsert True -> 得到半平面交对应的凸包
    void getConvex(Polygon &con){
        p[st] = hp[que[st]].crosspoint(hp[que[ed]]);
        for(int j = st, i = 0; j <= ed; i++, j++){
            con.p.push_back(p[j]);
        }
    }
    //压入新的半平面
    void push(Halfplane tmp){
        hp[n++] = tmp;
    }
};

//【已测试：70D - Professor's task】
/* 动态凸包（set维护）：需选用前三个点确定BASIC点，而后进行维护，不支持删除 */
Point DBASIC;
bool argcmpB(const Point &A, const Point &B){
    Point p1 = A - DBASIC, p2 = B - DBASIC;
    //此处可以换用叉积版本维护
    long double len1 = p1.len2(), len2 = p2.len2();
    long double ang1 = p1.angle(), ang2 = p2.angle();
    if(sgn(ang1 - ang2) == 0)    return len1 < len2;
    return ang1 < ang2;
}
struct DConvexHull{
    set<Point, decltype(&argcmpB)> Set{&argcmpB};
    void init(const Point &A, const Point &B, const Point &C){
        //A,B,C为任意三点，处理完毕后直接调用Insert即可
        DBASIC = {(A.x + B.x + C.x) / 3,
                  (A.y + B.y + C.y) / 3};
        Set.insert(A); Set.insert(B); Set.insert(C);
    }
    set<Point>::iterator Pre(set<Point>::iterator it){
        if(it == Set.begin())    it = Set.end();
        return --it;
    }
    set<Point>::iterator Nxt(set<Point>::iterator it){
        ++it;
        return it == Set.end() ? Set.begin() : it;
    }
    //询问点是否在凸包内
    bool Query(Point v){
        auto it = Set.lower_bound(v);
        if(it == Set.end())    it = Set.begin();
        Point v1 = v - (*Pre(it));
        Point v2 = (*it) - (*Pre(it));
        return sgn(v1 ^ v2) <= 0;
    }
    //往凸包中插入一个新的点
    void Insert(Point v){
        if(Query(v))    return;
        Set.insert(v);
        auto it = Nxt(Set.find(v));
        while(Set.size() > 3 && sgn((v - (*Nxt(it))) ^ ((*it) - (*Nxt(it)))) <= 0){
            Set.erase(it);
            it = Nxt(Set.find(v));
        }
        it = Pre(Set.find(v));
        while(Set.size() > 3 && sgn((v - (*it)) ^ ((*it) - (*Pre(it)))) >= 0){
            Set.erase(it);
            it = Pre(Set.find(v));
        }
    }
};

signed main(void){
    cout << "Helloworld!\n";
}

CCPC桂林 F. Illuminations II

vegetable1024的博客

作者 vegetable1024

2022年3月25日 14:07

题意

有两个给定的凸包，小凸包严格在大凸包内部，现在要在大凸包的边界上随机放置一个光源，问小凸包的边被光源照到的期望长度为？

思路

这题是比较直观的。由于要算的是期望长度，因此可以单独计算小凸包中每一条边对最终期望的贡献。

对于小凸包上的每一条边$M$，如果能被光源找到，则光源一定被放置在$M$的延长线与大凸包相交的线段的某一侧。

那么，只要我们能够快速求出可能放置光源位置的长度，就可以求出对应的概率，而该概率再乘上$M$的长度就是这一条边对最后期望的贡献。

可以利用凸包本身的单调性求解：维护一个双指针$L,R$，确定当前的延长线与大凸包的哪些边相交。设小凸包当前边为$\vec{u}$，大凸包的两个交点分别为$P_L,P_R$，那么有两种情况：

$\vec{u}$与$\vec{P_LP_R}$同向，则得知当前要求的是从$L$到$R$这一段的距离
$\vec{u}$与$\vec{P_LP_R}$异向，则得知当前要求的是从$R$到$L$这一段的距离。

而后，可以预处理前缀和快速求解这一段完整的边的长度，再依据这两个分类讨论不完整的线段的长度加在一起。

由于可能出现$L>R$的情况，前缀和讨论比较复杂，因此可以事先将存着长度数组复制一遍化环为链，当出现$L>R$的情况时，只要求$[L,R+siz]$这段的区间和就可以了。

另外，在维护双指针的时候，直线与线段相交有三种情况，如果维护不好其中有些情况可能会重合。其中一个比较重要的点是，我们可以规定若直线与线段是非规范相交的，则交点必然是该线段的起点，这样可以使每一情况都是独立的。

代码

#pragma GCC optimize("-fdelete-null-pointer-checks,inline-functions-called-once,-funsafe-loop-optimizations,-fexpensive-optimizations,-foptimize-sibling-calls,-ftree-switch-conversion,-finline-small-functions,inline-small-functions,-frerun-cse-after-loop,-fhoist-adjacent-loads,-findirect-inlining,-freorder-functions,no-stack-protector,-fpartial-inlining,-fsched-interblock,-fcse-follow-jumps,-fcse-skip-blocks,-falign-functions,-fstrict-overflow,-fstrict-aliasing,-fschedule-insns2,-ftree-tail-merge,inline-functions,-fschedule-insns,-freorder-blocks,-fwhole-program,-funroll-loops,-fthread-jumps,-fcrossjumping,-fcaller-saves,-fdevirtualize,-falign-labels,-falign-loops,-falign-jumps,unroll-loops,-fsched-spec,-ffast-math,Ofast,inline,-fgcse,-fgcse-lm,-fipa-sra,-ftree-pre,-ftree-vrp,-fpeephole2",3)
#pragma GCC target("avx", "sse2")
#define IO ios::sync_with_stdio(false);cin.tie(nullptr);cout.tie(nullptr)

#include <bits/stdc++.h>
#define maxn 800005
using namespace std;

const long double eps = 1e-8;
int sgn(long double x){
    if(fabs(x) < eps)    return 0;
    return x > eps ? 1 : -1;
}
template<typename T> struct TP{
    T x, y;
    friend TP operator +(const TP &a, const TP &b){
        return {a.x + b.x, a.y + b.y};
    }
    friend TP operator -(const TP &a, const TP &b){
        return {a.x - b.x, a.y - b.y};
    }
    friend T operator *(const TP &a, const TP &b){
        return a.x * b.x + a.y * b.y;
    }
    friend T operator ^(const TP &a, const TP &b){
        return a.x * b.y - a.y * b.x;
    }
    friend bool operator ==(const TP &a, const TP &b){
        return fabs(a.x - b.x) < eps && fabs(a.y - b.y) < eps;
    }
    int toleft(const TP &b){
        auto p = (*this) ^ b;
        return (p > eps) - (p < -eps);
    }
    T len(){
        return sqrt(x * x + y * y);
    }
    T distance(TP &b){
        return ((*this) - b).len();
    }
    void print(){
        cout << "[Point] " << x << " " << y << '\n';
    }
}; using Point = TP<long double>;

template<typename T> struct TL{
    TP<T> s, e;
    int relation(TP<T> &p){
        int c = sgn((p - s) ^ (e - s));
        if(c < 0)    return 1;
        if(c > 0)    return 2;
        return 3;
    }
    int linecrossseg(TL &v){
        int d1 = sgn((e - s) ^ (v.s - s));
        int d2 = sgn((e - s) ^ (v.e - s));
        if((d1 ^ d2) == -2)    return 2;
        return (d1 == 0 || d2 == 0);
    }   
    TP<T> crosspoint(TL &v){
        auto a1 = (v.e - v.s) ^ (s - v.s);
        auto a2 = (v.e - v.s) ^ (e - v.s);
        return {(s.x * a2 - e.x * a1) / (a2 - a1), (s.y * a2 - e.y * a1) / (a2 - a1)};
    }
}; using Line = TL<long double>;

struct ConvexHull{
    int siz;
    long double alen;
    vector<Point> p; vector<Line> l;
    vector<long double> len;
    void build(vector<Point> &o){
        auto nxt = [&](int i){ return i == p.size() - 1 ? 0 : i + 1; };
        auto pre = [&](int i){ return i == 0 ? p.size() - 1 : i - 1; };
        siz = o.size();
        for(int i = 0; i < siz; i++)    p.push_back(o[i]);
        for(int i = 0; i < siz; i++){
            l.push_back({o[i], o[nxt(i)]});
            len.push_back(o[i].distance(o[nxt(i)]));
            alen += len.back();
        }
    }
}B, S;

int n, m;
vector<Point> iB, iS;
long double querySum(int L, int R, vector<long double> &sum){
    if(R < L)    return 0;
    if(L == 0)   return sum[R];
    return sum[R] - sum[L - 1];
}
long double solve(){
    auto nxt = [&](int i, vector<Line> &p){ return i == p.size() - 1 ? 0 : i + 1; };
    auto pre = [&](int i, vector<Line> &p){ return i == 0 ? p.size() - 1 : i - 1; };
    long double exp = 0;
    vector<long double> sum, val = B.len; vector<Line> out = B.l, in = S.l;
    val.insert(val.end(), val.begin(), val.end());
    sum.resize(val.size());
    partial_sum(val.begin(), val.end(), sum.begin());
    int siz = in.size(), sizo = out.size(), L = 0, R = 0;
    for(int i = 0; i < siz; i++){
        long double result = 0;
        while(!in[i].linecrossseg(out[L]) || 
              (in[i].linecrossseg(out[L]) == 1 && in[i].crosspoint(out[L]) == out[L].e))             L = nxt(L, out);
        while(!in[i].linecrossseg(out[R]) || 
              (in[i].linecrossseg(out[R]) == 1 && in[i].crosspoint(out[R]) == out[R].e) || L == R)   R = nxt(R, out);
        Point crossL = in[i].crosspoint(out[L]);
        Point crossR = in[i].crosspoint(out[R]);
        if((crossR - crossL) * (in[i].e - in[i].s) > 0){
            int FL = L, FR = R;
            if(FL > FR)    FR += sizo;
            long double nowLen = querySum(FL + 1, FR - 1, sum) + 
                                 crossL.distance(out[L].e) +
                                 crossR.distance(out[R].s);
            result += nowLen * S.len[i];
        } else{
            int FL = R, FR = L;
            if(FL > FR)    FR += sizo;
            long double nowLen = querySum(FL + 1, FR - 1, sum) + 
                                 crossR.distance(out[R].e) +
                                 crossL.distance(out[L].s);
            result += nowLen * S.len[i];
        }
        exp += result;
    }
    return exp / B.alen;
}

signed main(void)
{
    scanf("%d %d", &n, &m);
    for(int i = 1; i <= n; i++){
        long double x, y;
        scanf("%Lf %Lf", &x, &y);
        iB.push_back({x, y});
    }
    for(int i = 1; i <= m; i++){
        long double x, y;
        scanf("%Lf %Lf", &x, &y);
        iS.push_back({x, y});
    }
    B.build(iB); S.build(iS);
    printf("%.15Lf\n", solve());
    return 0;
}

HACK数据

这里放一些数据，方便读者自测。

Case#1: Input

Case#1: Output

1	3.888185834356963

Case#2: Input

Case#2: Output

1	4.221320343559643

Case#3: Input

Case#3: Output

1	4.396949086650624

Case#4: Input

本组数据为GYM的Test#3。

8 4
-1000000000 -1
-1 -1000000000
1 -1000000000
1000000000 -1
1000000000 1
1 1000000000
-1 1000000000
-1000000000 1
-1 -1
1 -1
1 1
-1 1

Case#4: Output

1	3.999999997171573

POJ3347 - Kadj Squares

vegetable1024的博客

作者 vegetable1024

2022年2月2日 21:59

题意

给出指定变长的正方形，按输入顺序，将这些正方形以$45°$倾斜角，某一顶点在$x$轴，且后放的正方形与$x$轴相交的顶点的横坐标要大于前放的正方形的限制下依次放置。

问最后由上往下看，能看到多少正方形，其编号是？

思路

一个朴素的想法就是求出所有的正方形的坐标，然后就可以转化为一个类似线段覆盖（将如此放置的正方形的左端点当成线段左端点，对正方形的右端点同理）的问题。

而后，对这些线段的纵坐标从大到小进行排序。若某一线段可以从上往下看到，那么它不会被排序后位于它之前的线段所覆盖，$O(n^2)$进行枚举即可。

蠢方法

为了便于表述，先定义：
边$L$：正方形中，$x$轴相接的顶点和正方形的右端点所连接组成的线段；
边$T$：正方形中，纵坐标最大的点，即上端点，和正方形的右端点所连接组成的线段；
垂线段：正方形的上顶点和下顶点相连接得到的线段

我最初的做法是先求出所有正方形与$x$轴相接的点的坐标，具体过程如下：

首先确定第一个点，设长度数组为$Len$，容易得到第一个点的坐标为$(\frac{Len[1]}{\sqrt{2}}, 0)$
确定下一个正方形$S$与先前的某一正方形$S’$是如何相接的，有三种情况：
① $S$与$S’$的边$L$相接；设$S’$的右端点为$P$，若$P.y - \frac{S.Len}{\sqrt{2}} \geq 0$，则为该种情况；
② $S$与$S’$的边$T$相接；将边$T$延长，得到与$x$轴的交点$G$，而后得到正方形$S$的左端点$G’=(G.x - \frac{S.len}{\sqrt{2}}, G.y+\frac{S.len}{\sqrt{2}})$，若线段$GG’$与先前所有正方形的垂线段都不规范相交，则为该种情况；
③ $S$的左端点直接贴在$y$轴上，不与任何一个先前的正方形有相交之处；若前两种情况均不满足，则为该种情况。
依次求出所有正方形与$x$轴相接点的坐标，存入数组中。

由此即可得到所有正方形的坐标，容易进一步推导得到正方形的左右端点，即转化为线段覆盖问题。

代码：

// 3347.cpp
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <stack>
#include <vector>

#define maxn 100005
using namespace std;
int n;
double l[maxn];
const double eps = 1e-8;

int sgn(double x) {
    if (fabs(x) < eps)
        return 0;

    return x > 0 ? 1 : -1;
}

double sqrtV(double len) {
    return len / sqrt(2.0);
}

struct Point {
    int id;
    double x, y;

    Point() {}

    Point(double _x, double _y) {
        x = _x;
        y = _y;
    }

    Point(double _x, double _y, int _id) {
        x = _x;
        y = _y;
        id = _id;
    }

    Point(Point _p, int _id) {
        x = _p.x;
        y = _p.y;
        id = _id;
    }

    bool operator<(const Point &A) {
        return sgn(x - A.x) == 0 ? y < A.y : x < A.x;
    }

    bool operator==(const Point &A) {
        return sgn(x - A.x) == 0 && sgn(y - A.y) == 0;
    }

    double operator^(const Point &A) {
        return x * A.y - y * A.x;
    }

    double operator*(const Point &A) {
        return x * A.x + y * A.y;
    }

    Point operator-(const Point &A) {
        return Point(x - A.x, y - A.y);
    }

    void print() {
        printf("[Point] x:%.2lf, y:%.2lf id:%d\n", x, y, id);
    }
};

//ty = 1: top to mid, ty = 0: mid to low
struct Line {
    int id, ty;
    Point s, e;

    Line() {}

    Line(Point _s, Point _e) {
        s = _s;
        e = _e;
    }

    Line(Point _s, Point _e, int _id, int _ty) {
        s = _s;
        e = _e;
        id = _id;
        ty = _ty;
    }

    // 2 -> specification intersect
    int segcrossseg(Line v) {
        int d1 = sgn((e - s) ^ (v.s - s));
        int d2 = sgn((e - s) ^ (v.e - s));
        int d3 = sgn((v.e - v.s) ^ (s - v.s));
        int d4 = sgn((v.e - v.s) ^ (e - v.s));

        if ((d1 ^ d2) == -2 && (d3 ^ d4) == -2)
            return 2;

        return (d1 == 0 && sgn((v.s - s) * (v.s - e)) <= 0) ||
               (d2 == 0 && sgn((v.e - s) * (v.e - e)) <= 0) ||
               (d3 == 0 && sgn((s - v.s) * (s - v.e)) <= 0) ||
               (d4 == 0 && sgn((e - v.s) * (e - v.e)) <= 0);
    }
};

struct Segment {
    int id;
    double y;
    double lf, rt;
} SM[maxn];

bool cmp(const Segment &A, const Segment &B) {
    return A.y > B.y;
}

vector<Point> Si;
//mid point
Point calc_m(Point s, double len) {
    double mv = sqrtV(len);
    return Point(s.x + mv, s.y + mv);
}

//top point
Point calc_t(Point s, double len) {
    double mv = sqrtV(len) * 2;
    return Point(s.x, s.y + mv);
}

bool check(Point s, double len) {
    Line newV = Line(s, Point(s.x - sqrtV(len), s.y + sqrtV(len)));

    for (int i = 0; i < Si.size(); i++) {
        Point low = Si[i], up = Point(Si[i].x, 2 * sqrtV(l[Si[i].id]));

        if (newV.segcrossseg(Line(low, up)) == 2)
            return false;
    }

    return true;
}

int main(void) {
    while (cin >> n && n) {
        Si.clear();

        for (int i = 1; i <= n; i++)
            cin >> l[i];

        double init_x = sqrtV(l[1]), init_y = 0;
        //MUST Line(top, mid) or Line(mid, low)
        Si.push_back(Point(init_x, init_y, 1));
        int now = 2;

        while (now <= n) {
            int flag = 0, nsi = now - 1;

            while (!flag && nsi > 0) {
                nsi--;
                Point nowSi = Si[nsi];
                double nowLen = l[nsi + 1];

                //edge: mid to low
                if (sgn(calc_m(nowSi, nowLen).y - sqrtV(l[now])) >= 0) {
                    Point bas = calc_m(nowSi, nowLen);
                    Point ret = Point(bas.x + sqrtV(l[now]), bas.y - sqrtV(l[now]));
                    double delta = ret.y;
                    Si.push_back(Point(ret.x - delta, ret.y - delta, now));
                    now++;
                    flag = 1;
                }
                //edge: top to mid
                else if (check(Point(nowSi.x + 2 * sqrtV(nowLen), 0), l[now])) {
                    Point bas = nowSi;
                    double rl = sqrtV(nowLen);
                    Si.push_back(Point(bas.x + 2 * rl, 0, now));
                    now++;
                    flag = 1;
                }
            }

            //single
            if (!flag) {
                Si.push_back(Point(sqrtV(l[now]), 0, now));
                now++;
            }
        }

        for (int i = 0; i < n; i++) {
            SM[i].lf = Si[i].x - sqrtV(l[Si[i].id]);
            SM[i].rt = Si[i].x + sqrtV(l[Si[i].id]);
            SM[i].y = l[Si[i].id];
            SM[i].id = Si[i].id;
        }

        sort(SM, SM + n, cmp);
        vector<int> ans;
        double maxr = -1e9;

        for (int i = 0; i < n; i++) {
            double nlf = SM[i].lf, nrt = SM[i].rt;

            for (int j = 0; j < i; j++) {
                if (sgn(SM[j].y - SM[i].y) > 0) {
                    if (nlf >= SM[j].lf)
                        nlf = max(nlf, SM[j].rt);

                    if (nrt <= SM[j].rt)
                        nrt = min(nrt, SM[j].lf);
                }
            }

            //cout << nlf << " " << nrt << endl;
            if (sgn(nrt - nlf) > 0) {
                ans.push_back(SM[i].id);
            }
        }

        sort(ans.begin(), ans.end());

        for (int i = 0; i < ans.size(); i++)
            cout << ans[i] << " ";

        cout << '\n';
    }

    return 0;
}

更优的思路

我们不需要求出正方形的下端点坐标，只需要左端点和右端点的坐标就可以转化问题了。

而求解左端点的坐标，可以使用一些技巧推导。

初始时，正方形$S$的左端点坐标设为$0$，因为正方形只能放置在第一象限。
若正方形$S$贴在正方形$S’$边上的话，则有：
$S.lf.x = max(S.lf.x, S’.rt.x-|S.len-S’.len|)$
这里是放大了原来的正方形，规避了小数，且不会影响最后的结果。
遍历所有已经求出的正方形，由于正方形间不能重叠且正方形底部端点的$x$坐标的值严格单调递增，因此选择$S.lf.x$中最大的一个。

最后，可以通过比较左右端点坐标求出哪些会被覆盖（代码见：CSDN - POJ3347 Kadj Squares），也可以在处理的时候顺便处理纵坐标，对其排序便于求解覆盖关系。

推导左端点的示意图如下：

使用这种求解方法可以大大减小码量。

测试数据（仅供参考）

Input

4
3 5 1 4
3
2 1 2
3
2 2 1
4
3 3 3 3
5
1 3 9 3 1
9
6 1 2 1 1 3 5 1 1
5
4 1 3 5 2
5
4 3 1 7 5
11
5 3 3 1 4 1 1 4 2 8 2
15
4 3 15 6 10 9 15 5 1 10 10 14 5 10 8
20
21 21 4 25 22 1 7 23 26 9 9 17 16 27 1 19 10 20 10 2
50
5 18 4 19 7 23 7 18 12 22 3 17 30 7 2 19 10 1 23 24 14 2 2 28 23 6 7 22 21 8 2 17 25 5 6 23 19 4 2 30 17 4 8 8 2 1 26 11 18 3
36
23 25 9 29 8 24 15 1 21 18 28 30 12 13 24 16 21 1 28 5 24 7 7 19 20 29 30 5 2 9 16 25 4 24 15 5
21
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 20
18
2 5 7 4 3 20 3 4 6 1 3 25 6 7 8 10 2 1
0

Output

1 2 4 
1 3 
1 2 3 
1 2 3 4 
3 
1 4 6 7 
1 3 4 5 
1 2 4 5 
1 2 3 5 8 10 
1 3 5 6 7 10 11 12 14 15 
1 2 4 5 8 9 10 11 12 13 14 16 17 18 19 
2 4 6 8 9 10 12 13 16 19 20 21 24 25 28 29 32 33 36 37 40 41 43 47 49 
1 2 4 6 7 9 10 11 12 13 14 15 16 17 19 21 24 25 26 27 31 32 34 35 
1 2 3 4 5 6 7 8 9 10 11 21 
1 2 3 6 12 14 15 16

Codeforces Round #746 (Div. 2) - E. Bored Bakry

vegetable1024的博客

作者 vegetable1024

2021年10月6日 18:22

题意

给定一个数组$\lbrace a_n \rbrace$，求出其中最长的$Good$序列的长度，若无则输出$0$。

当一段序列$\lbrace a_l,a_{l+1},\dots,a_{r-1},a_{r} \rbrace$满足$a_l \& a_{l+1} \& \dots a_{r-1} \& a_r > a_l \oplus a_{l+1}\oplus\dots\oplus a_{r-1} \oplus a_r$时，称为$Good$序列。

思路

刚开始看到这个题的时候，想的是线段树维护与运算的结果，然后再二分？

首先把$X \& Y$拆成$X \& Y=X|Y-X \oplus Y$，异或运算可以用二十个前缀和来维护，而或运算可以开二十颗线段树维护区间最大值来实现（二十个是因为$2^{20} > 1e6$）。

然后考虑，如果要维护一个严格大于的关系，被选中的序列中，决定大小关系的那一位上的$1$一定要是连续的，只要被选中的序列中有某一位是$0$，那么相与之后得到的结果该位也一定是$0$了，因此貌似可以二分。

时间复杂度理论上是$O(nlog^2n)$，对这道题应该是过不去的。

以上均为口胡的猜想，还未写代码实现过…

后面没写出来，去评论区学习了更为精妙的$O(NlogN)$的写法…

首先我们分析一下，什么情况下$\&$运算得到的结果要优于$\oplus$运算。

假设$\&$运算得到的结果的前$k-1$位和$\oplus$得到的结果的前$k-1$位相同，而第$k$位$\&$运算得到$1$，而$\oplus$得到$0$，其余位随意，那么核心就变成了什么情况下第$k$位$\&$运算能得到$1$，而$\oplus$运算能得到$0$。

第K位

对于$\&$运算，任一元素的对应位为$0$，最后$\&$出来的结果的对应位也为$0$，因此若要使得第$k$位为$1$，序列中所有元素的第$k$位也要为$1$。在第$k$位均为$1$的情况下，为了使$\oplus$运算得到$0$，参与运算的元素个数要为偶数。

前K-1位

由上文对于第$K$位情况的分析，我们知道元素个数为偶数。

那么，假设$\&$运算后得到的第$M(1\leq M\leq K-1)$位结果为$1$，就意味着偶数个$1$进行$\&$运算，此时$\oplus$运算得到的第$M$位结果为$0$，与假设前$k-1$位相同冲突，因此$\&$运算得到的第$M$位的结果只能为$0$。

也就是说，前$K-1$位的$\oplus$运算后得到的结果$XOR$，与$\&$运算后得到的结果$AND$相同，且均为$0$。

思路归纳

由此，题目中要找的序列条件可以拆分为：

序列中的元素的前$k-1$位的异或和为$0$
序列中的元素的第$k$位均为$1$
序列中的元素个数一定为偶数

这样，问题就清晰多了，我们考虑枚举位置$k$。由于$2^{20}>1e6$，我们只需要枚举$20$个就可以了。

一般情况下，寻找前$k-1$位的异或和为$0$这一过程所需的时间复杂度是$O(N)$的，考虑优化。

首先对于连续一段数的异或和，我们考虑使用前缀和优化，设$pre[r]$为$1-r$中元素的异或和，则区间$[L,R]$的异或和可表示为$pre[R] \oplus pre[L-1]$。若$[L,R]$为一个合法序列，则有$pre[R] \oplus pre[L-1]=0$。

更进一步的，我们有：$pre[R]=pre[L-1]$。

有了这个条件，我们考虑对某一位置$R$，快速找到最近的$L$，满足$pre[R]=pre[L-1]$。

由于题目的值域很小$(1 \leq val \leq 1e6)$，我们可以开一个桶$last[val]=pos$，代表异或和$val$出现的合法位置$pos$。这样，目前枚举到第$R$个位置的元素时，满足异或和为$0$的序列$[L,R]$的位置$L=last[pre[R]]$。

此外，合法性主要是指枚举的第$k$位在当前序列中必须有连续的$1$，若出现了$0$则需要更新$last$，再更新答案。我们可以借助一个变量$l$，代表距离当前枚举到的位置$r$最近的$0$的位置，以此来保障求的解的合法性。

满足了前两个条件后，还需要满足序列中的元素个数一定为偶数。可以设两个桶，即$last[2][val]$，一个用来放奇数位的答案，一个用来放偶数位的答案。

整理完成之后可以编写代码如下。

Code

主要参考：https://codeforces.com/blog/entry/95583?#comment-847284

/* vegetable1024 | liushan.site | Maybe Lovely? */

#pragma GCC optimize("-fdelete-null-pointer-checks,inline-functions-called-once,-funsafe-loop-optimizations,-fexpensive-optimizations,-foptimize-sibling-calls,-ftree-switch-conversion,-finline-small-functions,inline-small-functions,-frerun-cse-after-loop,-fhoist-adjacent-loads,-findirect-inlining,-freorder-functions,no-stack-protector,-fpartial-inlining,-fsched-interblock,-fcse-follow-jumps,-fcse-skip-blocks,-falign-functions,-fstrict-overflow,-fstrict-aliasing,-fschedule-insns2,-ftree-tail-merge,inline-functions,-fschedule-insns,-freorder-blocks,-fwhole-program,-funroll-loops,-fthread-jumps,-fcrossjumping,-fcaller-saves,-fdevirtualize,-falign-labels,-falign-loops,-falign-jumps,unroll-loops,-fsched-spec,-ffast-math,Ofast,inline,-fgcse,-fgcse-lm,-fipa-sra,-ftree-pre,-ftree-vrp,-fpeephole2",3)
#pragma GCC target("avx", "sse2")
#define IO ios::sync_with_stdio(false);cin.tie(nullptr);cout.tie(nullptr)

#include <bits/stdc++.h>
#define maxn 1000005
#define int long long
using namespace std;

clock_t START_TIME, END_TIME;

//DEBUG TEMPLATE
template<typename T>
void Print(T value){
    std::cout << value << '\n';
}
template<typename Head, typename... Rail>
void Print(Head head, Rail... rail){
    std::cout << head << ", ";
    Print(rail...);
}

int read(){
    int x = 0; char ch = getchar();
    while(ch < '0' || ch > '9')    ch = getchar();
    while(ch >= '0' && ch <= '9')  x = x * 10 + ch - '0', ch = getchar();
    return x;
}

int n, ans, a[maxn], pre[maxn], last[2][maxn];
signed main(void)
{
    n = read();
    //START_TIME = clock();
    for(int i = 1; i <= n; i++)    a[i] = read();
    for(int k = 20; k >= 1; k--){
        memset(last, -1, sizeof(last));
        for(int i = 1; i <= n; i++){
            pre[i] = pre[i - 1] ^ (a[i] >> k);
        }
        last[0][0] = 0;
        int l = 0;
        for(int r = 1; r <= n; r++){
            if(a[r] & (1 << (k - 1))){
                int lans = last[r & 1][pre[r]];
                if(lans >= l)
                    ans = max(ans, r - lans);
            }
            else    l = r;
            if(last[r & 1][pre[r]] < l){
                last[r & 1][pre[r]] = r;
            }
        }
    }
    printf("%lld\n", ans);
    //END_TIME = clock();
    //printf("%.5lf ms\n", (double)(END_TIME - START_TIME) / CLOCKS_PER_SEC * 1000);
    return 0;
}

乱七八糟

寄了，早知道这把去写$C$和$D$了，想$B$想麻了（

2021ICPC网络赛 - L. Euler Function

vegetable1024的博客

作者 vegetable1024

2021年9月27日 16:20

题目概览

思路

考场上作为队伍$ds$选手，写这个题时候把队友演了，虽然知道是处理$25$个线段树，但是一直没调出来（

考后补题，整理如下。

前置知识：

首先观察到最开始时候$1 \leq x_i \leq 100$，而且$100$以内素数只有$25$个，考虑这是一个突破点。

假设说要对区间$[L,R]$内所有数都乘上一个数字$W$，分两种情况：

假设数字为$W=p_1^{k_1} \times p_2^{k_2} \times \dots p_l^{k_l}$，且区间$[L,R]$内所有数均含有质因子$p_1,p_2,\dots,p_l$时
我们考虑分开算每个质因子的贡献：
如：$[2 \times 2 \times 3,2 \times 3,2 \times 3 \times 3]$，即$[12,6,18]$，我们可以拆为：
- $p=2:[\phi(4)=2, \phi(2)=1, \phi(2)=1]$
- $p=3:[\phi(3)=2,\phi(3)=2, \phi(9)=6]$
由于任意两个质因数间一定互质，且$gcd(n,m)=1$时，$\phi(nm)=\phi(n)\phi(m)$，原数组的欧拉函数值的和为各个质因子对应的欧拉函数的积。
即：$[\phi(12)=\phi(4) \times \phi(3), \phi(6)=\phi(2) \times \phi(3), \phi(18)=\phi(2) \times \phi(9)]$
又因为区间$[L,R]$内均含有所乘的数的质因子$p_1,p_2,\dots,p_l$，且当$n=p^k$时，$\phi(n)=(p-1)p^{k-1}$，那么在计算区间乘的贡献时，其实只要对$[L,R]$的区间和乘上一个$W$就可以了，用上面的数组$[12,6,18]$举例如下：
假设$W=6=2 \times 3$，结果数组为$[12 \times 6,6 \times 6,18 \times 6]$，即：$[72,36,108]$，对应的欧拉函数为$[24,12,36]$；
拆开来看，有：
- $p=2:[\phi(4 \times 2) = \phi(4) \times 2 = 4, \phi(2 \times 2) = \phi(2) \times 2 = 2, \phi(2 \times 2)=\phi(2) \times 2 = 2]$
- $p=3:[\phi(3 \times 3)=\phi(3) \times 3=6,\phi(3 \times 3)=\phi(3) \times 3=6, \phi(9 \times 3)=\phi(9) \times 3=18]$
即$[\phi(8)\times\phi(9)=24,\phi(4)\times\phi(9)=12,\phi(4)\times\phi(27)=3]$
也就是说，这种情况下，你只需要对整个区间的和乘上$W$，得到的就是最终答案。

如果区间$[L,R]$内存在数字并不都含有$W$的全部质因数呢？
只要对区间内这些并不含有$W$的全部质因数的数字暴力更新就可以了。
由于$1\leq x_i,w_i \leq 100$，质因数一定在$100$以内，也就是说最多只有$25$个质因子。
考虑最坏情况，$25$个质因子暴力更新，所需要的也只是$25n$次操作而已，之后只要用线段树维护区间乘和区间和就可以了。

代码

/* vegetable1024 | Maybe Lovely? */

#include <bits/stdc++.h>
#define maxn 800005
#define int long long
#define lson(x) (x << 1)
#define rson(x) (x << 1 | 1)
using namespace std;

//DEBUG TEMPLATE
template<typename T>
void Print(T value){
    std::cout << value << '\n';
}
template<typename Head, typename... Rail>
void Print(Head head, Rail... rail){
    std::cout << head << ", ";
    Print(rail...);
}

int read(){
    int x = 0; char ch = getchar();
    while(ch < '0' || ch > '9')    ch = getchar();
    while(ch >= '0' && ch <= '9')  x = x * 10 + ch - '0', ch = getchar();
    return x;
}

int qpow(int a, int b, int m){
    if(b == 0)    return 1;
    int tmp = qpow(a, b / 2, m);
    tmp = (tmp * tmp) % m;
    if(b & 1)    tmp = (tmp * a) % m;
    return tmp;
}

const int mod = 998244353;
//维护区间欧拉函数的和，支持区间乘法
struct SegTreeSum{
    struct Node{
        int lf, rt, sum, lazy;
    } t[maxn];
    void Pushdown(int now){
        if(t[now].lazy > 1){
            t[lson(now)].lazy = (t[now].lazy * t[lson(now)].lazy) % mod;
            t[rson(now)].lazy = (t[now].lazy * t[rson(now)].lazy) % mod;
            t[lson(now)].sum = (t[now].lazy * t[lson(now)].sum) % mod;
            t[rson(now)].sum = (t[now].lazy * t[rson(now)].sum) % mod;
            t[now].lazy = 1;
        }
    }
    void Pushup(int now){
        t[now].sum = (t[lson(now)].sum + t[rson(now)].sum) % mod;
    }
    void Build(int now, int lf, int rt){
        t[now].lf = lf;
        t[now].rt = rt;
        if(lf == rt){
            t[now].sum = t[now].lazy = 1;
            return;
        }
        t[now].lazy = 1;
        int mid = (lf + rt) / 2;
        Build(lson(now), lf, mid);
        Build(rson(now), mid + 1, rt);
        Pushup(now);
    }
    void Update(int now, int lf, int rt, int k){
        if(t[now].lf >= lf && t[now].rt <= rt){
            t[now].sum = (t[now].sum * k) % mod;
            t[now].lazy = (t[now].lazy * k) % mod;
            return;
        }
        if(t[now].lf > rt || t[now].rt < lf){
            return;
        }
        Pushdown(now);
        Update(lson(now), lf, rt, k);
        Update(rson(now), lf, rt, k);
        Pushup(now);
    }
    int Query(int now, int lf, int rt){
        if(t[now].lf >= lf && t[now].rt <= rt){
            return t[now].sum;
        }
        if(t[now].lf > rt || t[now].rt < lf){
            return 0;
        }
        Pushdown(now);
        return (Query(lson(now), lf, rt) + Query(rson(now), lf, rt)) % mod;
    }
} smt;
//计算每个质因子的出现次数最大值和最小值，判断是否要暴力更新
struct SegTreeRM{
    struct Node{
        int lf, rt, minn, maxx, lazy;
    } t[maxn];
    void Pushup(int now){
        t[now].maxx = max(t[lson(now)].maxx, t[rson(now)].maxx);
        t[now].minn = min(t[lson(now)].minn, t[rson(now)].minn);
    }
    void Pushdown(int now){
        if(t[now].lazy){
            t[lson(now)].lazy += t[now].lazy; t[rson(now)].lazy += t[now].lazy;
            t[lson(now)].maxx += t[now].lazy; t[rson(now)].maxx += t[now].lazy;
            t[lson(now)].minn += t[now].lazy; t[rson(now)].minn += t[now].lazy;
            t[now].lazy = 0;
        }
    }
    void Build(int now, int lf, int rt){
        t[now].lf = lf;
        t[now].rt = rt;
        if(lf == rt){
            t[now].maxx = t[now].minn = t[now].lazy = 0;
            return;
        }
        int mid = (lf + rt) / 2;
        Build(lson(now), lf, mid);
        Build(rson(now), mid + 1, rt);
        Pushup(now);
    }
    void Update(int now, int lf, int rt, int p, int c){
        if(t[now].lf >= lf && t[now].rt <= rt){
            if(t[now].maxx == 0){
                t[now].maxx += c; t[now].minn += c; t[now].lazy += c;
                smt.Update(1, t[now].lf, t[now].rt, ((p - 1) * qpow(p, c - 1, mod)) % mod);
            }
            else if(t[now].minn > 0){
                t[now].maxx += c; t[now].minn += c; t[now].lazy += c;
                smt.Update(1, t[now].lf, t[now].rt, qpow(p, c, mod));
            }
            else{
                Pushdown(now);
                Update(lson(now), lf, rt, p, c);
                Update(rson(now), lf, rt, p, c);
                Pushup(now);
            }
            return;
        }
        if(t[now].lf > rt || t[now].rt < lf){
            return;
        }
        Pushdown(now);
        Update(lson(now), lf, rt, p, c);
        Update(rson(now), lf, rt, p, c);
        Pushup(now);
    }
} tr[30];

int n, m, xi[maxn];
bool check(int val){
    if(val == 1)    return false;
    for(int i = 2; i * i <= val; i++)
        if(val % i == 0)    return false;
    return true;
}
int cnt, pri[maxn];
void initp(int lim){
    for(int i = 1; i <= lim; i++)
        if(check(i))    pri[++cnt] = i;
    smt.Build(1, 1, n);
    for(int i = 1; i <= cnt; i++)    tr[i].Build(1, 1, n);
}
void De(int l, int r, int w){
    for(int i = 1; i <= cnt && w > 1; i++){
        int p = pri[i], c = 0;
        while(w % pri[i] == 0){
            w /= pri[i];
            c++;
        }
        if(c)    tr[i].Update(1, l, r, p, c);
    }
}
signed main(void)
{
    n = read(); m = read(); initp(100);
    for(int i = 1; i <= n; i++)    xi[i] = read();
    for(int i = 1; i <= n; i++)    De(i, i, xi[i]);
    for(int i = 1; i <= m; i++)
    {
        int ty = read();
        if(ty == 0){
            int l = read(), r = read(), w = read();
            De(l, r, w);
        }
        else{
            int l = read(), r = read();
            printf("%lld\n", smt.Query(1, l, r));
        }
    }
    return 0;
}