Horner scheme问题

本文介绍了数值分析中用于高效评估多项式的Horner算法。该算法不仅可用于手动近似多项式方程的根,还可用于不同进制数系统之间的转换,并且在矩阵计算中也有广泛应用。文中提供了一个简单的C++程序示例,演示了如何使用Horner算法来计算多项式的值。

Description

In numerical analysis, the Horner scheme or Horner algorithm, named after William George Horner, is an algorithm for the efficient evaluation of polynomials in monomial form. Horner’s method describes a manual process by which one may approximate the roots of a polynomial equation. The Horner scheme can also be viewed as a fast algorithm for dividing a polynomial by a linear polynomial with Ruffini’s rule.

Application

The Horner scheme is often used to convert between different positional numeral systems — in which case x is the base of the number system, and the ai coefficients are the digits of the base-x representation of a given number — and can also be used if x is a matrix, in which case the gain in computational efficiency is even greater.

History

Even though the algorithm is named after William George Horner, who described it in 1819, the method was already known to Isaac Newton in 1669, and even earlier to the Chinese mathematician Ch’in Chiu-Shao in the 13th century. TASK: write a program to calculate sum of Polynomial by Horner scheme.

Input

tow lines. The first line have tow numbers,n and x, n<=20, x<=10 The second line have n+1 numbers, a0,a1…an.

Output

The sum of Polynomial

Sample Input

5 2

0 1 2 3 4 5

Sample Output

258

#include<iostream>
using namespace std;
int qinjiushao(int a[], int n, int x)
{
    int sum = a[n];
    for (int i = n - 1; i >= 0; i--)
    {
        sum = sum*x + a[i];
    }
    return sum;
}
int main()
{
    int arr[1005], n, x;
    cin >> n >> x;
    for (int i = 0; i < n+1;i++)
        cin >> arr[i];
    cout << qinjiushao(arr, n, x) << endl;
    return 0;
}

