projection set に対する数値解の誤差（その２）

．．．とは言え，sep の評価なしでは，いくら R を小さくしても，前回の冒頭で述べたように，近接根に対応した数値解をカウントしてしまう可能性を排除出来ないかのように見えます．

しかし，我々が扱うのは一般の系ではなく，CAD の projection set であり，特に pscs を含む Collins-Hong タイプの projection（定義は，例えば，D.Wilson. Advances in Cylindrical Algebraic Decomposition の 2.3.1 Collins’Algorithm，2.4.1 Alternative Projection Operators をご参照ください）を利用すれば，sep の情報なしに，重根に対応した数値解全体のみを正確にカウントすることが出来ます．以下，そのことをお話しましょう．

まず，pscs（principal subresultant coefficient sequence）とは，多項式 f(x), g(x) に対して，それらの Sylvester 行列

syl(f1,f2,v):=block([d1:hipow(f1,v),d2:hipow(f2,v),c1,c2,z0,i,k],
           c1:makelist(coeff(f1,v,i),i,d1,0,-1),
           c2:makelist(coeff(f2,v,i),i,d2,0,-1),z0:makelist(0,d1+d2),
           [makelist(firstn(append(firstn(z0,k-1),c1,z0),d1+d2),k,d2),
            makelist(firstn(append(firstn(z0,k-1),c2,z0),d1+d2),k,d1)])$

の一部を除いた行列の行列式の列（普通の resultant はその初項です）

pscs(f01,f02,v):=block(
            [uratmx:ratmx,s1,s2,d,d1,d2,f1:expand(f01),f2:expand(f02),j],
            ratmx:on,[s1,s2]:syl(f1,f2,v),d1:hipow(f1,v),d2:hipow(f2,v),
            d:makelist(determinant(apply(matrix,
                                         append(
                                          map(lambda([e],firstn(e,d1+d2-2*j)),
                                              firstn(s1,d2-j)),
                                          map(lambda([e],firstn(e,d1+d2-2*j)),
                                              firstn(s2,d1-j))))),j,0,
                       min(d1,d2)-1),ratmx:uratmx,d)$

のことで

　( pscs の初項から続く 0 の個数 ) = ( gcd( f(x), g(x) ) の次数 )

という性質（Polynomial greatest common divisor - Wikipedia）を持ちます．例えば，上の（素朴な）コードを読み込んで，実行すると

(%i24) pscs( (x-1) * (x-2) , (x-3) * (x-4), x );
(%o24) [12,-4]

(%i25) pscs( (x-1) * (x-2)^3 * (x-3)^2 ,  (x-2)^2 * (x-3)^3 * (x-4), x );
(%o25) [0,0,0,0,12,-4]

(%i26) factor(pscs( (x-1) * (x-a)^3 * (x-b)^2 ,  (x-a)^2 * (x-b)^3 * (x-4), x ));
(%o26) [0,0,0,0,-3*(a-4)*(b-1)*(b-a),-(b-a+3)]

といった具合になります．

さて，上の性質を g(x) = f'(x) の場合に用いると

　gcd( f(x), f'(x) ) の次数

の値が pscs から特定できる訳ですが，この値が f(x) の

　(重根の個数) - (相異なる重根の個数)

であることは，積の微分公式により簡単に判ります．例で確認してみると

(%i31) pscs( f : (x-1) * (x-2)^3 * (x-3)^5 , diff(f,x) , x );
(%o31) [0,0,0,0,0,0,-60,-38]

のように 8 - 2 = 6 となっています．

ところが，この差は前回の冒頭で述べた方法における，重根に対応する数値解の検出カウントに他なりません．例えば，上の重根 2, 2, 2, 3, 3, 3, 3, 3 に対応した数値解 s1, s2, s3, t1, t2, t3, t4, t5 は
2*R >= | s1 - s2 |, 2*R >= | s2 - s3 |,
2*R >= | t1 - t2 |, 2*R >= | t2 - t3 |, 2*R >= | t3 - t4 |, 2*R >= | t4 - t5 |
として，合計 (3-1) + (5-1) = 6 回カウントされ，もし，検出カウントがこれを超えるなら，そこには近接根に対応する数値解のカウントが含まれているので，R を小さくして retry すれば，有限回で近接根に対応する数値解を排除でき，同時に，カウントの長さを見れば，各重根の重複度（各クラスタに属する数値解の個数）も判る，と言う仕組みです．