DDART - Editorial

PROBLEM LINK:

Practice
Div-1 Contest

Author: Kasra Mazaheri
Tester: Yash Chandnani
Editorialist: Michael Nematollahi

DIFFICULTY:

HARD

PREREQUISITES:

Geometry, Convex Hull Trick

PROBLEM:

Consider the \mathbb{R}^2 plane. You’re given Q queries. Each query gives you either a circle or a point on the plane.
If a query is a point, you have to determine if it lies inside or on the circumference of all of the circles provided in the previous queries.
You have to answer the queries online.

It is guaranteed that there’s a point P that lies on the circumference of all given circles and that the first three queries are circles.
It is also guaranteed that no two given circles are exactly the same (have the same centre and radius).

QUICK EXPLANATION:

After reading the first two circles, you should be able to determine the two (or one) possible choices for P.
For a point Q at angle \beta w.r.t P to be inside a circle centred at C at angle \alpha w.r.t P and radius r, we must have 2 \times r \times cos(\alpha - \beta) \ge |P-Q|.
Rotate the plane so that the angle of the first circle is 0. Now, if the angle of a point is not in range [-\pi/2, \pi/2], it’s obviously not in the first circle. Otherwise, because two curves of the form a \times cos(\alpha - \beta) intersect at most once in range (-\pi/2, \pi/2), you can maintain a lower hull of these curves efficiently to answer the queries.

EXPLANATION:

Let’s start with finding out what P is.

We know that the first two queries are going to be circles, and that they will not be the same. Hence, they intersect at at most two points. Since we know that P lies on the circumference of all given circles, P has to be one of them.
We assume the first one is P and we’ll keep using it as P until it can no longer be P (in other words, until a given circle does not pass through it).
So from now on, we can assume that we know what P is.

To make thinking about the problem easier, let’s make P the origin by translating the plane about P. From now on we assume P is the origin.

Let C_0 denote the first circle. Rotate the plane so that C_0 is transformed onto the positive side of the X axis.

Let Q be a point with a direction angle \beta. If \beta is not in range [-\pi/2, \pi/2], Q does not lie in the first circle and the answer to the query is “no”. Let’s assume that the angle of Q is in range [-\pi/2, \pi/2].

Let C be the centre of a circle with radius r and let \alpha be the direction angle of the vector represented by C. The point Q lies inside this circle iff 2 \times r \times cos(\alpha - \beta) \ge |Q|. We can interpret this as, "the point (\beta, |Q|) should be below (or on) the curve 2 \times r \times cos(\alpha - \beta)".
Because we know that \beta is in range [-\pi/2, \pi/2] and because two curves of the form above intersect at most once in range (-\pi/2, \pi/2), we can maintain a data structure containing the lower hull of the curves corresponding to the circles received so far. This is similar to maintaining a dynamic convex hull trick.

Maintaining the lower hull would take O(Q \times log(Q)) amortized and checking the value of the lower hull at an angle would take O(log(Q)). So the overall complexity would be O(Q \times log(Q)).

Refer to the tester’s code to see an implementation.

SOLUTIONS:

