公式

O - 面積クエリ / Area Queries 解説 by physics0523


この手の問題では、まずクエリを先読みすることを考えましょう。
理由: 点を偏角ソートすることになるという予想ができますが、点を追加しながら偏角ソートするといった方針を取ろうとすると計算量が大きくなり非効率です。点の偏角ソートの情報は予め持っておきたいので、クエリ先読みをして予め追加される各点を偏角ソートしておきます。

原点 \(O\) 、点 \(A(x_A,y_A)\) 、点 \(B(x_B,y_B)\) がなす三角形の面積の \(2\) 倍を考えましょう。
半直線 \(OA\) を基準として半直線 \(OB\) の偏角 \(\theta\)

  • \(0^{ \circ } \le \theta < 180^{ \circ }\) であれば \(x_Ay_B - x_By_A\)
  • そうでなければ \(x_By_A - x_Ay_B\)

となります。

特定の点 \(A(x_A,y_A)\)\(1\) つ追加した時に面積の \(2\) 倍の総和がいくら増えるかを考えましょう。
先ほど同様に半直線 \(OA\) を基準として \(P\) 中に追加されている点 \(B(x_B,y_B)\) について半直線 \(OB\) の偏角 \(\theta\)

  • \(0^{ \circ } \le \theta < 180^{ \circ }\) である点 \(B\) 全てについて、三角形 \(OAB\) の面積の \(2\) 倍の総和は \(\sum_B (x_Ay_B - x_By_A)\)
    • なお、 \(\sum_B (x_Ay_B - x_By_A) = x_A \sum_B y_B - y_A \sum_B x_B\) と変形できる。
  • そうでない点 \(B\) 全てについて、三角形 \(OAB\) の面積の \(2\) 倍の総和は \(\sum (x_By_A - x_Ay_B)\)
    • こちらも \(\sum (x_By_A - x_Ay_B) = y_A \sum_B x_B - x_A \sum_B y_B\) と変形できる。

よって、結局以下の問題が解ければよいことになります。
ある点 \(A\) が与えられるので、半直線 \(OA\) を基準として、現在 \(P\) に含まれる点 \(B(x_B,y_B)\) であって半直線 \(OB\) の偏角 \(\theta\)\(0^{ \circ } \le \theta < 180^{ \circ }\) であるものについて、 \(\sum_B x_B\)\(\sum_B y_B\) を求めよ。
これは、追加されることになる点全ての集合 \(S\) とそれらを \(O\) を中心に \(180^{\circ}\) 回転させた点の集合 \(S'\) との和集合 \(T\) を予め偏角ソートしておくと、一点更新区間和クエリで対応できます。これは、 ACL の segtree を利用して処理できます。

実装例(C++):

#include<bits/stdc++.h>
#include<atcoder/all>

using namespace std;
using namespace atcoder;

#define mod 998244353

long long op(long long a, long long b) { return (a+b)%mod; }
long long e() { return 0; }

using namespace std;
using pi=pair<int,int>;

typedef struct{
  string typ;
  int x,y;
}query;

int main(){
  ios::sync_with_stdio(false);
  cin.tie(nullptr);
  int q;
  cin >> q;
  vector<query> qv(q);
  set<pi> plis;
  for(int i=0;i<q;i++){
    cin >> qv[i].typ >> qv[i].x >> qv[i].y;

    int x=qv[i].x,y=qv[i].y;
    if(qv[i].typ=="+"){
      plis.insert({x,y});
      plis.insert({-x,-y});
    }
  }

  vector<pi> pv;
  for(auto &nx : plis){pv.push_back(nx);}
  sort(pv.begin(), pv.end(), [](const auto &p1, const auto &p2) {return atan2l(p1.second, p1.first) < atan2l(p2.second, p2.first);});

  int n=pv.size();
  map<pi,int> mp;
  for(int i=0;i<n;i++){mp[pv[i]]=i;}
  segtree<long long, op, e> segx(2*pv.size()),segy(2*pv.size());
  long long res=0;
  for(int i=0;i<q;i++){
    long long x=qv[i].x,y=qv[i].y;
    x+=mod;y+=mod;
    x%=mod;y%=mod;
    
    int tg=mp[{qv[i].x,qv[i].y}],otg=mp[{-qv[i].x,-qv[i].y}];
    int ntg=tg+n;
    if(tg>otg){otg+=n;}
    long long ce;
    long long sx,sy;
    if(qv[i].typ=="+"){
      sx=segx.prod(tg,otg);
      sy=segy.prod(tg,otg);
      res+=((x*sy)%mod);res%=mod;
      res+=(mod-(y*sx)%mod);res%=mod;

      sx=segx.prod(otg,ntg);
      sy=segy.prod(otg,ntg);
      res+=(mod-(x*sy)%mod);res%=mod;
      res+=((y*sx)%mod);res%=mod;

      ce=1;
    }
    else{
      sx=segx.prod(tg,otg);
      sy=segy.prod(tg,otg);
      res+=(mod-(x*sy)%mod);res%=mod;
      res+=((y*sx)%mod);res%=mod;

      sx=segx.prod(otg,ntg);
      sy=segy.prod(otg,ntg);
      res+=((x*sy)%mod);res%=mod;
      res+=(mod-(y*sx)%mod);res%=mod;

      ce=mod-1;
    }

    segx.set(tg,(segx.get(tg)+ce*x)%mod);
    segx.set(n+tg,(segx.get(n+tg)+ce*x)%mod);

    segy.set(tg,(segy.get(tg)+ce*y)%mod);
    segy.set(n+tg,(segy.get(n+tg)+ce*y)%mod);
    cout << res << '\n';
  }
  return 0;
}

投稿日時:
最終更新: