(((silence)))

// 内部使用，计算阶次
static inline double f ( double x, int n )
{
    double db=1;
    if (n==0) {
        db = 1;
    } else if (n>0) {
        for(;n--;db*=x);
    } else {
        for(;n++;db/=x);
    }
    return db;
}

// 最小二乘法拟合一元多项式曲线
// 返回值：
//     成功返回对应的多项式系数
//     失败返回空向量
NA::VecDP leastSquare ( NA::VecDP const& vx, // x 坐标
                        NA::VecDP const& vy, // 对应的 y 坐标
                        int power )          // 阶次
{
    int dim = power+1;
    int len = vx.size();
    NA_ASSERT( len == vy.size() );
    NA::MatDP mx(0.0, dim, dim), my(0.0, dim, 1);
    int i, j, k;
    double db;
    // init {{{
    /*
       mx =
                                          2
                  1           X          X
            -----------------------------------
           |                                  2
        1  | sum(1 x 1)  sum(1 x X)  sum(1 x X )
           |
           |                                  2
        X  | sum(X x 1)  sum(X x X)  sum(X x X )
           |
         2 |      2           2           2   2
        X  | sum(X x 1)  sum(X x X)  sum(X x X )
           |
     */
    for (i=0; i<dim; ++i) {
        for (j=i; j<dim; ++j) {
            db=0;
            for (k=0; k<len; ++k) {
                db += f(vx[k], i) * f(vx[k], j);
            }
            mx[i][j]=mx[j][i]=db;
        }
    }

    /*
       my =

                   Y
           ----------------
        1  |  sum(X x Y)
           |
           |
        X  |  sum(X x Y)
           |
         2 |       2
        X  |  sum(X x Y)
           |
     */
    for (j=0; j<dim; ++j) {
        db = 0.0;
        for (k=0; k<len; ++k) {
            db += vy[k] * f(vx[k], j);
        }
        my[j][0]=db;
    }
    // }}}

    // have a look {{{
#if DEBUG
    for (i=0; i<dim; ++i) {
        cout<<"| ";
        for (j=0; j<dim; ++j)
            cout<<mx[i][j]<<"\t";
        cout<<"|\t"<<my[i][0]<<endl;
    }
#endif
    // }}}

    if (!gaussianBackSubstitution (mx, my)) { // 直接用高斯消元法求解
        return NA::VecDP(0);
    } else {
        return NA::VecDP(my.bits(), dim);
    }
}