#ifndef ROOT_FINDER_HPP #define ROOT_FINDER_HPP #define _USE_MATH_DEFINES #include <cfloat> #include <cmath> #include <set> #include <Eigen/Eigen> namespace RootFinderParam { constexpr size_t highestOrder = 64; } namespace RootFinderPriv { inline int polyMod(double *u, double *v, double *r, int lu, int lv) // Modulus of u(x)/v(x) // The leading coefficient of v, i.e., v[0], must be 1.0 or -1.0 // The length of u, v, and r are lu, lv, and lu, respectively { int orderu = lu - 1; int orderv = lv - 1; memcpy(r, u, lu * sizeof(double)); if (v[0] < 0.0) { for (int i = orderv + 1; i <= orderu; i += 2) { r[i] = -r[i]; } for (int i = 0; i <= orderu - orderv; i++) { for (int j = i + 1; j <= orderv + i; j++) { r[j] = -r[j] - r[i] * v[j - i]; } } } else { for (int i = 0; i <= orderu - orderv; i++) { for (int j = i + 1; j <= orderv + i; j++) { r[j] = r[j] - r[i] * v[j - i]; } } } int k = orderv - 1; while (k >= 0 && fabs(r[orderu - k]) < DBL_EPSILON) { r[orderu - k] = 0.0; k--; } return (k <= 0) ? 1 : (k + 1); } inline double polyEval(double *p, int len, double x) // Evaluate the polynomial p(x), which has len coefficients // Note: Horner scheme should not be employed here !!! // Horner scheme has bad numerical stability despite of its efficiency. // These errors are particularly troublesome for root-finding algorithms. // When the polynomial is evaluated near a zero, catastrophic // cancellation (subtracting two nearby numbers) is guaranteed to occur. // Therefore, Horner scheme may slow down some root-finding algorithms. { double retVal = 0.0; if (len > 0) { if (fabs(x) < DBL_EPSILON) { retVal = p[len - 1]; } else if (x == 1.0) { for (int i = len - 1; i >= 0; i--) { retVal += p[i]; } } else { double xn = 1.0; for (int i = len - 1; i >= 0; i--) { retVal += p[i] * xn; xn *= x; } } } return retVal; } inline std::set<double> solveCub(double a, double b, double c, double d) // Calculate all roots of a*x^3 + b*x^2 + c*x + d = 0 { std::set<double> roots; constexpr double cos120 = -0.50; constexpr double sin120 = 0.866025403784438646764; if (fabs(d) < DBL_EPSILON) { // First solution is x = 0 roots.insert(0.0); // Converting to a quadratic equation d = c; c = b; b = a; a = 0.0; } if (fabs(a) < DBL_EPSILON) { if (fabs(b) < DBL_EPSILON) { // Linear equation if (fabs(c) > DBL_EPSILON) roots.insert(-d / c); } else { // Quadratic equation double discriminant = c * c - 4.0 * b * d; if (discriminant >= 0) { double inv2b = 1.0 / (2.0 * b); double y = sqrt(discriminant); roots.insert((-c + y) * inv2b); roots.insert((-c - y) * inv2b); } } } else { // Cubic equation double inva = 1.0 / a; double invaa = inva * inva; double bb = b * b; double bover3a = b * (1.0 / 3.0) * inva; double p = (3.0 * a * c - bb) * (1.0 / 3.0) * invaa; double halfq = (2.0 * bb * b - 9.0 * a * b * c + 27.0 * a * a * d) * (0.5 / 27.0) * invaa * inva; double yy = p * p * p / 27.0 + halfq * halfq; if (yy > DBL_EPSILON) { // Sqrt is positive: one real solution double y = sqrt(yy); double uuu = -halfq + y; double vvv = -halfq - y; double www = fabs(uuu) > fabs(vvv) ? uuu : vvv; double w = (www < 0) ? -pow(fabs(www), 1.0 / 3.0) : pow(www, 1.0 / 3.0); roots.insert(w - p / (3.0 * w) - bover3a); } else if (yy < -DBL_EPSILON) { // Sqrt is negative: three real solutions double x = -halfq; double y = sqrt(-yy); double theta; double r; double ux; double uyi; // Convert to polar form if (fabs(x) > DBL_EPSILON) { theta = (x > 0.0) ? atan(y / x) : (atan(y / x) + M_PI); r = sqrt(x * x - yy); } else { // Vertical line theta = M_PI / 2.0; r = y; } // Calculate cube root theta /= 3.0; r = pow(r, 1.0 / 3.0); // Convert to complex coordinate ux = cos(theta) * r; uyi = sin(theta) * r; // First solution roots.insert(ux + ux - bover3a); // Second solution, rotate +120 degrees roots.insert(2.0 * (ux * cos120 - uyi * sin120) - bover3a); // Third solution, rotate -120 degrees roots.insert(2.0 * (ux * cos120 + uyi * sin120) - bover3a); } else { // Sqrt is zero: two real solutions double www = -halfq; double w = (www < 0.0) ? -pow(fabs(www), 1.0 / 3.0) : pow(www, 1.0 / 3.0); // First solution roots.insert(w + w - bover3a); // Second solution, rotate +120 degrees roots.insert(2.0 * w * cos120 - bover3a); } } return roots; } 注释上述代码
09-19
inline double polyVal(const Eigen::VectorXd &coeffs, double x, bool numericalStability = true) // Evaluate the polynomial at x, i.e., coeffs(x) // Horner scheme is faster yet less stable // Stable one should be used when coeffs(x) is close to 0.0 { double retVal = 0.0; int order = (int)coeffs.size() - 1; if (order >= 0) { if (fabs(x) < DBL_EPSILON) { retVal = coeffs(order); } else if (x == 1.0) { retVal = coeffs.sum(); } else { if (numericalStability) { double xn = 1.0; for (int i = order; i >= 0; i--) { retVal += coeffs(i) * xn; xn *= x; } } else { int len = coeffs.size(); for (int i = 0; i < len; i++) { retVal = retVal * x + coeffs(i); } } } } return retVal; } inline int countRoots(const Eigen::VectorXd &coeffs, double l, double r) // Count the number of distinct roots of coeffs(x) inside (l, r), leveraging Sturm theory // Boundary values, i.e., coeffs(l) and coeffs(r), must be nonzero { int nRoots = 0; int originalSize = coeffs.size(); int valid = originalSize; for (int i = 0; i < originalSize; i++) { if (fabs(coeffs(i)) < DBL_EPSILON) { valid--; } else { break; } } if (valid > 0 && fabs(coeffs(originalSize - 1)) > DBL_EPSILON) { Eigen::VectorXd monicCoeffs(valid); monicCoeffs << 1.0, coeffs.segment(originalSize - valid + 1, valid - 1) / coeffs(originalSize - valid); // Build the Sturm sequence int len = monicCoeffs.size(); int order = len - 1; double sturmSeqs[(RootFinderParam::highestOrder + 1) * (RootFinderParam::highestOrder + 1)]; int szSeq[RootFinderParam::highestOrder + 1] = {0}; // Explicit ini as zero (gcc may neglect this in -O3) int num = 0; for (int i = 0; i < len; i++) { sturmSeqs[i] = monicCoeffs(i); sturmSeqs[i + 1 + len] = (order - i) * sturmSeqs[i] / order; } szSeq[0] = len; szSeq[1] = len - 1; num += 2; bool remainderConstant = false; int idx = 0; while (!remainderConstant) { szSeq[idx + 2] = RootFinderPriv::polyMod(&(sturmSeqs[(idx + 1) * len - szSeq[idx]]), &(sturmSeqs[(idx + 2) * len - szSeq[idx + 1]]), &(sturmSeqs[(idx + 3) * len - szSeq[idx]]), szSeq[idx], szSeq[idx + 1]); remainderConstant = szSeq[idx + 2] == 1; for (int i = 1; i < szSeq[idx + 2]; i++) { sturmSeqs[(idx + 3) * len - szSeq[idx + 2] + i] /= -fabs(sturmSeqs[(idx + 3) * len - szSeq[idx + 2]]); } sturmSeqs[(idx + 3) * len - szSeq[idx + 2]] /= -fabs(sturmSeqs[(idx + 3) * len - szSeq[idx + 2]]); num++; idx++; } // Count numbers of sign variations at two boundaries double yl, lastyl, yr, lastyr; lastyl = RootFinderPriv::polyEval(&(sturmSeqs[len - szSeq[0]]), szSeq[0], l); lastyr = RootFinderPriv::polyEval(&(sturmSeqs[len - szSeq[0]]), szSeq[0], r); for (int i = 1; i < num; i++) { yl = RootFinderPriv::polyEval(&(sturmSeqs[(i + 1) * len - szSeq[i]]), szSeq[i], l); yr = RootFinderPriv::polyEval(&(sturmSeqs[(i + 1) * len - szSeq[i]]), szSeq[i], r); if (lastyl == 0.0 || lastyl * yl < 0.0) { ++nRoots; } if (lastyr == 0.0 || lastyr * yr < 0.0) { --nRoots; } lastyl = yl; lastyr = yr; } } return nRoots; } inline std::set<double> solvePolynomial(const Eigen::VectorXd &coeffs, double lbound, double ubound, double tol, bool isolation = true) // Calculate roots of coeffs(x) inside (lbound, rbound) // // Closed-form solutions are employed for reduced_order < 5 // isolation = true: // Sturm' theory and some geometrical property are employed to bracket each root // Safe-Newton is employed to shrink the interval efficiently // isolation = false: // Eigen values of polynomial companion matrix are calculated // // Requirement: leading coefficient must be nonzero // coeffs(lbound) != 0, coeffs(rbound) != 0, lbound < rbound { std::set<double> rts; int valid = coeffs.size(); for (int i = 0; i < coeffs.size(); i++) { if (fabs(coeffs(i)) < DBL_EPSILON) { valid--; } else { break; } } int offset = 0; int nonzeros = valid; if (valid > 0) { for (int i = 0; i < valid; i++) { if (fabs(coeffs(coeffs.size() - i - 1)) < DBL_EPSILON) { nonzeros--; offset++; } else { break; } } } if (nonzeros == 0) { rts.insert(INFINITY); rts.insert(-INFINITY); } else if (nonzeros == 1 && offset == 0) { rts.clear(); } else { Eigen::VectorXd ncoeffs(std::max(5, nonzeros)); ncoeffs.setZero(); ncoeffs.tail(nonzeros) << coeffs.segment(coeffs.size() - valid, nonzeros); if (nonzeros <= 5) { rts = RootFinderPriv::solveQuart(ncoeffs(0), ncoeffs(1), ncoeffs(2), ncoeffs(3), ncoeffs(4)); } else { if (isolation) { rts = RootFinderPriv::isolateRealRoots(ncoeffs, lbound, ubound, tol); } else { rts = RootFinderPriv::eigenSolveRealRoots(ncoeffs, lbound, ubound, tol); } } if (offset > 0) { rts.insert(0.0); } } for (auto it = rts.begin(); it != rts.end();) { if (*it > lbound && *it < ubound) { it++; } else { it = rts.erase(it); } } return rts; } } // namespace RootFinder #endif 注释上述代码并保留原有注释
最新发布
09-19
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值