LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@kumanomi1000

陰関数でベジェ曲線の交点を求める ～ その7：実装へ反映

Posted at

重根への対応

分離区間がしきい値より小さい場合に2分割を反復する処理を打ち切るよう、VAG法を修正する必要があります。以下はVAG法で打ち切りが必要な例です。

inter_23.png $t=0.2, u=0.8$ で重複数2の重根

inter_22.png $t=0.2, u=0.2$ で重複数3の重根

二分法やRidder's法は重根に対応していないようなので、VAG法の結果の分離区間の中点を根とすることにします。本来、実根分離が目的のヴァンサンの定理／VAG法を、求根の収束計算に使うことになります。信頼性や計算効率の面でどうかと思われますが、そもそも重根は例外的なレアケースであり、さしあたって良しとします。

Wikipediaの『ヴァンサンの定理』のVAG法の疑似コードを、重根に対応するため以下の通り改造しました。（赤色の箇所）

1 $var \leftarrow (x + 1)^{deg(p)} p(\frac{a + bx}{1 + x})$ の符号変化の数
2 if $var = 0$ then RETURN $\varnothing$;
3 if $var = 1$ then RETURN $\set{(a, b)}$;
4 $m \leftarrow \frac{1}{2}(a + b)$
+ if $b-a \lt M$ then RETURN $\set{(m-A, m+A)}$;
5 if $p(m) \neq 0$ then
6   RETURN VAG($p$, $(a, m)$) $\cup$ VAG($p$, $(m, b)$)
7 else
+   $var \leftarrow (x + 1)^{deg(p)} p(\frac{m-A + (m+A)x}{1 + x})$ の符号変化の数
+   if $var = 1$ then
8     RETURN VAG($p$, $(a, m)$) $\cup$ $\set{[m, m]}$ $\cup$ VAG($p$, $(m, b)$)
+   else
+     RETURN VAG($p$, $(a, m-A)$) $\cup$ $\set{[m-A, m+A]}$ $\cup$ VAG($p$, $(m+A, b)$)
+   end
9 end

  • $M$は分離区間長の最小値です。区間がこれより短い場合、根が2個以上あっても分離区間として出力し、処理を打ち切ります。

  • $A$ は重根の場合の分離区間長の$\frac{1}{2}$です。

このように2つのパラメータを用いるようにした理由は、分離区間長の確保と中点の精度を両立する必要があるからです。
以下は『その5』のパターン④（カスプ点どうしが接している）の重複数4の根 $u=0.75$ に関して、関数 $(x + 1)^{deg(p)} p(\frac{a + bx}{1 + x})$ の各分離区間の係数と符号変化の数をSage Mathを用いて求めたものです。期待通りの符号変化の数を得るためには分離区間を長め（$A$を$10^{-3}$以上）に設定する必要があります。係数の値に鑑みて、浮動小数演算における計算誤差の影響によるものと推測されますが、解析が必要です。
ちなみに重複数3や2の他のパターンでは、より小さな$A$を設定しても期待通りの符号変化の数を得ることができます。重複数4が最悪ケースだとすると、カスプが排除できれば$A$をより小さな値に設定できるかもしれません。

$A$ $(x + 1)^{deg(p)} p(\frac{a + bx}{1 + x})$
の係数と符号変化の数
$10^{-6}$ VAG_10_1e-6.png
$10^{-5}$ VAG_10_1e-5.png
$10^{-4}$ VAG_10_1e-4.png
$10^{-3}$ VAG_10_1e-3.png

なお『実装への反映例』で$M$はABMIN、$A$はABADJです。

2次以下への対応

Sage Mathを用いて、曲線0が2次や1次の場合の陰伏方程式$F_0(u)$の係数を求める式を求めます。（$F_1(t)$は同様なので省略）

曲線0が2次の場合の$F_0(u)$の係数を求める式を求める
var('u a1 a2 a3 b1 b2 b3 c0 c1 c2 c3 d0 d1 d2 d3')
x=c0*u^3+c1*u^2+c2*u+c3
y=d0*u^3+d1*u^2+d2*u+d3
F0=a3^2*b1^2-a2*a3*b1*b2+a1*a3*b2^2+a2^2*b1*b3-2*a1*a3*b1*b3-a1*a2*b2*b3+a1^2*b3^2+(-2*a3*b1^2+a2*b1*b2-a1*b2^2+2*a1*b1*b3)*x+(-a2^2*b1+2*a1*a3*b1+a1*a2*b2-2*a1^2*b3)*y+b1^2*x^2-2*a1*b1*x*y+a1^2*y^2
cf0=F0.coefficients(u,sparse=False)
cf0[6].simplify_full()
cf0[5].simplify_full()
cf0[4].simplify_full()
cf0[3].simplify_full()
cf0[2].simplify_full()
cf0[1].simplify_full()
cf0[0].simplify_full()
実行結果
b1^2*c0^2 - 2*a1*b1*c0*d0 + a1^2*d0^2
2*b1^2*c0*c1 - 2*a1*b1*c1*d0 - 2*(a1*b1*c0 - a1^2*d0)*d1
b1^2*c1^2 + 2*b1^2*c0*c2 - 2*a1*b1*c2*d0 - 2*a1*b1*c1*d1 + a1^2*d1^2 - 2*(a1*b1*c0 - a1^2*d0)*d2
2*b1^2*c1*c2 + 2*b1^2*c0*c3 - 2*a1*b1*c2*d1 - (2*a3*b1^2 - a2*b1*b2 + a1*b2^2 - 2*a1*b1*b3)*c0 + (a1*a2*b2 - 2*a1^2*b3 - 2*a1*b1*c3 - (a2^2 - 2*a1*a3)*b1)*d0 - 2*(a1*b1*c1 - a1^2*d1)*d2 - 2*(a1*b1*c0 - a1^2*d0)*d3
b1^2*c2^2 + 2*b1^2*c1*c3 - 2*a1*b1*c2*d2 + a1^2*d2^2 - (2*a3*b1^2 - a2*b1*b2 + a1*b2^2 - 2*a1*b1*b3)*c1 + (a1*a2*b2 - 2*a1^2*b3 - 2*a1*b1*c3 - (a2^2 - 2*a1*a3)*b1)*d1 - 2*(a1*b1*c1 - a1^2*d1)*d3
2*b1^2*c2*c3 - (2*a3*b1^2 - a2*b1*b2 + a1*b2^2 - 2*a1*b1*b3)*c2 + (a1*a2*b2 - 2*a1^2*b3 - 2*a1*b1*c3 - (a2^2 - 2*a1*a3)*b1)*d2 - 2*(a1*b1*c2 - a1^2*d2)*d3
a3^2*b1^2 - a2*a3*b1*b2 + a1*a3*b2^2 + a1^2*b3^2 + b1^2*c3^2 + a1^2*d3^2 - (a1*a2*b2 - (a2^2 - 2*a1*a3)*b1)*b3 - (2*a3*b1^2 - a2*b1*b2 + a1*b2^2 - 2*a1*b1*b3)*c3 + (a1*a2*b2 - 2*a1^2*b3 - 2*a1*b1*c3 - (a2^2 - 2*a1*a3)*b1)*d3
曲線0が1次の場合の$F_0(u)$の係数を求める式を求める
var('u a2 a3 b2 b3 c0 c1 c2 c3 d0 d1 d2 d3')
x=c0*u^3+c1*u^2+c2*u+c3
y=d0*u^3+d1*u^2+d2*u+d3
F0=-a3*b2+a2*b3+b2*x-a2*y
cf0=F0.coefficients(u,sparse=False)
cf0[3].simplify_full()
cf0[2].simplify_full()
cf0[1].simplify_full()
cf0[0].simplify_full()
実行結果
b2*c0 - a2*d0
b2*c1 - a2*d1
b2*c2 - a2*d2
-a3*b2 + a2*b3 + b2*c3 - a2*d3

座標値の正規化

座標値を10倍や100倍にすると計算誤差が増大する傾向が見られたため、座標値を0～1の範囲になるよう正規化してから交点の根を求めるようにしました。
正規化はアスペクト比が維持される（xとyの倍率が同じになる）ようにしました。曲線が水平または垂直に細長い場合に、xまたはyの倍率が極端な値にならないようにするためです。

実装への反映例（『その4』からの更新ファイルのみ）

bzintersection.cpp
#include <cmath>
#include <QPointF>
#include <QList>
using namespace std;

#define EZERO   1e-14   // これより0に近い値は0とする
#define ABADJ   1e-3    // VAG法 重根の分離区間長/2
#define ABMIN   1e-5    // VAG法 分離区間長の最小
#define DSTMIN  1e-4    // 交点候補の突き合わせ判定の距離しきい値

#define PW3(x)	((x)*(x)*(x))
#define PW2(x)	((x)*(x))

typedef struct{
    double a;	// 区間from
    double b;	// 区間to
    double pa;	// aの関数値
    double pb;	// bの関数値
    int nr;     // 区間に含まれる実根の数
}Interval;

// 3次ベジェ曲線の係数を求める
void coefBezier(
  QList<QPointF> &cp,	// IN  ベジェ曲線 制御点
  double *a,			// OUT x係数
  double *b)			// OUT y係数
{
    a[0]=-cp[0].x()+3*cp[1].x()-3*cp[2].x()+cp[3].x();
    a[1]=3*cp[0].x()-6*cp[1].x()+3*cp[2].x();
    a[2]=-3*cp[0].x()+3*cp[1].x();
    a[3]=cp[0].x();
    b[0]=-cp[0].y()+3*cp[1].y()-3*cp[2].y()+cp[3].y();
    b[1]=3*cp[0].y()-6*cp[1].y()+3*cp[2].y();
    b[2]=-3*cp[0].y()+3*cp[1].y();
    b[3]=cp[0].y();
}

// 3次パラメトリック曲線の座標を求める
QPointF pointBezier(
  double *a,			// IN  x係数
  double *b,			// IN  y係数
  double t)				// IN  パラメータ
{
    return QPointF(a[0]*PW3(t)+a[1]*PW2(t)+a[2]*t+a[3],
                   b[0]*PW3(t)+b[1]*PW2(t)+b[2]*t+b[3]);
}