Setter's Solution
// ItnoE
#include<bits/stdc++.h>
using namespace std;
struct CHT
{
    typedef long double ld;
    typedef pair < ld , ld > Line;
    set < pair < Line , ld > > S;
    set < pair < ld , Line > > I;
    ld INF = 1e18;
    inline void Add(Line X)
    {
        ld t = -INF;
        auto it = S.lower_bound({X, -INF});
        while (it != S.begin())
        {
            it --;
            ld r = Intersection(it->first, X);
            if (r <= it->second)
                I.erase({it->second, it->first}), it = S.erase(it);
            else
                {t = r; break;}
        }
        it = S.lower_bound({X, -INF});
        while (it != S.end())
        {
            Line Y = it->first;
            I.erase({it->second, Y});
            it = S.erase(it);
            ld r = Intersection(X, Y);
            if (r <= t)
            {
                r = -INF;
                if (it != S.begin())
                    it --, r = Intersection(it->first, Y);
                S.insert({Y, r}); I.insert({r, Y});
                return ;
            }
            if (it != S.end() && it->second <= r)
                continue;
            S.insert({Y, r});
            I.insert({r, Y});
            break;
        }
        S.insert({X, t});
        I.insert({t, X});
    }
    inline ld GetMax(ld X)
    {
        if (!I.size()) return (-INF);
        auto it = I.upper_bound({X, {INF, INF}}); it --;
        return (it->second.first * X + it->second.second);
    }
    inline ld Intersection(Line X, Line Y)
    {
        if (X.first == Y.first && X.second <= Y.second)
            return (-INF);
        if (X.first == Y.first)
            return (INF);
        return ((X.second - Y.second) / (Y.first - X.first));
    }
};
namespace Geometry
{
    #define x real()
    #define y imag()
    typedef long double ld;
    typedef complex < ld > point;
    const ld PI = acos(-1.0), DG = PI / 180, eps = 1e-8, INF = 1e47;
    struct Line
    {
        point a, b, ab;
        Line (point _a, point _b) : a(_a), b(_b), ab(b - a) {}
    };
    struct Circle
    {
        ld r;
        point o;
        Circle (point _o, ld _r) : o(_o), r(_r) {}
    };
    inline bool Equal(ld a, ld b)
    {
        return (abs(a - b) < eps);
    }
    inline tuple < ld , ld , ld > LineEq(Line l) // ax + by + c = 0.
    {
        return make_tuple(l.ab.y, -l.ab.x, l.a.y * l.ab.x - l.a.x * l.ab.y);
    }
    inline ld Dot(point a, point b)
    {
        return (a.x * b.x + a.y * b.y);
    }
    inline ld Cross(point a, point b)
    {
        return (a.x * b.y - a.y * b.x);
    }
    inline ld Dist(point a, point b)
    {
        return (sqrt(Dot(a - b, a - b)));
    }
    inline point UnitVector(point a)
    {
        return (a / abs(a));
    }
    inline Line UnitVector(Line l)
    {
        return Line(l.a, l.a + UnitVector(l.ab));
    }
    inline ld Dist(point a, Line l)
    {
        return (Cross(l.ab, a - l.a) / abs(l.ab));
    }
    inline point Image(point a, Line l)
    {
        return (l.a + (l.ab) * Dot(a - l.a, l.ab) / (abs(l.ab) * abs(l.ab)));
    }
    inline ld Angle(Line L1, Line L2) // The first line is the upper one.
    {
        return (arg(L1.ab / L2.ab));
    }
    inline ld Angle(point a, point b, point c) // Angle of abc (a is upper).
    {
        return (arg((a - b) / (c - b)));
    }
    inline point Rotate(point a, ld teta)
    {
        /*  [ cos , -sin ] [x]
            [ sin , cos  ] [y] */
        ld sin_teta = sin(teta), cos_teta = cos(teta);
        return point(cos_teta * a.x - sin_teta * a.y, sin_teta * a.x + cos_teta * a.y);
    }
    inline Line Rotate(Line l, ld teta)
    {
        return Line(l.a, l.a + Rotate(l.ab, teta));
    }
    inline point Intersection(Line L1, Line L2)
    {
        ld cross = Cross(L1.ab, L2.ab);
        if (cross == 0)
            return point(INF, INF);
        return (L1.a + L1.ab * Cross(L2.a - L1.a, L2.ab) / cross);
    }
    inline Line AngleBisector(point a, point b, point c)
    {
        return Line(b, b + UnitVector(a - b) + UnitVector(c - b));
    }
    inline Line LineBisection(Line l)
    {
        return Line((l.a + l.b) * (ld)0.5, (l.a + l.b) * (ld)0.5 + point(-l.ab.y, l.ab.x));
    }
    inline Circle CircumscribedCircle(point a, point b, point c)
    {
        point o = Intersection(LineBisection(Line(a, b)), LineBisection(Line(a, c)));
        if (o.x == o.y && o.x == INF)
            return Circle(o, -1);
        return Circle(o, Dist(o, a));
    }
    inline Circle Incircle(point a, point b, point c)
    {
        point o = Intersection(AngleBisector(a, b, c), AngleBisector(b, a, c));
        if (o.x == o.y && o.x == INF)
            return Circle(o, -1);
        return Circle(o, Dist(o, Line(a, b)));
    }
    inline vector < point > Intersection(Circle c1, Circle c2)
    {
        vector < point > R;
        ld d = Dist(c1.o, c2.o);
        ld cos_teta = (c1.r - c2.r + d * d) / (2.0 * sqrt(c1.r) * d);
        if (cos_teta < -1.0 && cos_teta > -1.0 - eps)
            cos_teta = -1.0;
        if (cos_teta > 1.0 && cos_teta < 1.0 + eps)
            cos_teta = 1.0;
        ld teta = acos(cos_teta);
        Line l = Line(c1.o, c1.o + UnitVector(c2.o - c1.o) * sqrt(c1.r));
        R = {Rotate(l, teta).b, Rotate(l, -teta).b};
        return (R);
    }
}
using namespace Geometry;
vector < point > P;
CHT D, U;
inline void Add(point T)
{
    T -= P[0]; T *= 2;
    ld Ang = arg(T);
    ld d = abs(T);
    d = (ld)1.0 / d;
    T = polar(d, Ang);
    Line L = Line(T, point(0.0, 0.0));
    L = Rotate(L, PI * 0.5);
 
    ld m = L.ab.y / L.ab.x;
    ld b = L.a.y - L.a.x * m;
 
    if (b < 0) // The point should be below the line
        D.Add({- m, - b});
    else // The point should be above the line
        U.Add({m, b});
}
const int N = 1000006;
int n, X[N], Y[N];
long long R[N];
inline void Handle(point T, long long rr)
{
    bool cg = 0;
    ld r = rr; r = sqrt(r);
    while (P.size() && !Equal(Dist(P[0], T), r))
        P.erase(P.begin()), cg = 1;
    assert(P.size());
    if (!cg) return ;
    U.S.clear(); U.I.clear();
    D.S.clear(); D.I.clear();
    for (int i = 1; i < n; i ++)
        Handle(point(X[i], Y[i]), R[i]), Add(point(X[i], Y[i]));
}
int main()
{
    int q, last = 0;
    scanf("%d", &q);
 
    {
        q -= 3;
        int t;
        for (int i = 1; i <= 3; i ++)
            scanf("%d%d%d%lld", &t, &X[i], &Y[i], &R[i]);
 
        long double r1, r2, r3;
        r1 = R[1]; r1 = sqrt(r1);
        r2 = R[2]; r2 = sqrt(r2);
        r3 = R[3]; r3 = sqrt(r3);
 
        vector < point > P1 = Intersection(Circle(point(X[1], Y[1]), R[1]), Circle(point(X[2], Y[2]), R[2]));
        for (auto p1 : P1)
            if (Equal(Dist(p1, point(X[3], Y[3])), r3))
                P.push_back(p1);
 
        assert(P.size());
 
        for (int i = 1; i <= 3; i ++)
            Handle(point(X[i], Y[i]), R[i]), Add(point(X[i], Y[i]));
        n = 4;
    }
 
    for (; q; q --)
    {
        int t;
        scanf("%d%d%d", &t, &X[n], &Y[n]);
        if (last) swap(X[n], Y[n]);
 
        if (t == 1)
        {
            scanf("%lld", &R[n]);
            Handle(point(X[n], Y[n]), R[n]);
            Add(point(X[n], Y[n])); n ++;
        }
        else
        {
            point T = point(X[n], Y[n]) - P[0];
            if (abs(T) <= eps)
            {
                last = 1;
                printf("1\n");
                continue;
            }
            ld Ang = arg(T);
            ld d = abs(T);
            d = (ld)1.0 / d;
            T = polar(d, Ang);
 
            bool bu = U.GetMax(T.x) <= T.y;
            bool bd = - D.GetMax(T.x) >= T.y;
            if (bu && bd) printf("1\n"), last = 1;
            else printf("0\n"), last = 0;
        }
    }
    return 0;
} 
Tester's Solution
#include <bits/stdc++.h>
using namespace std;
 
 
void __print(int x) {cerr << x;}
void __print(long x) {cerr << x;}
void __print(long long x) {cerr << x;}
void __print(unsigned x) {cerr << x;}
void __print(unsigned long x) {cerr << x;}
void __print(unsigned long long x) {cerr << x;}
void __print(float x) {cerr << x;}
void __print(double x) {cerr << x;}
void __print(long double x) {cerr << x;}
void __print(char x) {cerr << '\'' << x << '\'';}
void __print(const char *x) {cerr << '\"' << x << '\"';}
void __print(const string &x) {cerr << '\"' << x << '\"';}
void __print(bool x) {cerr << (x ? "true" : "false");}
 
template<typename T, typename V>
void __print(const pair<T, V> &x) {cerr << '{'; __print(x.first); cerr << ','; __print(x.second); cerr << '}';}
template<typename T>
void __print(const T &x) {int f = 0; cerr << '{'; for (auto &i: x) cerr << (f++ ? "," : ""), __print(i); cerr << "}";}
void _print() {cerr << "]\n";}
template <typename T, typename... V>
void _print(T t, V... v) {__print(t); if (sizeof...(v)) cerr << ", "; _print(v...);}
#ifndef ONLINE_JUDGE
#define debug(x...) cerr << "[" << #x << "] = ["; _print(x)
#else
#define debug(x...)
#endif
 
#define rep(i, n)    for(int i = 0; i < (n); ++i)
#define repA(i, a, n)  for(int i = a; i <= (n); ++i)
#define repD(i, a, n)  for(int i = a; i >= (n); --i)
#define trav(a, x) for(auto& a : x)
#define all(x) x.begin(), x.end()
#define sz(x) (int)(x).size()
#define fill(a)  memset(a, 0, sizeof (a))
#define fst first
#define snd second
#define mp make_pair
#define pb push_back
typedef long double ld;
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
 
void pre(){
 
 
}
 
bool Q;
struct Line {
	mutable ld k, m, p;
	bool operator<(const Line& o) const {
		return Q ? p < o.p : k>o.k;
	}
};
 
