結果

問題 No.3033 エルハートの数え上げ
ユーザー Kude
提出日時 2025-02-21 22:56:43
言語 C++23
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 4 ms / 2,000 ms
コード長 12,439 bytes
コンパイル時間 4,185 ms
コンパイル使用メモリ 297,136 KB
実行使用メモリ 6,820 KB
最終ジャッジ日時 2025-02-21 22:56:49
合計ジャッジ時間 5,538 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 20
権限があれば一括ダウンロードができます

ソースコード

diff #
プレゼンテーションモードにする

// https://atcoder.jp/contests/ttpc2022/submissions/35187581
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const long long MOD = 998244353ll;
const long long INF = 2305843009213693951ll;
bool comp(pair<ll,ll>& a,pair<ll,ll>& b){
    if(a.first*b.second < b.first*a.second){
        return true;
    }
    return false;
}
pair<long long,long long> reduct(pair<long long,long long> x){
    long long g = gcd(x.first,x.second);
    x.first/=g;x.second/=g;
    if(x.second<0){
        x.first*=-1;
        x.second*=-1;
    }
    return x;
}
int main(void){
    int n, m;
    cin >> n >> m;
    int stats[31][31][31]{};
    for (int i = 0; i < m; i++) {
      ll a, b, c, d;
      cin >> a >> b >> c >> d;
      for (int x = -15; x <= 15; x++) for (int y = -15; y <= 15; y++) for (int z = -15; z <= 15; z++) {
        ll v = a * x + b * y + c * z + d;
        int& s = stats[x+15][y+15][z+15];
        if (v == 0) s++;
        else if (v < 0) s = -1e9;
      }
    }
    
    // ll N,K;
    // cin >> N >> K;
    vector<array<ll,3>> P;
    // for(int i=0;i<N;i++){
        // cin >> P[i][0] >> P[i][1] >> P[i][2];
    // }
    ll K = n;
    for (int x = -15; x <= 15; x++) for (int y = -15; y <= 15; y++) for (int z = -15; z <= 15; z++) {
      if (stats[x+15][y+15][z+15] >= 2) P.push_back({x, y, z});
    }
    ll N = P.size();
    //まず内部が空でない3-単体を探す
    array<ll,4> In = {0,1,-1,-1};
    for(int i=2;i<N;i++){
        array<ll,3> d1 = {P[1][0]-P[0][0],P[1][1]-P[0][1],P[1][2]-P[0][2]};
        array<ll,3> d2 = {P[i][0]-P[0][0],P[i][1]-P[0][1],P[i][2]-P[0][2]};
        array<ll,3> n = {d1[1]*d2[2]-d1[2]*d2[1],d1[2]*d2[0]-d1[0]*d2[2],d1[0]*d2[1]-d1[1]*d2[0]};
        if(n[0]!=0 || n[1]!=0 || n[2]!=0){
            In[2] = i;
            break;
        }
    }
    for(int i=2;i<N;i++){
        if(i==In[2])continue;
        array<ll,3> d1 = {P[1][0]-P[0][0],P[1][1]-P[0][1],P[1][2]-P[0][2]};
        array<ll,3> d2 = {P[In[2]][0]-P[0][0],P[In[2]][1]-P[0][1],P[In[2]][2]-P[0][2]};
        array<ll,3> d3 = {P[i][0]-P[0][0],P[i][1]-P[0][1],P[i][2]-P[0][2]};
        array<ll,3> n = {d1[1]*d2[2]-d1[2]*d2[1],d1[2]*d2[0]-d1[0]*d2[2],d1[0]*d2[1]-d1[1]*d2[0]};
        ll d = d3[0]*n[0] + d3[1]*n[1] + d3[2]*n[2];
        
        if(d!=0){
            In[3] = i;
            break;
        }
    }
    for(int i=0;i<4;i++){
        swap(P[i],P[In[i]]);
    }
    //同一平面上にない4点をとれた
    ll M = 4;
    vector<array<ll,3>> F = {
        {0,1,2},
        {0,1,3},
        {0,2,3},
        {1,2,3}
    };
    
    vector<array<ll,3>> Nor(M);
    //法線ベクトルを求める
    for(int i=0;i<M;i++){
        array<ll,3> d1 = {P[F[i][1]][0]-P[F[i][0]][0],P[F[i][1]][1]-P[F[i][0]][1],P[F[i][1]][2]-P[F[i][0]][2]};
        array<ll,3> d2 = {P[F[i][2]][0]-P[F[i][0]][0],P[F[i][2]][1]-P[F[i][0]][1],P[F[i][2]][2]-P[F[i][0]][2]};
        array<ll,3> n = {d1[1]*d2[2]-d1[2]*d2[1],d1[2]*d2[0]-d1[0]*d2[2],d1[0]*d2[1]-d1[1]*d2[0]};
        Nor[i] = n;
    }
    //法線ベクトルを外向きに統一する
    ll avex=0,avey=0,avez=0;
    for(int i=0;i<4;i++){
        avex += P[i][0];
        avey += P[i][1];
        avez += P[i][2];
    }
    for(int i=0;i<M;i++){
        ll x = 4*P[F[i][0]][0] - avex,y = 4*P[F[i][0]][1] - avey,z = 4*P[F[i][0]][2] - avez;
        ll p = Nor[i][0]*x + Nor[i][1]*y + Nor[i][2]*z;
        if(p<0){
            swap(F[i][1],F[i][2]);
            Nor[i][0]*=-1;Nor[i][1]*=-1;Nor[i][2]*=-1;
        }
    }
    
    //残りのN-4点を凸包につける
    //horizon: i-j辺が地平線にあるかどうか
    int horizon[200][200] = {};
    for(int i=4;i<N;i++){
        // 法線と(面の点の一つ-追加する点)の内積を求める
        // すべて正なら内部にある(零なら同じ平面にある)
        // 負があるならばその面を削除する
        //R:内積
        vector<ll> R(M);
        for(int j=0;j<M;j++){
            array<ll,3> d = {P[F[j][0]][0]-P[i][0],P[F[j][0]][1]-P[i][1],P[F[j][0]][2]-P[i][2]};
            R[j] = Nor[j][0]*d[0] + Nor[j][1]*d[1] + Nor[j][2]*d[2];
            if(R[j]<0){
                horizon[min(F[j][0],F[j][1])][max(F[j][0],F[j][1])] ^= 1;
                horizon[min(F[j][0],F[j][2])][max(F[j][0],F[j][2])] ^= 1;
                horizon[min(F[j][1],F[j][2])][max(F[j][1],F[j][2])] ^= 1;
            }
        }
        vector<array<ll,3>> Ft,Nort;
        
        for(int j=0;j<M;j++){
            if(R[j]<0){
                for(int k=0;k<3;k++){
                    ll u,v;
                    if(k==0){
                        u = F[j][0],v = F[j][1];
                    }else if(k==1){
                        u = F[j][0],v = F[j][2];
                    }else{
                        u = F[j][1],v = F[j][2];
                    }
                    if(horizon[min(u,v)][max(u,v)]!=0){
                        horizon[min(u,v)][max(u,v)] = 0;
                        ll d1 = i, d2 = u, d3 = v;
                        array<ll,3> du = {P[d2][0]-P[d1][0],P[d2][1]-P[d1][1],P[d2][2]-P[d1][2]},dv = {P[d3][0]-P[d1][0],P[d3][1]-P[d1][1],P[d3][2]
                            -P[d1][2]};
                        array<ll,3> d = {du[1]*dv[2]-du[2]*dv[1],du[2]*dv[0]-du[0]*dv[2],du[0]*dv[1]-du[1]*dv[0]};
                        ll p = d[0]*(4*P[d1][0]-avex) + d[1]*(4*P[d1][1]-avey) + d[2]*(4*P[d1][2]-avez);
                        if(p<0){
                            swap(d2,d3);
                            d[0]*=-1;d[1]*=-1;d[2]*=-1;
                        }
                        Ft.push_back({d1,d2,d3});
                        Nort.push_back(d);
                    }
                }
            }
        }
        ll now = 0;
        for(int j=0;j<M;j++){
            if(R[j]>=0){
                F[now] = F[j];
                Nor[now] = Nor[j];
                now++;
            }
        }
        M = now + Ft.size();
        F.resize(M);
        Nor.resize(M);
        for(int j=0;j<(int)Ft.size();j++){
            F[now+j] = Ft[j];
            Nor[now+j] = Nort[j];
        }
    }
    //凸包の体積の8倍
    ll Vol = 0;
    for(int i=0;i<M;i++){
        array<ll,3> d1 = P[F[i][0]],d2 = P[F[i][1]],d3 = P[F[i][2]];
        Vol += d1[0]*d2[1]*d3[2] + d1[1]*d2[2]*d3[0] + d1[2]*d2[0]*d3[1]
             - d1[2]*d2[1]*d3[0] - d1[0]*d2[2]*d3[1] - d1[1]*d2[0]*d3[2];
    }
    //凸包が求まったので、その格子点の個数を求める
    //境界つき内部と境界なし内部が欲しいので、L(内部)とB(境界)を求める
    ll L = 0,B = 0;
    ll Xmin = INF,Xmax = -INF,Ymin = INF,Ymax = -INF,Zmin = INF,Zmax = -INF;
    for(int i=0;i<N;i++){
        Xmin = min(Xmin,P[i][0]);
        Xmax = max(Xmax,P[i][0]);
        Ymin = min(Ymin,P[i][1]);
        Ymax = max(Ymax,P[i][1]);
        Zmin = min(Zmin,P[i][2]);
        Zmax = max(Zmax,P[i][2]);
    }
    //x,y座標を決め打ってzを求める
    //各面について、平面がz軸を含むか含まないかで場合分け
    //もしz軸を含むなら面と決め打ったx,yの共通部分は空または線分　その両端を求める
    //もしz軸を含まないなら共通部分は空または点　上か下かを調べる
    for(ll x=Xmin;x<=Xmax;x++){
        for(ll y=Ymin;y<=Ymax;y++){
            pair<ll,ll> H={-INT_MAX,1},D={INT_MAX,1};
            bool boundary = false;
            for(int i=0;i<M;i++){
                ll a = P[F[i][1]][0]-P[F[i][0]][0], b = P[F[i][1]][1]-P[F[i][0]][1];
                ll c = P[F[i][2]][0]-P[F[i][0]][0], d = P[F[i][2]][1]-P[F[i][0]][1];
                
                if(a*d == b*c){
                    //z軸を含む
                    vector<pair<ll,ll>> W = {
                        {P[F[i][0]][0],P[F[i][0]][1]},
                        {P[F[i][1]][0],P[F[i][1]][1]},
                        {P[F[i][2]][0],P[F[i][2]][1]}
                    };
                    sort(W.begin(),W.end());
                    
                    pair<ll,ll> v1 = {W[2].first-W[0].first,W[2].second-W[0].second};
                    pair<ll,ll> v2 = {x-W[0].first,y-W[0].second};
                    if(v1.first*v2.second-v1.second*v2.first!=0)continue;
                    if(v1.first*v2.first + v1.second*v2.second<0)continue;
                    if(v1.first*v1.first+v1.second*v1.second < v2.first*v2.first+v2.second*v2.second)continue;
                    boundary = true;
                    vector<pair<ll,ll>> E;
                    for(int j=0;j<3;j++){
                        array<ll,3> u = P[F[i][j]], v = P[F[i][(j+1)%3]];
                        ll s1 = x-u[0],s2 = v[0]-x;
                        if(s1==0 && s2==0){s1 = y-u[1];s2 = v[1]-y;}
                        if(s1<0 && s2<0){s1*=-1;s2*=-1;}
                        if(s1<0 || s2<0)continue;
                        if(s1+s2==0){
                            //線分がz軸に平行
                            if(H.first*1 < u[2]*H.second)H = {u[2],1};
                            if(D.first*1 > u[2]*D.second)D = {u[2],1};
                            if(H.first*1 < v[2]*H.second)H = {v[2],1};
                            if(D.first*1 > v[2]*D.second)D = {v[2],1};
                        }else{
                            pair<ll,ll> s = reduct({s2*u[2]+s1*v[2],s1+s2});
                            
                            //upper point
                            if(H.first*s.second < s.first*H.second)H = s;
                            //lower point
                            if(D.first*s.second > s.first*D.second)D = s;
                        }
                    }
                }else{
                    //z軸を含まない
                    
                    //(B-A)s+(C-A)t=((x,y,z)-A) を満たすs,t 実際の値はこれをad-bcで割ったもの
                    ll s = d*(x-P[F[i][0]][0])-c*(y-P[F[i][0]][1]);
                    ll t = -b*(x-P[F[i][0]][0])+a*(y-P[F[i][0]][1]);
                    ll det = a*d-b*c;
                    if(det<0){
                        s*=-1;t*=-1;det*=-1;
                    }
                    if(s<0 || t<0 || s+t>det)continue;
                    pair<ll,ll> z = reduct({P[F[i][0]][2]*(det-s-t)+P[F[i][1]][2]*s+P[F[i][2]][2]*t,det});
                    //upper point
                    if(H.first*z.second < z.first*H.second)H = z;
                    //lower point
                    if(D.first*z.second > z.first*D.second)D = z;
                }
            }
            if(H != pair<ll,ll>{-INT_MAX,1} && D !=pair<ll,ll>{INT_MAX,1}){
                ll Hz,Dz;
                if(H.first>=0){
                    Hz = H.first/H.second;
                }else{
                    Hz = -((-H.first+H.second-1)/H.second);
                }
                if(D.first>=0){
                    Dz = (D.first+D.second-1)/D.second;
                }else{
                    Dz = -((-D.first)/D.second);
                }
                L += Hz-Dz+1;
                H = reduct(H),D = reduct(D);
                if(boundary){
                    B += Hz-Dz+1;
                }else{
                    if(H.second == 1)B++;
                    if(D.second == 1 && H!=D)B++;
                }
            }
        }
    }
    Vol = (Vol*166374059)%MOD;
    L %= MOD;
    B %= MOD;
    //cout << Vol  << " " << L << " " << B << endl;
    /*
    (L_P( 1) - V)   (1  1 1)(c_2)
    (L_P( 0)    ) = (0  0 1)(c_1)
    (L_P(-1) + V)   (1 -1 1)(  1)
    */
    /*
    (c_2)       (1 -2  1)( L - Vol)
    (c_1) = 1/2 (1  0 -1)(    1   )
    (  1)       (0  2  0)(-L+B+Vol)
    */
    //これで内部・内部と境界の格子点の個数と体積が求まったのでエルハート多項式を求める
    //逆行列を埋め込む
    ll a = (L-Vol+MOD)%MOD,b = (-L+B+Vol+MOD)%MOD;
    vector<ll> C(4);
    C[0] = (0*a +       2*1 +       0*b)%MOD;
    C[1] = (1*a +       0*1 + (MOD-1)*b)%MOD;
    C[2] = (1*a + (MOD-2)*1 +       1*b)%MOD;
    C[3] = (2*Vol)%MOD;
    for(int i=0;i<4;i++){
        C[i] = (C[i]*499122177)%MOD;//2で割る
    }
    //エルハート多項式が求まったのでKを代入してK th dilateのdiscrete volumeを求める
    K %= MOD;
    if (K) K = MOD - K;
    ll V = 0;
    for(int i=0;i<4;i++){
        ll num = C[i];
        for(int j=0;j<i;j++){
            num = (num*K)%MOD;
        }
        V = (V + num)%MOD;
    }
    if (V) V = MOD - V;
    cout << V << endl;
    return 0;
}
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
0