https://kyopro.hateblo.jp/entry/2025/07/07/125501

問題. 面の数

3次元ユークリッド空間中に 2 つの平面 $H_1: z = 1$ と $H_2: z = 2$ がある。
実数列 $d = (d_1, \ldots, d_n)$ と $d' = (d'_1, \ldots, d'_m)$ が与えられる。平面 $H_1, H_2$ に頂点の内角が原点から見て反時計周り順にそれぞれ $d, d'$ となるように凸多角形を構成して、それら2 つの凸多角形を含む $n + m$ 個の頂点からなる凸多面体の面の数としてあり得るものを全列挙せよ。

制約： $3 \le n, m \le 50$ 、 $10^{-9} \le d_i, d_i' < 180$

解法．

凸多面体（凸多角形）の Minkowski-Weyl の定理から凸多面体には、有限点集合の凸包による表現と、線形不等式系としての表現の 2 つがある。この問題では後者の線形不等式系による表現を用いて考察をする。すなわち、凸多面体 $P \subseteq \mathbb{R}^3$ に対して、ある行列 $A \in \mathbb{R}^{k \times 3}$ とベクトル $b \in \mathbb{R}^k$ が存在して、 $P = \{x \in \mathbb{R}^3 \mid A x \le b \}$ となる。ただし、 $k \in \mathbb{Z}_{> 0}$ は線形不等式の数で凸多面体によって異なる。

まずは平面 $H_1, H_2$ にある凸多角形を考える。下の図では凸多角形とそれを表現する線形不等式系（平面なので直線）を描いている。直線の交点が頂点となり、頂点を結ぶ線分が辺である。凸多角形の辺を記述する直線に対して、その直線の法線ベクトルの向きが変わらないように平行移動しても内角は変わらない。右の辺の青色の頂点は直線を左に平行移動したとき赤色の頂点に、右に平行移動したとき緑色の頂点に移動して、それぞれ内角が変わらないことが確認できる。

次に、2 つの凸多角形を固定したときにできる凸多面体 $P$ の面がどのようになるのかを考える。 $P$ は3次元空間にあり構成方法から記述する線形不等式系の各線形不等式は平面となる。 $H_1$ と $H_2$ にある凸多角形をそれぞれ $P_1, P_2$ とする。 $P_2$ の任意の辺 $e$ を1つ考える。 $e$ を含む $P$ の面を記述する $H_2$ とは異なる平面のイメージは下図のようになる。今、 $e$ は赤色の頂点を結ぶ辺とする。候補となる平面と $H_1$ との共通部分は直線となっており、その直線の法線ベクトルの向きが z 軸を無視して $e$ の法線ベクトルの向きと等しくなる。そして、平面を外側から $P_1$ の頂点のいづれかと交差するまで回転する（ $e$ が蝶番、平面がドアで、蝶番で固定されているドアを動かすイメージ）。このとき $H_1$ 上では直線が移動しているように見える。
平面が $P_1$ と初めて接触するとき $P_1$ の頂点は高々 2 つである。これは入力として与えられる凸多角形のどの 3 頂点も同一直線上にないことから分かる。交点が 1 つの場合は図のように面は三角形となる。交点が 2 つの場合は四角形となり、 $H_1$ と $H_2$ 上にある 2 つの辺の法線ベクトルの向きは z 軸を無視して等しくなる。

$P_1$ と $P_2$ のどの辺も上で構成した三角形または四角形が $P$ の面に含まれている。また、2 つの凸多角形の配置を決めたときに、法線ベクトルの向きが等しくなる辺同士からなる面が四角形となり、それ以外は三角形となる。上部と下部を除いてすべての面が三角形となる凸多面体は片方の凸多角形を適当に回転すれば構成でき、そのときの面の数は $n + m + 2$ である（三角形の底辺を対応させると分かる）。
面が四角形となるような凸多面体を列挙するには、2 つの凸多角形の各辺のペアの法線ベクトルの向きが等しいと固定して、残りの辺で法線ベクトルの向きが等しくなるペアが何個あるか求めればよい。もし法線ベクトルの向きが等しい辺の組数が $k$ 個あれば面の数は $n + m + 2 - k$ 個となる。

実装するときの注意点として、各辺の法線ベクトルを求めるときに基準の法線ベクトルから順に加える回転角は内角ではなく「180° - 内角」を加える。
あと入力は倍精度浮動小数点数型に収まるが、入力時と法線ベクトルの角度を求めるときに加算する必要がありそのときに丸め誤差が生じる。ただし、 $d, d'$ の ulp も小さいし加算回数も高々 50 なので適当に評価すればよい。

計算時間： $O\left( n m (n + m) )\right)$

ソースコードを表示

#include <iostream>
#include <vector>
#include <algorithm>

// 加算時の丸目誤差を考慮
bool is_same(double lhs, double rhs) {
    return std::abs(lhs - rhs) < 1e-10;
}

int main() {
    int n, m;

    while (std::cin >> n, n) {
        // Input H1
        std::vector<double> in_ang_h1(n);
        for (auto &in_ang: in_ang_h1) {
            std::cin >> in_ang;
        }

        // Input H2
        std::cin >> m;
        std::vector<double> in_ang_h2(m);        
        for (auto &in_ang: in_ang_h2) {
            std::cin >> in_ang;
        }

        // 面の数の上限（面がすべて三角形の場合）
        const int MAX_NUM_FACE = n + m + 2;
        // exist_num_face[i] = true の場合は面の数が i となる凸多面体が存在する
        std::vector<bool> exist_num_face(MAX_NUM_FACE + 1, false);

        // 面の数が MAX_NUM_FACE を達成する凸多面体が存在する
        exist_num_face[MAX_NUM_FACE] = true;

        // 多角形 H1 と H2 の各辺の法線ベクトル
        std::vector<double> norm1(n + 1), norm2(m + 1);

        // 多角形 H1 と H2 のそれぞれ i 番目と j 番目の辺を対応させる
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < m; ++j) {
                // H1 の各辺の法線ベクトルを設定（i 番目の辺の法線ベクトルは 0°）
                for (int idx = 0; idx < n; ++idx) {
                    norm1[idx + 1] = norm1[idx] + 180.0 - in_ang_h1[(i + idx) % n];
                }
                // H2 の各辺の法線ベクトルを設定（j 番目の辺の法線ベクトルは 0°）
                for (int idx = 0; idx < m; ++idx) {
                    norm2[idx + 1] = norm2[idx] + 180.0 - in_ang_h2[(j + idx) % m];
                }

                // H1 と H2 の辺の法線ベクトルが等しい組数を求める
                int num_match_pair = 0;
                int idx2 = 0;
                for (int idx1 = 0; idx1 < n; ++idx1) {
                    while (idx2 < m && !is_same(norm1[idx1], norm2[idx2]) && norm2[idx2] < norm1[idx1]) {
                        ++idx2;
                    }

                    if (idx2 < m && is_same(norm1[idx1], norm2[idx2])) {
                        ++num_match_pair;
                    }
                }

                // 現在の凸多面体の面の数を設定
                exist_num_face[n + m + 2 - num_match_pair] = true;
            }
        }

        int num = std::count(exist_num_face.begin(), exist_num_face.end(), true); // 出力フォーマットのため
        for (int i = 1; i <= MAX_NUM_FACE; ++i) {
            if (exist_num_face[i]) {
                std::cout << i;
                if (--num) std::cout << " ";
            }
        }
        std::cout << std::endl;
    }

    return 0;
}

用語の定義をしっかりしようとすると量が倍になりそうなので結果説明が雰囲気だけになってしまった。実装が楽な幾何問題はうれしい。