// 制御点の座標値を[0,1]に正規化
void normCP(
  QList<QPointF> &cp0,  // IN  曲線0の制御点
  QList<QPointF> &cp1,  // IN  曲線1の制御点
  QList<QPointF> &ncp0, // OUT 曲線0の制御点（正規化）
  QList<QPointF> &ncp1) // OUT 曲線0の制御点（正規化）
{
    double xmin=cp0[0].x();
    double xmax=cp0[0].x();
    double ymin=cp0[0].y();
    double ymax=cp0[0].y();
    foreach(QPointF p,cp0.mid(1)){
        if(xmin>p.x()) xmin=p.x();
        if(xmax<p.x()) xmax=p.x();
        if(ymin>p.y()) ymin=p.y();
        if(ymax<p.y()) ymax=p.y();
    }
    foreach(QPointF p,cp1){
        if(xmin>p.x()) xmin=p.x();
        if(xmax<p.x()) xmax=p.x();
        if(ymin>p.y()) ymin=p.y();
        if(ymax<p.y()) ymax=p.y();
    }
    double m=max(xmax-xmin,ymax-ymin);
    ncp0.clear();
    ncp1.clear();
    foreach(QPointF p,cp0)
        ncp0<<QPointF((p.x()-xmin)/m,(p.y()-ymin)/m);
    foreach(QPointF p,cp1)
        ncp1<<QPointF((p.x()-xmin)/m,(p.y()-ymin)/m);
}