struct LineContainer : multiset<Line> {
	// (for doubles, use inf = 1/.0, div(a,b) = a/b)
	const ld inf = 1/.0;
	ld div(ld a, ld b) { // floored division
		return a / b; }
	bool isect(iterator x, iterator y,bool dbg=false) {
		if (y == end()) { x->p = inf; return false; }
		if (x->k==y->k) x->p = x->m > y->m ? inf : -inf;
		else {
			x->p = atan((x->m*cos(x->k)-y->m*cos(y->k))/(y->m*sin(y->k)-x->m*sin(x->k)));
			if(y->m*cos(x->p+0.1-y->k)>x->m*cos(x->p+0.1-x->k)){
				if(x->p<0) x->p+=acos(-1);
				else x->p-=acos(-1);
			}
		}
		return x->p >= y->p;
	}
	void add(ld k, ld m) {
		auto z = insert({k, m, 0}), y = z++, x = y;
		while (isect(y, z)) z = erase(z);
		if (x != begin() && isect(--x, y)) {
			isect(x, y = erase(y));
		}
		while ((y = x) != begin() && (--x)->p >= y->p)
			isect(x, erase(y));
	}
	ld query(ld x) {
		assert(!empty());
		Q = 1; auto l = *lower_bound({0,0,x}); Q = 0;
		return l.m*cos(x-l.k);
	}
};
 
