CAD in Mathematica

　今回は，Mathematica の CAD 関数についてのお話ですが，まずは用語の確認から，．．．，R^{n+1}(n≧1)の部分集合 X が Cylindrical であるとは，R^n の部分集合 Y と Y 上で定義された無限大を含めた実数値関数 f,g が存在して

X={(x_1,...,x_n,x_{n+1}) ; (x_1,...,x_n)∈Y ∧ f(x_1,...,x_n)≦または＜x_{n+1}≦または＜g(x_1,...,x_n) }

と表せることであり，R の部分集合 X が Cylindrical であるとは，X が R の区間であることでした．そして，R^n(n≧1)の部分集合 X を有限個の Cylindrical な集合，しかも，そこに現れる全ての Y も Cylindrical であり，以下同様であるように，分割したものを X の Cylindrical Decomposition と呼び，その f,g の全てが代数関数であるものを Cylindrical Algebraic Decomposition と呼ぶのでした．
　さて，前回の Maple 上の CylindricalAlgebraicDecompose は，多項式のリストに対して，全空間を，そのリストの各多項式の符号を保つ Cylindrical な部分集合に分解する，いわゆる full CAD でしたが，Mathematica の CAD 関数 CylindricalDecomposition（CylindricalDecomposition—Wolfram Language Documentation）は，論理式のリストに対して，そのリストの各論理式の連言の解集合を，その連言の真理値を保つ Cylindrical な部分集合に分解する，いわゆる partial CAD です（ QEPCAD B や RedLog の CAD もこれです）．
　使い方は簡単で

CylindricalDecomposition[ {論理式のリスト} , {自由変数，束縛変数のリスト} ]

と入力するだけで，その論理式の連言と等価な式を，Cylindrical な各集合を定める条件の選言として出力します．変数の消去（射影）は，リストで指定した右側の変数から順に実行されます．例えば

In[1]:= CylindricalDecomposition[ {a*x + b == 0}, {a, b, x} ]

Out[1]= (a < 0 && x == -(b/a)) ||
(a == 0 && b == 0) ||
(a > 0 && x == -(b/a))

In[2]:= CylindricalDecomposition[ {a*x + b == 0}, {b, a, x} ]

Out[2]= (b < 0 && ( (a < 0 && x == -(b/a)) || (a > 0 && x == -(b/a)))) ||
(b == 0 && ( (a < 0 && x == 0) || a == 0 || (a > 0 && x == 0))) ||
(b >0 && ( (a < 0 && x == -(b/a)) || (a > 0 && x == -(b/a))))

といった具合です．
　この CylindricalDecomposition の正式実装は Mathematica 5 ですが，Mathematica 6 において，指定した段階の射影や open CAD も出力できる GenericCylindricalDecomposition（http://reference.wolfram.com/mathematica/ref/GenericCylindricalDecomposition.ja.html）が追加されました．
　その書式は

GenericCylindricalDecomposition[ {論理式のリスト} , {自由変数，束縛変数のリスト} , {自由変数，束縛変数のリスト} ]

であり，論理式に現れる変数は，第2，3引数のリストの少なくとも一方に属さねばなりません．まず，第2引数のリストが空ならば，出力は CylindricalDecomposition と同じです．逆に，第3引数のリストが空ならば，出力は

{ 解集合の内部の CAD , 解集合の境界の拡大集合の条件 }

であり，空ならそれぞれ False となります．このとき，入力された式を A，出力を { B , C } とおくと

¬C → ( A ←→ B )

が成り立っています．例えば，

In[3]:= CylindricalDecomposition[{x^2 + y^2 + z^2 <= 1}, {x, y, z}]

Out[3]= (x == -1 && y == 0 &&
z == 0) || (-1 < x <
1 && ( ( y == -Sqrt[1 - x^2] &&
z == -Sqrt[1 - x^2 - y^2]) || (-Sqrt[1 - x^2] < y < Sqrt[
1 - x^2] && -Sqrt[1 - x^2 - y^2] <= z <= Sqrt[
1 - x^2 - y^2]) || (y == Sqrt[1 - x^2] &&
z == -Sqrt[1 - x^2 - y^2]))) || (x == 1 && y == 0 && z == 0)

In[4]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 <=
1}, {}, {x, y, z}]

