(((silence)))

// 简单的求导
static inline double df ( NA::functionDP_DP f, double x )
{
    const double eps=1.0e-4;
    return (f(x+eps)-f(x)) / eps;
}

// 牛顿迭代法求解非线性方程
// f:   函数
// y:   对应的y值
// x0:  x 的初始值
// eps: x 的精确度
// eta: y 的精确度
// N:   允许的最大迭代次数
// 返回值: 在 [a,b] 区间内使 |f(x)-y|<eps 成立的 x 值。
double newtonSolve(
                   NA::functionDP_DP f,
                   double y,
                   double x0,
                   double eps,
                   double eta,
                   int N
) {
    double x=x0;
    for( ; fabs(f(x)-y)>=eta; ) {
        x-=f(x)/df(f,x);
        if ( fabs(x-x0)<eps ) {
            return x; // throw ( NA::Error("newtonSolve: No solution") ); ??
        }
        x0=x;
    }
    return x;
}