template <class T>
struct Point {
	typedef Point P;
	T x, y;
	explicit Point(T x=0, T y=0) : x(x), y(y) {}
	bool operator<(P p) const { return tie(x,y) < tie(p.x,p.y); }
	bool operator==(P p) const { return tie(x,y)==tie(p.x,p.y); }
	P operator+(P p) const { return P(x+p.x, y+p.y); }
	P operator-(P p) const { return P(x-p.x, y-p.y); }
	P operator*(T d) const { return P(x*d, y*d); }
	P operator/(T d) const { return P(x/d, y/d); }
	T dot(P p) const { return x*p.x + y*p.y; }
	T cross(P p) const { return x*p.y - y*p.x; }
	T cross(P a, P b) const { return (a-*this).cross(b-*this); }
	T dist2() const { return x*x + y*y; }
	ld dist() const { return sqrt((ld)dist2()); }
	// angle to x-axis in interval [-pi, pi]
	ld angle() const { return atan2(y, x); }
	P unit() const { return *this/dist(); } // makes dist()=1
	P perp() const { return P(-y, x); } // rotates +90 degrees
	P normal() const { return perp().unit(); }
	// returns point rotated 'a' radians ccw around the origin
	P rotate(ld a) const {
		return P(x*cos(a)-y*sin(a),x*sin(a)+y*cos(a)); }
};
typedef Point<ld> P;
bool circleIntersection(P a, P b, ld r1, ld r2,
				pair<P, P>& out) {
	P delta = b - a;
	assert(delta.x || delta.y || r1 != r2);
	if (!delta.x && !delta.y) return false;
	ld r = r1 + r2, d2 = delta.dist2();
	ld p = (d2 + r1*r1 - r2*r2) / (2.0 * d2);
	ld h2 = r1*r1 - p*p*d2;
	if (d2 > r*r || h2 < 0) return false;
	P mid = a + delta*p, per = delta.perp() * sqrt(h2 / d2);
	out = {mid + per, mid - per};
	return true;
}
void solve(){
 
 
}
vector<pair<P,ld>> v;
pair<P,P> z;
bool fg1=1,fg2=1;
int main() {
	cin.sync_with_stdio(0); cin.tie(0);
	cin.exceptions(cin.failbit);
	pre();
	int n;cin>>n;
	cerr<<setprecision(20);
	LineContainer L1,L2;
	int gg=0;
	int cnt=1;
	rep(i,n){
		int t;
		ld x,y;
		ll r;cin>>t>>x>>y;
		if(gg) swap(x,y);
		if(t==1) {
			cin>>r;
			v.pb({P(x,y),sqrt(r)});
		}
		if(i==1){
			circleIntersection(v[0].fst,v[1].fst,v[0].snd,v[1].snd,z);
			ld d = v[0].snd*2;
			L1.add(0,d);
			L2.add(0,d);
		}
		if(t==1&&i&&fg1){
			ld d = v.back().snd*2;
			ld ang = (v.back().fst-z.fst).angle() - (v[0].fst-z.fst).angle();
			if(ang>acos(-1)) ang-=2*acos(-1);
			else if(ang<-acos(-1)) ang+=2*acos(-1);
			if(fabs((v.back().fst-z.fst).dist()-v.back().snd)>1e-5) fg1=0;
			else L1.add(ang,d);
		}
		if(t==1&&i&&fg2){
			ld d = v.back().snd*2;
			ld ang = (v.back().fst-z.snd).angle() - (v[0].fst-z.snd).angle();
			if(ang>acos(-1)) ang-=2*acos(-1);
			else if(ang<-acos(-1)) ang+=2*acos(-1);
			if(fabs((v.back().fst-z.snd).dist()-v.back().snd)>1e-5) fg2=0;
			else L2.add(ang,d);
		}
		if(t==2){
			P m(x,y);
			ld ang;
			if(fg1) ang = (m-z.fst).angle() - (v[0].fst-z.fst).angle();
			else ang = (m-z.snd).angle() - (v[0].fst-z.snd).angle();
			if(ang>acos(-1)) ang-=2*acos(-1);
			else if(ang<-acos(-1)) ang+=2*acos(-1);
			if((fg1&&(m-z.fst).dist()<1e-5)||(!fg1&&(m-z.snd).dist()<1e-5)) cout<<1<<'\n',gg=1;
			else if(ang>asin(1)||ang<-asin(1)) cout<<0<<'\n',gg=0;
			else if((fg1&&L1.query(ang)>=(m-z.fst).dist())||(!fg1&&L2.query(ang)>=(m-z.snd).dist())) cout<<1<<'\n',gg=1;
			else cout<<0<<'\n',gg=0;
		}
	}
	return 0;
}
2 Likes

You can also use inversion (centered at the common point of all circles). After inversion circles become half-planes, and the problem is the same: we either add a new half-plane, or ask if a given point belongs to the intersection of all half-planes added so far.

Now this is no longer a dynamic hull of curves, but a dynamic hull of ordinary straight lines, which is much easier to think about.

5 Likes

Hi Krismaz,

what do you mean with inversion? Do you have a resource that explains the concept?

Thanks

Geometric inversion. Wikipedia have plenty of resources.

3 Likes

In short, it’s an easy-to-compute operation that transforms points on the plane; and it has some nice properties. For this problem, an important property is that you can transform a set of points that forms a circle through inversion to get a set of points that forms a line.

You can read about it for example here.

1 Like

Thanks i will look into it.

How did you derive this condition. A diagram will be helpful here.