Out[4]= (x == -1 && y == 0 &&
z == 0) || (-1 < x <
1 && ( ( y == -Sqrt[1 - x^2] &&
z == -Sqrt[1 - x^2 - y^2]) || (-Sqrt[1 - x^2] < y < Sqrt[
1 - x^2] && -Sqrt[1 - x^2 - y^2] <= z <= Sqrt[
1 - x^2 - y^2]) || (y == Sqrt[1 - x^2] &&
z == -Sqrt[1 - x^2 - y^2]))) || (x == 1 && y == 0 && z == 0)

In[5]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 <= 1}, {x,
y, z}, {}]

Out[5]= {-1 < x < 1 && -Sqrt[1 - x^2] < y < Sqrt[
1 - x^2] && -Sqrt[1 - x^2 - y^2] < z < Sqrt[1 - x^2 - y^2], -1 +
x^2 + y^2 + z^2 == 0}

In[6]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 < 1}, {x,
y, z}, {}]

Out[6]= {-1 < x < 1 && -Sqrt[1 - x^2] < y < Sqrt[
1 - x^2] && -Sqrt[1 - x^2 - y^2] < z < Sqrt[1 - x^2 - y^2],
False }

In[7]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 == 1}, {x,
y, z}, {}]

Out[7]= { False, -1 + x^2 + y^2 + z^2 == 0 }

といった具合です．対外的には GenericCylindricalDecomposition は，その open CAD 出力を以って，CylindricalDecomposition の "内部による近似版"，"高速版" という位置付けのようで，確かに

In[8]:=
GenericCylindricalDecomposition[{x^2 + y^2 + z^2 < 1 &&
a*x + b*y + c*z > 1}, {a, b, c, x, y, z}, {}]; // Timing

Out[8]= {0.796, Null}

ですが，CylindricalDecomposition は 1 時間以内には停止しませんので，体積や可視化などという恐ろしく贅沢な CAD の使い方にはピッタリなのかも知れません．
　一方，本来の QE の立場から注目すべきは，何れの引数のリストも空でない入力に対する出力です．この場合の出力も 2 つの成分からなるリストですが，例えば

In[9]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 <= 1},
{x, y}, {z}]

Out[9]= {-1 < x < 1 && -Sqrt[1 - x^2] < y < Sqrt[
1 - x^2] && -Sqrt[1 - x^2 - y^2] <= z <= Sqrt[1 - x^2 - y^2],
1 + x == 0 || -1 + x == 0 || -1 + x^2 + y^2 == 0}

のように，{x,y},{z} なら y,x の評価が open，同様に {x},{y,z} なら x のみが open で評価されます．さらに

In[10]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 <= 1},
{x, y}, {y, z}]

Out[10]= {-1 < x < 1 && -Sqrt[1 - x^2 - y^2] <= z <= Sqrt[
1 - x^2 - y^2], -1 + x == 0 || 1 + x == 0}

In[11]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 <= 1},
{x, y}, {x, y, z}]

Out[11]= {-Sqrt[1 - x^2 - y^2] <= z <= Sqrt[1 - x^2 - y^2], False}

In[12]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 <= 1},
{z}, {x, y, z}]

Out[12]= {(x == -1 &&
y == 0) || (-1 < x < 1 && -Sqrt[1 - x^2] <= y <= Sqrt[
1 - x^2]) || (x == 1 && y == 0), False}

In[13]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 <= 1},
{x, y, z}, {z}]

Out[13]= {-1 < x < 1 && -Sqrt[1 - x^2] < y < Sqrt[1 - x^2],
1 + x == 0 || -1 + x == 0 || -1 + x^2 + y^2 == 0}

のようにすれば，2つのリストに共通に属する変数の評価が出力されません．従って，xy 空間への射影（つまり，z の消去）は {x, y}, {x, y, z}，そこから x 空間への射影（つまり，z,y の消去）は {y, z}, {x, y, z} とすれば，それぞれ

In[14]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 <= 1}, {x,
y}, {x, y, z}]

Out[14]= {-Sqrt[1 - x^2 - y^2] <= z <= Sqrt[1 - x^2 - y^2], False}

In[15]:= GenericCylindricalDecomposition[{x^2 + y^2 + z^2 <= 1}, {y,
z}, {x, y, z}]

Out[15]= {-1 <= x <= 1, False}

のように取り出せます．
　なお，一連の関数の作者である Wolfram Research の Adam Strzebonski 氏が公開しておられる CAFView（http://members.wolfram.com/adams/Talks/Talk1007.nb）を用いると，CylindricalDecomposition の出力を

のように，セルの tree に対応した形で表示できます．