// 3次パラメトリック曲線どうしの交点の陰伏方程式F0(u)の係数を求める
void coefImp0(
  double *a,	// IN  曲線0 (x,y)=f(t)の係数列 a[0],b[0]が3次項の係数
  double *b,
  double *c,	// IN  曲線1 (x,y)=g(u)の係数列 c[0],d[0]が3次項の係数
  double *d,
  double *cf)	// OUT F0(u)の係数列  cf[0]が9次項の係数
{
    if(abs(a[0])>=EZERO||abs(b[0])>=EZERO){
        // 曲線0は3次
        cf[0]=-PW3(b[0])*PW3(c[0])+3*a[0]*PW2(b[0])*PW2(c[0])*d[0]-3*PW2(a[0])*b[0]*c[0]*PW2(d[0])+PW3(a[0])*PW3(d[0]);
        cf[1]=-3*PW3(b[0])*PW2(c[0])*c[1]+6*a[0]*PW2(b[0])*c[0]*c[1]*d[0]-3*PW2(a[0])*b[0]*c[1]*PW2(d[0])+3*(a[0]*PW2(b[0])*PW2(c[0])-2*PW2(a[0])*b[0]*c[0]*d[0]+PW3(a[0])*PW2(d[0]))*d[1];
        cf[2]=-3*PW3(b[0])*c[0]*PW2(c[1])-3*PW3(b[0])*PW2(c[0])*c[2]-3*PW2(a[0])*b[0]*c[2]*PW2(d[0])-3*(PW2(a[0])*b[0]*c[0]-PW3(a[0])*d[0])*PW2(d[1])+3*(a[0]*PW2(b[0])*PW2(c[1])+2*a[0]*PW2(b[0])*c[0]*c[2])*d[0]+6*(a[0]*PW2(b[0])*c[0]*c[1]-PW2(a[0])*b[0]*c[1]*d[0])*d[1]+3*(a[0]*PW2(b[0])*PW2(c[0])-2*PW2(a[0])*b[0]*c[0]*d[0]+PW3(a[0])*PW2(d[0]))*d[2];
        cf[3]=-PW3(b[0])*PW3(c[1])-6*PW3(b[0])*c[0]*c[1]*c[2]-3*PW3(b[0])*PW2(c[0])*c[3]-3*PW2(a[0])*b[0]*c[1]*PW2(d[1])+PW3(a[0])*PW3(d[1])+(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*PW2(c[0])+(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*PW2(d[0])+(6*a[0]*PW2(b[0])*c[1]*c[2]+6*a[0]*PW2(b[0])*c[0]*c[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[0])*d[0]+3*(a[0]*PW2(b[0])*PW2(c[1])+2*a[0]*PW2(b[0])*c[0]*c[2]-2*PW2(a[0])*b[0]*c[2]*d[0])*d[1]+6*(a[0]*PW2(b[0])*c[0]*c[1]-PW2(a[0])*b[0]*c[1]*d[0]-(PW2(a[0])*b[0]*c[0]-PW3(a[0])*d[0])*d[1])*d[2]+3*(a[0]*PW2(b[0])*PW2(c[0])-2*PW2(a[0])*b[0]*c[0]*d[0]+PW3(a[0])*PW2(d[0]))*d[3];
        cf[4]=-3*PW3(b[0])*PW2(c[1])*c[2]-3*PW3(b[0])*c[0]*PW2(c[2])-6*PW3(b[0])*c[0]*c[1]*c[3]-3*PW2(a[0])*b[0]*c[2]*PW2(d[1])+2*(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*c[0]*c[1]-3*(PW2(a[0])*b[0]*c[0]-PW3(a[0])*d[0])*PW2(d[2])+(3*a[0]*PW2(b[0])*PW2(c[2])+6*a[0]*PW2(b[0])*c[1]*c[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[1])*d[0]+(6*a[0]*PW2(b[0])*c[1]*c[2]+6*a[0]*PW2(b[0])*c[0]*c[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[0]+2*(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*d[0])*d[1]+3*(a[0]*PW2(b[0])*PW2(c[1])+2*a[0]*PW2(b[0])*c[0]*c[2]-2*PW2(a[0])*b[0]*c[2]*d[0]-2*PW2(a[0])*b[0]*c[1]*d[1]+PW3(a[0])*PW2(d[1]))*d[2]+6*(a[0]*PW2(b[0])*c[0]*c[1]-PW2(a[0])*b[0]*c[1]*d[0]-(PW2(a[0])*b[0]*c[0]-PW3(a[0])*d[0])*d[1])*d[3];
        cf[5]=-3*PW3(b[0])*c[1]*PW2(c[2])+(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*PW2(c[1])+2*(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*c[0]*c[2]+(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*PW2(d[1])-3*(PW2(a[0])*b[0]*c[1]-PW3(a[0])*d[1])*PW2(d[2])-3*(PW3(b[0])*PW2(c[1])+2*PW3(b[0])*c[0]*c[2])*c[3]+(6*a[0]*PW2(b[0])*c[2]*c[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[2])*d[0]+(3*a[0]*PW2(b[0])*PW2(c[2])+6*a[0]*PW2(b[0])*c[1]*c[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[1])*d[1]+(6*a[0]*PW2(b[0])*c[1]*c[2]+6*a[0]*PW2(b[0])*c[0]*c[3]-6*PW2(a[0])*b[0]*c[2]*d[1]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[0]+2*(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*d[0])*d[2]+3*(a[0]*PW2(b[0])*PW2(c[1])+2*a[0]*PW2(b[0])*c[0]*c[2]-2*PW2(a[0])*b[0]*c[2]*d[0]-2*PW2(a[0])*b[0]*c[1]*d[1]+PW3(a[0])*PW2(d[1])-2*(PW2(a[0])*b[0]*c[0]-PW3(a[0])*d[0])*d[2])*d[3];
        cf[6]=-PW3(b[0])*PW3(c[2])-3*PW3(b[0])*c[0]*PW2(c[3])-3*PW2(a[0])*b[0]*c[2]*PW2(d[2])+PW3(a[0])*PW3(d[2])+2*(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*c[1]*c[2]-3*(PW2(a[0])*b[0]*c[0]-PW3(a[0])*d[0])*PW2(d[3])-(3*PW2(a[3])*PW3(b[0])-2*a[2]*a[3]*PW2(b[0])*b[1]+2*a[1]*a[3]*b[0]*PW2(b[1])-2*a[0]*a[3]*PW3(b[1])+PW2(a[0])*PW3(b[2])+3*PW2(a[0])*b[0]*PW2(b[3])-(a[0]*a[1]*b[1]-(PW2(a[1])-2*a[0]*a[2])*b[0])*PW2(b[2])+(a[0]*a[2]*PW2(b[1])+(PW2(a[2])-4*a[1]*a[3])*PW2(b[0])-(a[1]*a[2]-6*a[0]*a[3])*b[0]*b[1])*b[2]+(2*a[0]*a[1]*PW2(b[1])+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*b[3])*c[0]-2*(3*PW3(b[0])*c[1]*c[2]-(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*c[0])*c[3]+(PW2(a[0])*a[2]*PW2(b[2])+3*PW3(a[0])*PW2(b[3])+3*a[0]*PW2(b[0])*PW2(c[3])+(PW3(a[2])-3*a[1]*a[2]*a[3]+3*a[0]*PW2(a[3]))*PW2(b[0])-(a[1]*PW2(a[2])-(2*PW2(a[1])+a[0]*a[2])*a[3])*b[0]*b[1]+(a[0]*PW2(a[2])-2*a[0]*a[1]*a[3])*PW2(b[1])+((PW2(a[1])*a[2]-2*a[0]*PW2(a[2])-a[0]*a[1]*a[3])*b[0]-(a[0]*a[1]*a[2]-3*PW2(a[0])*a[3])*b[1])*b[2]-2*(PW2(a[0])*a[1]*b[2]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*b[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[3])*d[0]+(6*a[0]*PW2(b[0])*c[2]*c[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[2])*d[1]+(3*a[0]*PW2(b[0])*PW2(c[2])+6*a[0]*PW2(b[0])*c[1]*c[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[1]+2*(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*d[1])*d[2]+(6*a[0]*PW2(b[0])*c[1]*c[2]+6*a[0]*PW2(b[0])*c[0]*c[3]-6*PW2(a[0])*b[0]*c[2]*d[1]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[0]+2*(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*d[0]-6*(PW2(a[0])*b[0]*c[1]-PW3(a[0])*d[1])*d[2])*d[3];
        cf[7]=-3*PW3(b[0])*c[1]*PW2(c[3])+(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*PW2(c[2])+(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*PW2(d[2])-3*(PW2(a[0])*b[0]*c[1]-PW3(a[0])*d[1])*PW2(d[3])-(3*PW2(a[3])*PW3(b[0])-2*a[2]*a[3]*PW2(b[0])*b[1]+2*a[1]*a[3]*b[0]*PW2(b[1])-2*a[0]*a[3]*PW3(b[1])+PW2(a[0])*PW3(b[2])+3*PW2(a[0])*b[0]*PW2(b[3])-(a[0]*a[1]*b[1]-(PW2(a[1])-2*a[0]*a[2])*b[0])*PW2(b[2])+(a[0]*a[2]*PW2(b[1])+(PW2(a[2])-4*a[1]*a[3])*PW2(b[0])-(a[1]*a[2]-6*a[0]*a[3])*b[0]*b[1])*b[2]+(2*a[0]*a[1]*PW2(b[1])+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*b[3])*c[1]-(3*PW3(b[0])*PW2(c[2])-2*(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*c[1])*c[3]+(PW2(a[0])*a[2]*PW2(b[2])+3*PW3(a[0])*PW2(b[3])+3*a[0]*PW2(b[0])*PW2(c[3])+(PW3(a[2])-3*a[1]*a[2]*a[3]+3*a[0]*PW2(a[3]))*PW2(b[0])-(a[1]*PW2(a[2])-(2*PW2(a[1])+a[0]*a[2])*a[3])*b[0]*b[1]+(a[0]*PW2(a[2])-2*a[0]*a[1]*a[3])*PW2(b[1])+((PW2(a[1])*a[2]-2*a[0]*PW2(a[2])-a[0]*a[1]*a[3])*b[0]-(a[0]*a[1]*a[2]-3*PW2(a[0])*a[3])*b[1])*b[2]-2*(PW2(a[0])*a[1]*b[2]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*b[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[3])*d[1]+(6*a[0]*PW2(b[0])*c[2]*c[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[2])*d[2]+(3*a[0]*PW2(b[0])*PW2(c[2])+6*a[0]*PW2(b[0])*c[1]*c[3]-6*PW2(a[0])*b[0]*c[2]*d[2]+3*PW3(a[0])*PW2(d[2])+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[1]+2*(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*d[1])*d[3];
        cf[8]=-3*PW3(b[0])*c[2]*PW2(c[3])+2*(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*c[2]*c[3]-3*(PW2(a[0])*b[0]*c[2]-PW3(a[0])*d[2])*PW2(d[3])-(3*PW2(a[3])*PW3(b[0])-2*a[2]*a[3]*PW2(b[0])*b[1]+2*a[1]*a[3]*b[0]*PW2(b[1])-2*a[0]*a[3]*PW3(b[1])+PW2(a[0])*PW3(b[2])+3*PW2(a[0])*b[0]*PW2(b[3])-(a[0]*a[1]*b[1]-(PW2(a[1])-2*a[0]*a[2])*b[0])*PW2(b[2])+(a[0]*a[2]*PW2(b[1])+(PW2(a[2])-4*a[1]*a[3])*PW2(b[0])-(a[1]*a[2]-6*a[0]*a[3])*b[0]*b[1])*b[2]+(2*a[0]*a[1]*PW2(b[1])+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*b[3])*c[2]+(PW2(a[0])*a[2]*PW2(b[2])+3*PW3(a[0])*PW2(b[3])+3*a[0]*PW2(b[0])*PW2(c[3])+(PW3(a[2])-3*a[1]*a[2]*a[3]+3*a[0]*PW2(a[3]))*PW2(b[0])-(a[1]*PW2(a[2])-(2*PW2(a[1])+a[0]*a[2])*a[3])*b[0]*b[1]+(a[0]*PW2(a[2])-2*a[0]*a[1]*a[3])*PW2(b[1])+((PW2(a[1])*a[2]-2*a[0]*PW2(a[2])-a[0]*a[1]*a[3])*b[0]-(a[0]*a[1]*a[2]-3*PW2(a[0])*a[3])*b[1])*b[2]-2*(PW2(a[0])*a[1]*b[2]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*b[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[3])*d[2]+(6*a[0]*PW2(b[0])*c[2]*c[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[2]+2*(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*d[2])*d[3];
        cf[9]=PW3(a[3])*PW3(b[0])-a[2]*PW2(a[3])*PW2(b[0])*b[1]+a[1]*PW2(a[3])*b[0]*PW2(b[1])-a[0]*PW2(a[3])*PW3(b[1])+PW2(a[0])*a[3]*PW3(b[2])-PW3(a[0])*PW3(b[3])-PW3(b[0])*PW3(c[3])+PW3(a[0])*PW3(d[3])-(a[0]*a[1]*a[3]*b[1]-(PW2(a[1])-2*a[0]*a[2])*a[3]*b[0])*PW2(b[2])+(PW2(a[0])*a[1]*b[2]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*PW2(b[3])+(3*a[3]*PW3(b[0])-a[2]*PW2(b[0])*b[1]+a[1]*b[0]*PW2(b[1])-a[0]*PW3(b[1])-3*a[0]*PW2(b[0])*b[3]-(2*a[1]*PW2(b[0])-3*a[0]*b[0]*b[1])*b[2])*PW2(c[3])+(PW2(a[0])*a[1]*b[2]-3*PW3(a[0])*b[3]-3*PW2(a[0])*b[0]*c[3]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*PW2(d[3])+(a[0]*a[2]*a[3]*PW2(b[1])+(PW2(a[2])*a[3]-2*a[1]*PW2(a[3]))*PW2(b[0])-(a[1]*a[2]*a[3]-3*a[0]*PW2(a[3]))*b[0]*b[1])*b[2]-(PW2(a[0])*a[2]*PW2(b[2])+(PW3(a[2])-3*a[1]*a[2]*a[3]+3*a[0]*PW2(a[3]))*PW2(b[0])-(a[1]*PW2(a[2])-(2*PW2(a[1])+a[0]*a[2])*a[3])*b[0]*b[1]+(a[0]*PW2(a[2])-2*a[0]*a[1]*a[3])*PW2(b[1])+((PW2(a[1])*a[2]-2*a[0]*PW2(a[2])-a[0]*a[1]*a[3])*b[0]-(a[0]*a[1]*a[2]-3*PW2(a[0])*a[3])*b[1])*b[2])*b[3]-(3*PW2(a[3])*PW3(b[0])-2*a[2]*a[3]*PW2(b[0])*b[1]+2*a[1]*a[3]*b[0]*PW2(b[1])-2*a[0]*a[3]*PW3(b[1])+PW2(a[0])*PW3(b[2])+3*PW2(a[0])*b[0]*PW2(b[3])-(a[0]*a[1]*b[1]-(PW2(a[1])-2*a[0]*a[2])*b[0])*PW2(b[2])+(a[0]*a[2]*PW2(b[1])+(PW2(a[2])-4*a[1]*a[3])*PW2(b[0])-(a[1]*a[2]-6*a[0]*a[3])*b[0]*b[1])*b[2]+(2*a[0]*a[1]*PW2(b[1])+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*b[3])*c[3]+(PW2(a[0])*a[2]*PW2(b[2])+3*PW3(a[0])*PW2(b[3])+3*a[0]*PW2(b[0])*PW2(c[3])+(PW3(a[2])-3*a[1]*a[2]*a[3]+3*a[0]*PW2(a[3]))*PW2(b[0])-(a[1]*PW2(a[2])-(2*PW2(a[1])+a[0]*a[2])*a[3])*b[0]*b[1]+(a[0]*PW2(a[2])-2*a[0]*a[1]*a[3])*PW2(b[1])+((PW2(a[1])*a[2]-2*a[0]*PW2(a[2])-a[0]*a[1]*a[3])*b[0]-(a[0]*a[1]*a[2]-3*PW2(a[0])*a[3])*b[1])*b[2]-2*(PW2(a[0])*a[1]*b[2]+(PW3(a[1])-3*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3])*b[0]-(a[0]*PW2(a[1])-2*PW2(a[0])*a[2])*b[1])*b[3]+(2*a[0]*a[1]*PW2(b[1])+6*PW2(a[0])*b[0]*b[3]+3*(a[1]*a[2]-2*a[0]*a[3])*PW2(b[0])-(2*PW2(a[1])+a[0]*a[2])*b[0]*b[1]+(a[0]*a[1]*b[0]-3*PW2(a[0])*b[1])*b[2])*c[3])*d[3];
    }else if(abs(a[1])>=EZERO||abs(b[1])>=EZERO){
        // 曲線0は2次
        cf[0]=0;
        cf[1]=0;
        cf[2]=0;
        cf[3]=PW2(b[1])*PW2(c[0])-2*a[1]*b[1]*c[0]*d[0]+PW2(a[1])*PW2(d[0]);
        cf[4]=2*PW2(b[1])*c[0]*c[1]-2*a[1]*b[1]*c[1]*d[0]-2*(a[1]*b[1]*c[0]-PW2(a[1])*d[0])*d[1];
        cf[5]=PW2(b[1])*PW2(c[1])+2*PW2(b[1])*c[0]*c[2]-2*a[1]*b[1]*c[2]*d[0]-2*a[1]*b[1]*c[1]*d[1]+PW2(a[1])*PW2(d[1])-2*(a[1]*b[1]*c[0]-PW2(a[1])*d[0])*d[2];
        cf[6]=2*PW2(b[1])*c[1]*c[2]+2*PW2(b[1])*c[0]*c[3]-2*a[1]*b[1]*c[2]*d[1]-(2*a[3]*PW2(b[1])-a[2]*b[1]*b[2]+a[1]*PW2(b[2])-2*a[1]*b[1]*b[3])*c[0]+(a[1]*a[2]*b[2]-2*PW2(a[1])*b[3]-2*a[1]*b[1]*c[3]-(PW2(a[2])-2*a[1]*a[3])*b[1])*d[0]-2*(a[1]*b[1]*c[1]-PW2(a[1])*d[1])*d[2]-2*(a[1]*b[1]*c[0]-PW2(a[1])*d[0])*d[3];
        cf[7]=PW2(b[1])*PW2(c[2])+2*PW2(b[1])*c[1]*c[3]-2*a[1]*b[1]*c[2]*d[2]+PW2(a[1])*PW2(d[2])-(2*a[3]*PW2(b[1])-a[2]*b[1]*b[2]+a[1]*PW2(b[2])-2*a[1]*b[1]*b[3])*c[1]+(a[1]*a[2]*b[2]-2*PW2(a[1])*b[3]-2*a[1]*b[1]*c[3]-(PW2(a[2])-2*a[1]*a[3])*b[1])*d[1]-2*(a[1]*b[1]*c[1]-PW2(a[1])*d[1])*d[3];
        cf[8]=2*PW2(b[1])*c[2]*c[3]-(2*a[3]*PW2(b[1])-a[2]*b[1]*b[2]+a[1]*PW2(b[2])-2*a[1]*b[1]*b[3])*c[2]+(a[1]*a[2]*b[2]-2*PW2(a[1])*b[3]-2*a[1]*b[1]*c[3]-(PW2(a[2])-2*a[1]*a[3])*b[1])*d[2]-2*(a[1]*b[1]*c[2]-PW2(a[1])*d[2])*d[3];
        cf[9]=PW2(a[3])*PW2(b[1])-a[2]*a[3]*b[1]*b[2]+a[1]*a[3]*PW2(b[2])+PW2(a[1])*PW2(b[3])+PW2(b[1])*PW2(c[3])+PW2(a[1])*PW2(d[3])-(a[1]*a[2]*b[2]-(PW2(a[2])-2*a[1]*a[3])*b[1])*b[3]-(2*a[3]*PW2(b[1])-a[2]*b[1]*b[2]+a[1]*PW2(b[2])-2*a[1]*b[1]*b[3])*c[3]+(a[1]*a[2]*b[2]-2*PW2(a[1])*b[3]-2*a[1]*b[1]*c[3]-(PW2(a[2])-2*a[1]*a[3])*b[1])*d[3];
    }else{
        // 曲線0は1次
        cf[0]=0;
        cf[1]=0;
        cf[2]=0;
        cf[3]=0;
        cf[4]=0;
        cf[5]=0;
        cf[6]=b[2]*c[0]-a[2]*d[0];
        cf[7]=b[2]*c[1]-a[2]*d[1];
        cf[8]=b[2]*c[2]-a[2]*d[2];
        cf[9]=-a[3]*b[2]+a[2]*b[3]+b[2]*c[3]-a[2]*d[3];
    }
}

// 3次パラメトリック曲線どうしの交点の陰伏方程式F1(t)の係数を求める
void coefImp1(
  double *a,	// IN  曲線0 (x,y)=f(t)の係数列 a[0],b[0]が3次項の係数
  double *b,
  double *c,	// IN  曲線1 (x,y)=g(u)の係数列 c[0],d[0]が3次項の係数
  double *d,
  double *cf)	// OUT F1(t)の係数列  cf[0]が9次項の係数
{
    if(abs(c[0])>=EZERO||abs(d[0])>=EZERO){
        // 曲線1は3次
        cf[0]=PW3(b[0])*PW3(c[0])-3*a[0]*PW2(b[0])*PW2(c[0])*d[0]+3*PW2(a[0])*b[0]*c[0]*PW2(d[0])-PW3(a[0])*PW3(d[0]);
        cf[1]=3*PW2(b[0])*b[1]*PW3(c[0])-3*PW2(a[0])*a[1]*PW3(d[0])-3*(a[1]*PW2(b[0])+2*a[0]*b[0]*b[1])*PW2(c[0])*d[0]+3*(2*a[0]*a[1]*b[0]+PW2(a[0])*b[1])*c[0]*PW2(d[0]);
        cf[2]=3*(b[0]*PW2(b[1])+PW2(b[0])*b[2])*PW3(c[0])-3*(a[2]*PW2(b[0])+2*a[1]*b[0]*b[1]+a[0]*PW2(b[1])+2*a[0]*b[0]*b[2])*PW2(c[0])*d[0]+3*(2*a[0]*a[1]*b[1]+PW2(a[0])*b[2]+(PW2(a[1])+2*a[0]*a[2])*b[0])*c[0]*PW2(d[0])-3*(a[0]*PW2(a[1])+PW2(a[0])*a[2])*PW3(d[0]);
        cf[3]=-PW2(a[0])*c[0]*PW3(d[1])+(PW3(b[1])+6*b[0]*b[1]*b[2]+3*PW2(b[0])*b[3])*PW3(c[0])-(PW3(a[1])+6*a[0]*a[1]*a[2]+3*PW2(a[0])*a[3]-3*PW2(a[0])*c[3])*PW3(d[0])+3*(a[0]*b[0]*c[1]*c[2]-2*a[0]*b[0]*c[0]*c[3]+(2*a[0]*a[1]*b[2]+PW2(a[0])*b[3]+2*(a[1]*a[2]+a[0]*a[3])*b[0]+(PW2(a[1])+2*a[0]*a[2])*b[1])*c[0])*PW2(d[0])+(2*a[0]*b[0]*c[0]*c[1]+PW2(a[0])*c[1]*d[0])*PW2(d[1])+(PW2(b[0])*PW3(c[1])-3*PW2(b[0])*c[0]*c[1]*c[2]+3*PW2(b[0])*PW2(c[0])*c[3]-3*(a[3]*PW2(b[0])+2*a[2]*b[0]*b[1]+a[1]*PW2(b[1])+2*a[0]*b[0]*b[3]+2*(a[1]*b[0]+a[0]*b[1])*b[2])*PW2(c[0]))*d[0]-(PW2(b[0])*c[0]*PW2(c[1])-2*PW2(b[0])*PW2(c[0])*c[2]+PW2(a[0])*c[2]*PW2(d[0])+(2*a[0]*b[0]*PW2(c[1])+a[0]*b[0]*c[0]*c[2])*d[0])*d[1]+(PW2(b[0])*PW2(c[0])*c[1]+a[0]*b[0]*c[0]*c[1]*d[0]-2*PW2(a[0])*c[1]*PW2(d[0])-3*(a[0]*b[0]*PW2(c[0])-PW2(a[0])*c[0]*d[0])*d[1])*d[2]-3*(PW2(b[0])*PW3(c[0])-2*a[0]*b[0]*PW2(c[0])*d[0]+PW2(a[0])*c[0]*PW2(d[0]))*d[3];
        cf[4]=-2*a[0]*a[1]*c[0]*PW3(d[1])+3*(PW2(b[1])*b[2]+b[0]*PW2(b[2])+2*b[0]*b[1]*b[3])*PW3(c[0])-3*(PW2(a[1])*a[2]+a[0]*PW2(a[2])+2*a[0]*a[1]*a[3]-2*a[0]*a[1]*c[3])*PW3(d[0])+3*((a[1]*b[0]+a[0]*b[1])*c[1]*c[2]-2*(a[1]*b[0]+a[0]*b[1])*c[0]*c[3]+(2*a[0]*a[1]*b[3]+(PW2(a[2])+2*a[1]*a[3])*b[0]+2*(a[1]*a[2]+a[0]*a[3])*b[1]+(PW2(a[1])+2*a[0]*a[2])*b[2])*c[0])*PW2(d[0])+2*(a[0]*a[1]*c[1]*d[0]+(a[1]*b[0]+a[0]*b[1])*c[0]*c[1])*PW2(d[1])+(2*b[0]*b[1]*PW3(c[1])-6*b[0]*b[1]*c[0]*c[1]*c[2]+6*b[0]*b[1]*PW2(c[0])*c[3]-3*(2*a[3]*b[0]*b[1]+a[2]*PW2(b[1])+a[0]*PW2(b[2])+2*(a[2]*b[0]+a[1]*b[1])*b[2]+2*(a[1]*b[0]+a[0]*b[1])*b[3])*PW2(c[0]))*d[0]-(2*b[0]*b[1]*c[0]*PW2(c[1])-4*b[0]*b[1]*PW2(c[0])*c[2]+2*a[0]*a[1]*c[2]*PW2(d[0])+(2*(a[1]*b[0]+a[0]*b[1])*PW2(c[1])+(a[1]*b[0]+a[0]*b[1])*c[0]*c[2])*d[0])*d[1]+(2*b[0]*b[1]*PW2(c[0])*c[1]-4*a[0]*a[1]*c[1]*PW2(d[0])+(a[1]*b[0]+a[0]*b[1])*c[0]*c[1]*d[0]+3*(2*a[0]*a[1]*c[0]*d[0]-(a[1]*b[0]+a[0]*b[1])*PW2(c[0]))*d[1])*d[2]-6*(b[0]*b[1]*PW3(c[0])+a[0]*a[1]*c[0]*PW2(d[0])-(a[1]*b[0]+a[0]*b[1])*PW2(c[0])*d[0])*d[3];
        cf[5]=-(PW2(a[1])+2*a[0]*a[2])*c[0]*PW3(d[1])+3*(b[1]*PW2(b[2])+(PW2(b[1])+2*b[0]*b[2])*b[3])*PW3(c[0])-3*(a[1]*PW2(a[2])+(PW2(a[1])+2*a[0]*a[2])*a[3]-(PW2(a[1])+2*a[0]*a[2])*c[3])*PW3(d[0])+3*((a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*c[1]*c[2]-2*(a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*c[0]*c[3]+(2*a[2]*a[3]*b[0]+(PW2(a[2])+2*a[1]*a[3])*b[1]+2*(a[1]*a[2]+a[0]*a[3])*b[2]+(PW2(a[1])+2*a[0]*a[2])*b[3])*c[0])*PW2(d[0])+(2*(a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*c[0]*c[1]+(PW2(a[1])+2*a[0]*a[2])*c[1]*d[0])*PW2(d[1])+((PW2(b[1])+2*b[0]*b[2])*PW3(c[1])-3*(PW2(b[1])+2*b[0]*b[2])*c[0]*c[1]*c[2]+3*(PW2(b[1])+2*b[0]*b[2])*PW2(c[0])*c[3]-3*(a[3]*PW2(b[1])+a[1]*PW2(b[2])+2*(a[3]*b[0]+a[2]*b[1])*b[2]+2*(a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*b[3])*PW2(c[0]))*d[0]-((PW2(b[1])+2*b[0]*b[2])*c[0]*PW2(c[1])-2*(PW2(b[1])+2*b[0]*b[2])*PW2(c[0])*c[2]+(PW2(a[1])+2*a[0]*a[2])*c[2]*PW2(d[0])+(2*(a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*PW2(c[1])+(a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*c[0]*c[2])*d[0])*d[1]+((PW2(b[1])+2*b[0]*b[2])*PW2(c[0])*c[1]+(a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*c[0]*c[1]*d[0]-2*(PW2(a[1])+2*a[0]*a[2])*c[1]*PW2(d[0])-3*((a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*PW2(c[0])-(PW2(a[1])+2*a[0]*a[2])*c[0]*d[0])*d[1])*d[2]-3*((PW2(b[1])+2*b[0]*b[2])*PW3(c[0])-2*(a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*PW2(c[0])*d[0]+(PW2(a[1])+2*a[0]*a[2])*c[0]*PW2(d[0]))*d[3];
        cf[6]=-a[0]*PW2(c[0])*PW3(d[2])+(PW3(b[2])+6*b[1]*b[2]*b[3]+3*b[0]*PW2(b[3]))*PW3(c[0])-(PW3(a[2])+6*a[1]*a[2]*a[3]+3*a[0]*PW2(a[3])+3*a[0]*PW2(c[3])-6*(a[1]*a[2]+a[0]*a[3])*c[3])*PW3(d[0])+2*(a[0]*c[0]*c[3]-(a[1]*a[2]+a[0]*a[3])*c[0])*PW3(d[1])+(b[0]*PW3(c[2])+3*b[0]*c[0]*PW2(c[3])+3*(a[3]*b[0]+a[2]*b[1]+a[1]*b[2]+a[0]*b[3])*c[1]*c[2]+3*(PW2(a[3])*b[0]+2*a[2]*a[3]*b[1]+(PW2(a[2])+2*a[1]*a[3])*b[2]+2*(a[1]*a[2]+a[0]*a[3])*b[3])*c[0]-3*(b[0]*c[1]*c[2]+2*(a[3]*b[0]+a[2]*b[1]+a[1]*b[2]+a[0]*b[3])*c[0])*c[3])*PW2(d[0])+(b[0]*c[0]*PW2(c[2])-2*b[0]*c[0]*c[1]*c[3]+2*(a[3]*b[0]+a[2]*b[1]+a[1]*b[2]+a[0]*b[3])*c[0]*c[1]-2*(a[0]*c[1]*c[3]-(a[1]*a[2]+a[0]*a[3])*c[1])*d[0])*PW2(d[1])+(b[0]*PW2(c[0])*c[2]+a[0]*c[0]*c[1]*d[1]-(a[0]*PW2(c[1])-2*a[0]*c[0]*c[2])*d[0])*PW2(d[2])+3*(b[0]*PW3(c[0])-a[0]*PW2(c[0])*d[0])*PW2(d[3])+(2*(b[1]*b[2]+b[0]*b[3])*PW3(c[1])-6*(b[1]*b[2]+b[0]*b[3])*c[0]*c[1]*c[2]+6*(b[1]*b[2]+b[0]*b[3])*PW2(c[0])*c[3]-3*(2*a[3]*b[1]*b[2]+a[2]*PW2(b[2])+a[0]*PW2(b[3])+2*(a[3]*b[0]+a[2]*b[1]+a[1]*b[2])*b[3])*PW2(c[0]))*d[0]-(2*(b[1]*b[2]+b[0]*b[3])*c[0]*PW2(c[1])-4*(b[1]*b[2]+b[0]*b[3])*PW2(c[0])*c[2]-2*(a[0]*c[2]*c[3]-(a[1]*a[2]+a[0]*a[3])*c[2])*PW2(d[0])+(b[0]*c[1]*PW2(c[2])+2*(a[3]*b[0]+a[2]*b[1]+a[1]*b[2]+a[0]*b[3])*PW2(c[1])+(a[3]*b[0]+a[2]*b[1]+a[1]*b[2]+a[0]*b[3])*c[0]*c[2]-(2*b[0]*PW2(c[1])+b[0]*c[0]*c[2])*c[3])*d[0])*d[1]-(a[0]*c[0]*c[2]*PW2(d[1])-2*(b[1]*b[2]+b[0]*b[3])*PW2(c[0])*c[1]+(a[0]*PW2(c[2])-4*a[0]*c[1]*c[3]+4*(a[1]*a[2]+a[0]*a[3])*c[1])*PW2(d[0])-(b[0]*PW2(c[1])*c[2]-2*b[0]*c[0]*PW2(c[2])-b[0]*c[0]*c[1]*c[3]+(a[3]*b[0]+a[2]*b[1]+a[1]*b[2]+a[0]*b[3])*c[0]*c[1])*d[0]+(b[0]*c[0]*c[1]*c[2]-3*b[0]*PW2(c[0])*c[3]+3*(a[3]*b[0]+a[2]*b[1]+a[1]*b[2]+a[0]*b[3])*PW2(c[0])-(a[0]*c[1]*c[2]-6*a[0]*c[0]*c[3]+6*(a[1]*a[2]+a[0]*a[3])*c[0])*d[0])*d[1])*d[2]-(2*a[0]*c[0]*c[1]*PW2(d[1])+6*(b[1]*b[2]+b[0]*b[3])*PW3(c[0])+3*(a[0]*c[1]*c[2]-2*a[0]*c[0]*c[3]+2*(a[1]*a[2]+a[0]*a[3])*c[0])*PW2(d[0])+2*(b[0]*PW3(c[1])-3*b[0]*c[0]*c[1]*c[2]+3*b[0]*PW2(c[0])*c[3]-3*(a[3]*b[0]+a[2]*b[1]+a[1]*b[2]+a[0]*b[3])*PW2(c[0]))*d[0]-(2*b[0]*c[0]*PW2(c[1])-4*b[0]*PW2(c[0])*c[2]+(2*a[0]*PW2(c[1])+a[0]*c[0]*c[2])*d[0])*d[1]+(2*b[0]*PW2(c[0])*c[1]+a[0]*c[0]*c[1]*d[0]-3*a[0]*PW2(c[0])*d[1])*d[2])*d[3];
        cf[7]=-a[1]*PW2(c[0])*PW3(d[2])+3*(PW2(b[2])*b[3]+b[1]*PW2(b[3]))*PW3(c[0])-3*(PW2(a[2])*a[3]+a[1]*PW2(a[3])+a[1]*PW2(c[3])-(PW2(a[2])+2*a[1]*a[3])*c[3])*PW3(d[0])+(2*a[1]*c[0]*c[3]-(PW2(a[2])+2*a[1]*a[3])*c[0])*PW3(d[1])+(b[1]*PW3(c[2])+3*b[1]*c[0]*PW2(c[3])+3*(a[3]*b[1]+a[2]*b[2]+a[1]*b[3])*c[1]*c[2]+3*(PW2(a[3])*b[1]+2*a[2]*a[3]*b[2]+(PW2(a[2])+2*a[1]*a[3])*b[3])*c[0]-3*(b[1]*c[1]*c[2]+2*(a[3]*b[1]+a[2]*b[2]+a[1]*b[3])*c[0])*c[3])*PW2(d[0])+(b[1]*c[0]*PW2(c[2])-2*b[1]*c[0]*c[1]*c[3]+2*(a[3]*b[1]+a[2]*b[2]+a[1]*b[3])*c[0]*c[1]-(2*a[1]*c[1]*c[3]-(PW2(a[2])+2*a[1]*a[3])*c[1])*d[0])*PW2(d[1])+(b[1]*PW2(c[0])*c[2]+a[1]*c[0]*c[1]*d[1]-(a[1]*PW2(c[1])-2*a[1]*c[0]*c[2])*d[0])*PW2(d[2])+3*(b[1]*PW3(c[0])-a[1]*PW2(c[0])*d[0])*PW2(d[3])+((PW2(b[2])+2*b[1]*b[3])*PW3(c[1])-3*(PW2(b[2])+2*b[1]*b[3])*c[0]*c[1]*c[2]+3*(PW2(b[2])+2*b[1]*b[3])*PW2(c[0])*c[3]-3*(a[3]*PW2(b[2])+a[1]*PW2(b[3])+2*(a[3]*b[1]+a[2]*b[2])*b[3])*PW2(c[0]))*d[0]-((PW2(b[2])+2*b[1]*b[3])*c[0]*PW2(c[1])-2*(PW2(b[2])+2*b[1]*b[3])*PW2(c[0])*c[2]-(2*a[1]*c[2]*c[3]-(PW2(a[2])+2*a[1]*a[3])*c[2])*PW2(d[0])+(b[1]*c[1]*PW2(c[2])+2*(a[3]*b[1]+a[2]*b[2]+a[1]*b[3])*PW2(c[1])+(a[3]*b[1]+a[2]*b[2]+a[1]*b[3])*c[0]*c[2]-(2*b[1]*PW2(c[1])+b[1]*c[0]*c[2])*c[3])*d[0])*d[1]-(a[1]*c[0]*c[2]*PW2(d[1])-(PW2(b[2])+2*b[1]*b[3])*PW2(c[0])*c[1]+(a[1]*PW2(c[2])-4*a[1]*c[1]*c[3]+2*(PW2(a[2])+2*a[1]*a[3])*c[1])*PW2(d[0])-(b[1]*PW2(c[1])*c[2]-2*b[1]*c[0]*PW2(c[2])-b[1]*c[0]*c[1]*c[3]+(a[3]*b[1]+a[2]*b[2]+a[1]*b[3])*c[0]*c[1])*d[0]+(b[1]*c[0]*c[1]*c[2]-3*b[1]*PW2(c[0])*c[3]+3*(a[3]*b[1]+a[2]*b[2]+a[1]*b[3])*PW2(c[0])-(a[1]*c[1]*c[2]-6*a[1]*c[0]*c[3]+3*(PW2(a[2])+2*a[1]*a[3])*c[0])*d[0])*d[1])*d[2]-(2*a[1]*c[0]*c[1]*PW2(d[1])+3*(PW2(b[2])+2*b[1]*b[3])*PW3(c[0])+3*(a[1]*c[1]*c[2]-2*a[1]*c[0]*c[3]+(PW2(a[2])+2*a[1]*a[3])*c[0])*PW2(d[0])+2*(b[1]*PW3(c[1])-3*b[1]*c[0]*c[1]*c[2]+3*b[1]*PW2(c[0])*c[3]-3*(a[3]*b[1]+a[2]*b[2]+a[1]*b[3])*PW2(c[0]))*d[0]-(2*b[1]*c[0]*PW2(c[1])-4*b[1]*PW2(c[0])*c[2]+(2*a[1]*PW2(c[1])+a[1]*c[0]*c[2])*d[0])*d[1]+(2*b[1]*PW2(c[0])*c[1]+a[1]*c[0]*c[1]*d[0]-3*a[1]*PW2(c[0])*d[1])*d[2])*d[3];
        cf[8]=3*b[2]*PW2(b[3])*PW3(c[0])-a[2]*PW2(c[0])*PW3(d[2])-3*(a[2]*PW2(a[3])-2*a[2]*a[3]*c[3]+a[2]*PW2(c[3]))*PW3(d[0])-2*(a[2]*a[3]*c[0]-a[2]*c[0]*c[3])*PW3(d[1])+(b[2]*PW3(c[2])+3*b[2]*c[0]*PW2(c[3])+3*(a[3]*b[2]+a[2]*b[3])*c[1]*c[2]+3*(PW2(a[3])*b[2]+2*a[2]*a[3]*b[3])*c[0]-3*(b[2]*c[1]*c[2]+2*(a[3]*b[2]+a[2]*b[3])*c[0])*c[3])*PW2(d[0])+(b[2]*c[0]*PW2(c[2])-2*b[2]*c[0]*c[1]*c[3]+2*(a[3]*b[2]+a[2]*b[3])*c[0]*c[1]+2*(a[2]*a[3]*c[1]-a[2]*c[1]*c[3])*d[0])*PW2(d[1])+(b[2]*PW2(c[0])*c[2]+a[2]*c[0]*c[1]*d[1]-(a[2]*PW2(c[1])-2*a[2]*c[0]*c[2])*d[0])*PW2(d[2])+3*(b[2]*PW3(c[0])-a[2]*PW2(c[0])*d[0])*PW2(d[3])+(2*b[2]*b[3]*PW3(c[1])-6*b[2]*b[3]*c[0]*c[1]*c[2]+6*b[2]*b[3]*PW2(c[0])*c[3]-3*(2*a[3]*b[2]*b[3]+a[2]*PW2(b[3]))*PW2(c[0]))*d[0]-(2*b[2]*b[3]*c[0]*PW2(c[1])-4*b[2]*b[3]*PW2(c[0])*c[2]+2*(a[2]*a[3]*c[2]-a[2]*c[2]*c[3])*PW2(d[0])+(b[2]*c[1]*PW2(c[2])+2*(a[3]*b[2]+a[2]*b[3])*PW2(c[1])+(a[3]*b[2]+a[2]*b[3])*c[0]*c[2]-(2*b[2]*PW2(c[1])+b[2]*c[0]*c[2])*c[3])*d[0])*d[1]+(2*b[2]*b[3]*PW2(c[0])*c[1]-a[2]*c[0]*c[2]*PW2(d[1])-(4*a[2]*a[3]*c[1]+a[2]*PW2(c[2])-4*a[2]*c[1]*c[3])*PW2(d[0])+(b[2]*PW2(c[1])*c[2]-2*b[2]*c[0]*PW2(c[2])-b[2]*c[0]*c[1]*c[3]+(a[3]*b[2]+a[2]*b[3])*c[0]*c[1])*d[0]-(b[2]*c[0]*c[1]*c[2]-3*b[2]*PW2(c[0])*c[3]+3*(a[3]*b[2]+a[2]*b[3])*PW2(c[0])-(6*a[2]*a[3]*c[0]+a[2]*c[1]*c[2]-6*a[2]*c[0]*c[3])*d[0])*d[1])*d[2]-(6*b[2]*b[3]*PW3(c[0])+2*a[2]*c[0]*c[1]*PW2(d[1])+3*(2*a[2]*a[3]*c[0]+a[2]*c[1]*c[2]-2*a[2]*c[0]*c[3])*PW2(d[0])+2*(b[2]*PW3(c[1])-3*b[2]*c[0]*c[1]*c[2]+3*b[2]*PW2(c[0])*c[3]-3*(a[3]*b[2]+a[2]*b[3])*PW2(c[0]))*d[0]-(2*b[2]*c[0]*PW2(c[1])-4*b[2]*PW2(c[0])*c[2]+(2*a[2]*PW2(c[1])+a[2]*c[0]*c[2])*d[0])*d[1]+(2*b[2]*PW2(c[0])*c[1]+a[2]*c[0]*c[1]*d[0]-3*a[2]*PW2(c[0])*d[1])*d[2])*d[3];
        cf[9]=PW3(b[3])*PW3(c[0])-PW3(c[0])*PW3(d[3])-(PW3(a[3])-3*PW2(a[3])*c[3]+3*a[3]*PW2(c[3])-PW3(c[3]))*PW3(d[0])-(PW2(a[3])*c[0]-2*a[3]*c[0]*c[3]+c[0]*PW2(c[3]))*PW3(d[1])-(a[3]*PW2(c[0])-PW2(c[0])*c[3])*PW3(d[2])+(3*PW2(a[3])*b[3]*c[0]+3*a[3]*b[3]*c[1]*c[2]+b[3]*PW3(c[2])+3*b[3]*c[0]*PW2(c[3])-3*(2*a[3]*b[3]*c[0]+b[3]*c[1]*c[2])*c[3])*PW2(d[0])+(2*a[3]*b[3]*c[0]*c[1]+b[3]*c[0]*PW2(c[2])-2*b[3]*c[0]*c[1]*c[3]+(PW2(a[3])*c[1]-2*a[3]*c[1]*c[3]+c[1]*PW2(c[3]))*d[0])*PW2(d[1])+(b[3]*PW2(c[0])*c[2]-(a[3]*PW2(c[1])-2*a[3]*c[0]*c[2]-(PW2(c[1])-2*c[0]*c[2])*c[3])*d[0]+(a[3]*c[0]*c[1]-c[0]*c[1]*c[3])*d[1])*PW2(d[2])+(3*b[3]*PW3(c[0])+PW2(c[0])*c[1]*d[2]-(3*a[3]*PW2(c[0])-PW3(c[1])+3*c[0]*c[1]*c[2]-3*PW2(c[0])*c[3])*d[0]-(c[0]*PW2(c[1])-2*PW2(c[0])*c[2])*d[1])*PW2(d[3])-(3*a[3]*PW2(b[3])*PW2(c[0])-PW2(b[3])*PW3(c[1])+3*PW2(b[3])*c[0]*c[1]*c[2]-3*PW2(b[3])*PW2(c[0])*c[3])*d[0]-(PW2(b[3])*c[0]*PW2(c[1])-2*PW2(b[3])*PW2(c[0])*c[2]+(PW2(a[3])*c[2]-2*a[3]*c[2]*c[3]+c[2]*PW2(c[3]))*PW2(d[0])+(2*a[3]*b[3]*PW2(c[1])+a[3]*b[3]*c[0]*c[2]+b[3]*c[1]*PW2(c[2])-(2*b[3]*PW2(c[1])+b[3]*c[0]*c[2])*c[3])*d[0])*d[1]+(PW2(b[3])*PW2(c[0])*c[1]-(2*PW2(a[3])*c[1]+a[3]*PW2(c[2])+2*c[1]*PW2(c[3])-(4*a[3]*c[1]+PW2(c[2]))*c[3])*PW2(d[0])-(a[3]*c[0]*c[2]-c[0]*c[2]*c[3])*PW2(d[1])+(a[3]*b[3]*c[0]*c[1]+b[3]*PW2(c[1])*c[2]-2*b[3]*c[0]*PW2(c[2])-b[3]*c[0]*c[1]*c[3])*d[0]-(3*a[3]*b[3]*PW2(c[0])+b[3]*c[0]*c[1]*c[2]-3*b[3]*PW2(c[0])*c[3]-(3*PW2(a[3])*c[0]+a[3]*c[1]*c[2]+3*c[0]*PW2(c[3])-(6*a[3]*c[0]+c[1]*c[2])*c[3])*d[0])*d[1])*d[2]-(3*PW2(b[3])*PW3(c[0])+PW2(c[0])*c[2]*PW2(d[2])+(3*PW2(a[3])*c[0]+3*a[3]*c[1]*c[2]+PW3(c[2])+3*c[0]*PW2(c[3])-3*(2*a[3]*c[0]+c[1]*c[2])*c[3])*PW2(d[0])+(2*a[3]*c[0]*c[1]+c[0]*PW2(c[2])-2*c[0]*c[1]*c[3])*PW2(d[1])-2*(3*a[3]*b[3]*PW2(c[0])-b[3]*PW3(c[1])+3*b[3]*c[0]*c[1]*c[2]-3*b[3]*PW2(c[0])*c[3])*d[0]-(2*b[3]*c[0]*PW2(c[1])-4*b[3]*PW2(c[0])*c[2]+(2*a[3]*PW2(c[1])+a[3]*c[0]*c[2]+c[1]*PW2(c[2])-(2*PW2(c[1])+c[0]*c[2])*c[3])*d[0])*d[1]+(2*b[3]*PW2(c[0])*c[1]+(a[3]*c[0]*c[1]+PW2(c[1])*c[2]-2*c[0]*PW2(c[2])-c[0]*c[1]*c[3])*d[0]-(3*a[3]*PW2(c[0])+c[0]*c[1]*c[2]-3*PW2(c[0])*c[3])*d[1])*d[2])*d[3];
    }else if(abs(c[1])>=EZERO||abs(d[1])>=EZERO){
        // 曲線1は2次
        cf[0]=0;
        cf[1]=0;
        cf[2]=0;
        cf[3]=PW2(b[0])*PW2(c[1])-2*a[0]*b[0]*c[1]*d[1]+PW2(a[0])*PW2(d[1]);
        cf[4]=2*b[0]*b[1]*PW2(c[1])+2*a[0]*a[1]*PW2(d[1])-2*(a[1]*b[0]+a[0]*b[1])*c[1]*d[1];
        cf[5]=(PW2(b[1])+2*b[0]*b[2])*PW2(c[1])-2*(a[2]*b[0]+a[1]*b[1]+a[0]*b[2])*c[1]*d[1]+(PW2(a[1])+2*a[0]*a[2])*PW2(d[1]);
        cf[6]=-a[0]*c[1]*PW2(d[2])+2*(b[1]*b[2]+b[0]*b[3])*PW2(c[1])+2*(a[1]*a[2]+a[0]*a[3]-a[0]*c[3])*PW2(d[1])-(b[0]*PW2(c[2])-2*b[0]*c[1]*c[3]+2*(a[3]*b[0]+a[2]*b[1]+a[1]*b[2]+a[0]*b[3])*c[1])*d[1]+(b[0]*c[1]*c[2]+a[0]*c[2]*d[1])*d[2]-2*(b[0]*PW2(c[1])-a[0]*c[1]*d[1])*d[3];
        cf[7]=-a[1]*c[1]*PW2(d[2])+(PW2(b[2])+2*b[1]*b[3])*PW2(c[1])+(PW2(a[2])+2*a[1]*a[3]-2*a[1]*c[3])*PW2(d[1])-(b[1]*PW2(c[2])-2*b[1]*c[1]*c[3]+2*(a[3]*b[1]+a[2]*b[2]+a[1]*b[3])*c[1])*d[1]+(b[1]*c[1]*c[2]+a[1]*c[2]*d[1])*d[2]-2*(b[1]*PW2(c[1])-a[1]*c[1]*d[1])*d[3];
        cf[8]=2*b[2]*b[3]*PW2(c[1])-a[2]*c[1]*PW2(d[2])+2*(a[2]*a[3]-a[2]*c[3])*PW2(d[1])-(b[2]*PW2(c[2])-2*b[2]*c[1]*c[3]+2*(a[3]*b[2]+a[2]*b[3])*c[1])*d[1]+(b[2]*c[1]*c[2]+a[2]*c[2]*d[1])*d[2]-2*(b[2]*PW2(c[1])-a[2]*c[1]*d[1])*d[3];
        cf[9]=PW2(b[3])*PW2(c[1])+PW2(c[1])*PW2(d[3])+(PW2(a[3])-2*a[3]*c[3]+PW2(c[3]))*PW2(d[1])-(a[3]*c[1]-c[1]*c[3])*PW2(d[2])-(2*a[3]*b[3]*c[1]+b[3]*PW2(c[2])-2*b[3]*c[1]*c[3])*d[1]+(b[3]*c[1]*c[2]+(a[3]*c[2]-c[2]*c[3])*d[1])*d[2]-(2*b[3]*PW2(c[1])+c[1]*c[2]*d[2]-(2*a[3]*c[1]+PW2(c[2])-2*c[1]*c[3])*d[1])*d[3];
    }else{
        // 曲線1は1次
        cf[0]=0;
        cf[1]=0;
        cf[2]=0;
        cf[3]=0;
        cf[4]=0;
        cf[5]=0;
        cf[6]=-b[0]*c[2]+a[0]*d[2];
        cf[7]=-b[1]*c[2]+a[1]*d[2];
        cf[8]=-b[2]*c[2]+a[2]*d[2];
        cf[9]=-b[3]*c[2]+(a[3]-c[3])*d[2]+c[2]*d[3];
    }
}

// 関数値を求める
double p_(
  double *ca,	// IN  関数の係数列 ca[0]が9次項の係数
  double x)		// IN  x
{
    return ((((((((ca[0]*x+ca[1])*x+ca[2])*x+ca[3])*x+ca[4])*x+ca[5])*x+ca[6])*x+ca[7])*x+ca[8])*x+ca[9];
}

// 微分値を求める
double pd(
  double *ca,	// IN  関数の係数列 ca[0]が9次項の係数
  double x)		// IN  x
{
    return (((((((9*ca[0]*x+8*ca[1])*x+7*ca[2])*x+6*ca[3])*x+5*ca[4])*x+4*ca[5])*x+3*ca[6])*x+2*ca[7])*x+ca[8];
}

// 係数列の符号変化の回数を求める
int nSignVariation(
  double *a,	// IN  係数列
  int n)		// IN  係数列のサイズ-1
{
    int rv=0;
    int sign=0;
    for(int i=0;i<=n;i++){
        if((sign==1&&*a<0)||(sign==-1&&*a>0))
            ++rv;
        if(*a>0)
            sign=1;
        else if(*a<0)
            sign=-1;
        ++a;
    }
    return rv;
}

// Budanの定理のp(x+1)の係数を求める
void coefBudan(
  double *ca,	// IN  関数の係数列 ca[0]が9次項の係数
  double *cb)	// OUT p(x+1)の係数列  cb[0]が9次項の係数
{
    cb[0]=ca[0];
    cb[1]=9*ca[0]+ca[1];
    cb[2]=36*ca[0]+8*ca[1]+ca[2];
    cb[3]=84*ca[0]+28*ca[1]+7*ca[2]+ca[3];
    cb[4]=126*ca[0]+56*ca[1]+21*ca[2]+6*ca[3]+ca[4];
    cb[5]=126*ca[0]+70*ca[1]+35*ca[2]+15*ca[3]+5*ca[4]+ca[5];
    cb[6]=84*ca[0]+56*ca[1]+35*ca[2]+20*ca[3]+10*ca[4]+4*ca[5]+ca[6];
    cb[7]=36*ca[0]+28*ca[1]+21*ca[2]+15*ca[3]+10*ca[4]+6*ca[5]+3*ca[6]+ca[7];
    cb[8]=9*ca[0]+8*ca[1]+7*ca[2]+6*ca[3]+5*ca[4]+4*ca[5]+3*ca[6]+2*ca[7]+ca[8];
    cb[9]=ca[0]+ca[1]+ca[2]+ca[3]+ca[4]+ca[5]+ca[6]+ca[7]+ca[8]+ca[9];
}

// VAG法の係数を求める
void coefVAG(
  double *ca,	// IN  関数の係数列 ca[0]が9次項の係数
  double a,		// IN  区間の端点
  double b,
  double *cv)	// OUT VAG係数列  cv[0]が9次項の係数
{
    cv[0]=((((((((ca[0]*b+ca[1])*b+ca[2])*b+ca[3])*b+ca[4])*b+ca[5])*b+ca[6])*b+ca[7])*b+ca[8])*b+ca[9];
    cv[1]=((((((((9*a*ca[0]+ca[1])*b+8*a*ca[1]+2*ca[2])*b+7*a*ca[2]+3*ca[3])*b+6*a*ca[3]+4*ca[4])*b+5*a*ca[4]+5*ca[5])*b+4*a*ca[5]+6*ca[6])*b+3*a*ca[6]+7*ca[7])*b+2*a*ca[7]+8*ca[8])*b+a*ca[8]+9*ca[9];
    cv[2]=((((((((36*a*ca[0]+8*ca[1])*a+ca[2])*b+(28*a*ca[1]+14*ca[2])*a+3*ca[3])*b+(21*a*ca[2]+18*ca[3])*a+6*ca[4])*b+(15*a*ca[3]+20*ca[4])*a+10*ca[5])*b+(10*a*ca[4]+20*ca[5])*a+15*ca[6])*b+(6*a*ca[5]+18*ca[6])*a+21*ca[7])*b+(3*a*ca[6]+14*ca[7])*a+28*ca[8])*b+(a*ca[7]+8*ca[8])*a+36*ca[9];
    cv[3]=((((((((84*a*ca[0]+28*ca[1])*a+7*ca[2])*a+ca[3])*b+((56*a*ca[1]+42*ca[2])*a+18*ca[3])*a+4*ca[4])*b+((35*a*ca[2]+45*ca[3])*a+30*ca[4])*a+10*ca[5])*b+((20*a*ca[3]+40*ca[4])*a+40*ca[5])*a+20*ca[6])*b+((10*a*ca[4]+30*ca[5])*a+45*ca[6])*a+35*ca[7])*b+((4*a*ca[5]+18*ca[6])*a+42*ca[7])*a+56*ca[8])*b+((a*ca[6]+7*ca[7])*a+28*ca[8])*a+84*ca[9];
    cv[4]=((((((((126*a*ca[0]+56*ca[1])*a+21*ca[2])*a+6*ca[3])*a+ca[4])*b+(((70*a*ca[1]+70*ca[2])*a+45*ca[3])*a+20*ca[4])*a+5*ca[5])*b+(((35*a*ca[2]+60*ca[3])*a+60*ca[4])*a+40*ca[5])*a+15*ca[6])*b+(((15*a*ca[3]+40*ca[4])*a+60*ca[5])*a+60*ca[6])*a+35*ca[7])*b+(((5*a*ca[4]+20*ca[5])*a+45*ca[6])*a+70*ca[7])*a+70*ca[8])*b+(((a*ca[5]+6*ca[6])*a+21*ca[7])*a+56*ca[8])*a+126*ca[9];
    cv[5]=((((((((126*ca[0]*b+56*ca[1])*b+21*ca[2])*b+6*ca[3])*b+ca[4])*a+(((70*ca[1]*b+70*ca[2])*b+45*ca[3])*b+20*ca[4])*b+5*ca[5])*a+(((35*ca[2]*b+60*ca[3])*b+60*ca[4])*b+40*ca[5])*b+15*ca[6])*a+(((15*ca[3]*b+40*ca[4])*b+60*ca[5])*b+60*ca[6])*b+35*ca[7])*a+(((5*ca[4]*b+20*ca[5])*b+45*ca[6])*b+70*ca[7])*b+70*ca[8])*a+(((ca[5]*b+6*ca[6])*b+21*ca[7])*b+56*ca[8])*b+126*ca[9];
    cv[6]=((((((((84*ca[0]*b+28*ca[1])*b+7*ca[2])*b+ca[3])*a+((56*ca[1]*b+42*ca[2])*b+18*ca[3])*b+4*ca[4])*a+((35*ca[2]*b+45*ca[3])*b+30*ca[4])*b+10*ca[5])*a+((20*ca[3]*b+40*ca[4])*b+40*ca[5])*b+20*ca[6])*a+((10*ca[4]*b+30*ca[5])*b+45*ca[6])*b+35*ca[7])*a+((4*ca[5]*b+18*ca[6])*b+42*ca[7])*b+56*ca[8])*a+((ca[6]*b+7*ca[7])*b+28*ca[8])*b+84*ca[9];
    cv[7]=((((((((36*ca[0]*b+8*ca[1])*b+ca[2])*a+(28*ca[1]*b+14*ca[2])*b+3*ca[3])*a+(21*ca[2]*b+18*ca[3])*b+6*ca[4])*a+(15*ca[3]*b+20*ca[4])*b+10*ca[5])*a+(10*ca[4]*b+20*ca[5])*b+15*ca[6])*a+(6*ca[5]*b+18*ca[6])*b+21*ca[7])*a+(3*ca[6]*b+14*ca[7])*b+28*ca[8])*a+(ca[7]*b+8*ca[8])*b+36*ca[9];
    cv[8]=((((((((9*ca[0]*b+ca[1])*a+8*ca[1]*b+2*ca[2])*a+7*ca[2]*b+3*ca[3])*a+6*ca[3]*b+4*ca[4])*a+5*ca[4]*b+5*ca[5])*a+4*ca[5]*b+6*ca[6])*a+3*ca[6]*b+7*ca[7])*a+2*ca[7]*b+8*ca[8])*a+ca[8]*b+9*ca[9];
    cv[9]=((((((((ca[0]*a+ca[1])*a+ca[2])*a+ca[3])*a+ca[4])*a+ca[5])*a+ca[6])*a+ca[7])*a+ca[8])*a+ca[9];
}

// Budanの定理による実根の有無検定
int vBudan(     // RETURN p(x)とp(x+1)の係数の符号変化の数の差
  double *ca)	// IN  関数の係数列 ca[0]が9次項の係数
{
    double cb[10];
    coefBudan(ca,cb);
    return nSignVariation(ca,9)-nSignVariation(cb,9);
}

// VAG法により実根が1個含まれる区間を求める
// ※注意事項：
//   この関数は再起呼び出しされます。
//   ABMIN（区間長の最小値）の設定や処理系によってはスタックがオーバーフローする可能性があります。
//   必要に応じスタックサイズを調節するなどしてください。
void vVAG(
  double *ca,			// IN  関数の係数列 ca[0]が9次項の係数
  Interval &iv,			// IN  対象の区間
  vector<Interval> &ov)	// OUT 実根が1個含まれる区間
{
    Interval nv;
    double cv[10];
    coefVAG(ca,iv.a,iv.b,cv);
    int var=nSignVariation(cv,9);
    double m=(iv.a+iv.b)/2;
    if(var==0){
        return;
    }else if(var==1){
        nv.a=iv.a;
        nv.b=iv.b;
        nv.pa=p_(ca,iv.a);
        nv.pb=p_(ca,iv.b);
        nv.nr=var;
        ov.push_back(nv);
        return;
    }else if(iv.b-iv.a<ABMIN){
        // 重根
        // 区間をABADJ*2に拡張し、重複数を求める
        double a=m-ABADJ;
        double b=m+ABADJ;
        double pa=p_(ca,a);
        double pb=p_(ca,b);
        coefVAG(ca,a,b,cv);
        var=nSignVariation(cv,9);
        nv.a=a;
        nv.b=b;
        nv.pa=pa;
        nv.pb=pb;
        nv.nr=var;
        ov.push_back(nv);
        return;
    }
    double pm=p_(ca,m);
    if(abs(pm)<EZERO){
        double a=m-ABADJ;
        double b=m+ABADJ;
        double pa=p_(ca,a);
        double pb=p_(ca,b);
        nv.a=iv.a;
        nv.b=a;
        nv.pa=iv.pa;
        nv.pb=pa;
        vVAG(ca,nv,ov);

        coefVAG(ca,a,b,cv);
        var=nSignVariation(cv,9);
        nv.nr=var;
        if(var==1){
            nv.a=nv.b=m;
            nv.pa=nv.pb=pm;
        }else{
            // 重根
            nv.a=a;
            nv.b=b;
            nv.pa=pa;
            nv.pb=pb;
        }
        ov.push_back(nv);

        nv.a=b;
        nv.b=iv.b;
        nv.pa=pb;
        nv.pb=iv.pb;
        vVAG(ca,nv,ov);
    }else{
        nv.a=iv.a;
        nv.b=m;
        nv.pa=iv.pa;
        nv.pb=pm;
        vVAG(ca,nv,ov);

        nv.a=m;
        nv.b=iv.b;
        nv.pa=pm;
        nv.pb=iv.pb;
        vVAG(ca,nv,ov);
    }
}

// 二分法により実根が1個含まれる区間を1/2に縮小する
void vBisection(
  double *ca,		// IN  関数の係数列 ca[0]が9次項の係数
  Interval &v)		// IN/OUT 対象の区間
{
    double m=(v.a+v.b)/2;
    double pm=p_(ca,m);
    if(abs(pm)<EZERO){
        v.a=v.b=m;
        v.pa=v.pb=0.0;
        return;
    }else if(((v.pa>=EZERO||(abs(v.pa)<EZERO&&pd(ca,v.a)>0.0))&&pm<=-EZERO)||
             ((v.pa<=-EZERO||(abs(v.pa)<EZERO&&pd(ca,v.a)<0.0))&&pm>=EZERO)){
        v.b=m;
        v.pb=pm;
    }else{
        v.a=m;
        v.pa=pm;
    }
    return;
}

// Ridder's法により9次関数の実根を求める
// 出典：Kiusalaas, Jaan (2010). Numerical Methods in Engineering with Python (2nd ed.). Cambridge University Press. pp. 147–148. ISBN 978-0-521-19132-6
double ridder(      // RETURN 根の近似値
  double *ca,		// IN  関数の係数列 ca[0]が9次項の係数
  Interval &v,		// IN  対象の区間
  double tol)		// IN  収束判定しきい値
{
    double a=v.a;
    double b=v.b;
    double fa=v.pa;
    double fb=v.pb;
    if(fa==0.0)
        return a;
    if(fb==0.0)
        return b;
    if(fa*fb>0.0){
        return nan("");	// root is not bracketed
    }
    double xOld;
    for(int i=0;i<30;i++){
        // Compute the improved root x from Ridder's formula
        double c=0.5*(a+b);
        double fc=p_(ca,c);
        double s=sqrt(fc*fc-fa*fb);
        if(s==0)
            return nan("");
        double dx=(c-a)*fc/s;
        if(fa-fb<0.0)
            dx=-dx;
        double x=c+dx;
        double fx=p_(ca,x);
        // Test for convergence
        if(i>0&&abs(x-xOld)<tol)
            return x;
        xOld=x;
        // Re-bracket the root as tightly as possible
        if(fc*fx>0.0){
            if(fa*fx<0.0){
                b=x;
                fb=fx;
            }else{
                a=x;
                fa=fx;
            }
        }else{
            a=c;
            b=x;
            fa=fc;
            fb=fx;
        }
    }
    return nan("");	// Too many iterations
}

// 9次関数の区間[0,1}の実根を求める
void vRoots(
  double *ca,			// IN  関数の係数列 ca[0]が9次項の係数
  vector<double> &rt)	// OUT 実根
{
    rt.clear();

    Interval iv;			// 根を求める区間
    iv.a=0.0;
    iv.b=1.0;
    iv.pa=p_(ca,0.0);
    iv.pb=p_(ca,1.0);

    vector<Interval> rv;	// 1つの根が含まれる区間
    rv.reserve(9);

    if(abs(iv.pa)<EZERO){
        // 始点の関数値が0なので、0を根とする
        Interval nv;
        nv.a=nv.b=0.0;
        nv.pa=nv.pb=iv.pa;
        rv.push_back(nv);
    }

    // {0,1}の実根を求める
    int nb=vBudan(ca);	// Budanの定理による実根の有無検定
    if(nb==1){
        // {0,1}に1個の実根がある
        rv.push_back(iv);
    }else if(nb>=2){
        // {0,1}に1個以上の実根がある
        vVAG(ca,iv,rv);
    }
    if(rv.size()==0){
        // 交点なし、終了
        return;
    }
    // 二分法により実根が1個含まれる区間を1/8に縮小する
    // 重根の場合はスキップ
    for(Interval& v:rv){
        if(v.nr==1&&v.a<v.b){
            vBisection(ca,v);
            vBisection(ca,v);
            vBisection(ca,v);
        }
    }
    for(Interval& v:rv){
        if(v.a==v.b){
            rt.push_back(v.a);
        }else if(v.nr==1){
            // 重複なし  Ridder's法により実根を求める
            double r=ridder(ca,v,ABMIN);
            if(!isnan(r)){
                rt.push_back(r);
            }
        }else{
            // 重複あり  区間の中央を根とする
            rt.push_back((v.a+v.b)/2);
        }
    }
}

// 3次ベジェ曲線の交点座標を求める
bool vBezierIntersection(   // RETURN true:正常
  QList<QPointF> &cp0,	// IN  3次ベジェ曲線#0 制御点
  QList<QPointF> &cp1,	// IN  3次ベジェ曲線#1 制御点
  QList<QPointF> &op)	// OUT 交点座標
{
    op.clear();
    if(cp0.size()!=4||cp1.size()!=4)
        return false;

    // 制御点の正規化座標を求める
    QList<QPointF> ncp0;
    QList<QPointF> ncp1;
    normCP(cp0, cp1, ncp0, ncp1);

    // 正規化したベジェ曲線の係数を求める
    double na[4],nb[4],nc[4],nd[4];
    coefBezier(ncp0,na,nb);
    coefBezier(ncp1,nc,nd);

    // 交点の陰伏方程式F0(u),F1(t)の係数を求める
    double cf0[10],cf1[10];
    coefImp0(na,nb,nc,nd,cf0);
    coefImp1(na,nb,nc,nd,cf1);

    vector<double> rt0,rt1;		// 根
    rt0.reserve(9);
    rt1.reserve(9);

    // F0(u),F1(t)の根を求める
    vRoots(cf0,rt0);
    vRoots(cf1,rt1);

    QList<QPointF> nip0,nip1;		// 交点候補（正規化座標）
    for(double u:rt0)
        nip0<<pointBezier(nc,nd,u);
    for(double t:rt1)
        nip1<<pointBezier(na,nb,t);

    // 交点候補の座標を突き合わせる
    double c[4],d[4];
    coefBezier(cp1,c,d);
    int i0=0;
    foreach(QPointF p0,nip0){
        foreach(QPointF p1,nip1){
            double dx=p1.x()-p0.x();
            double dy=p1.y()-p0.y();
            if(sqrt(PW2(dx)+PW2(dy))<=DSTMIN){
                // 交点座標（正規化なし）を出力
                op<<pointBezier(c,d,rt0[i0]);
            }
        }
        ++i0;
    }
    return true;
}

コメントにも記載していますが、vVAGは再起呼び出しを用いるので反復回数（ネストの深さ）に応じてスタックフレームを確保する必要があります。この実装例でvVAGが使用するスタックを実測したところ4624バイトでした。MSVCではあまり気にしなくても大丈夫だと思いますが、gccではスタックサイズがきっちり計算・確保される傾向があり対処必要かもしれません。

関連記事

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?