Program Listing for File PolynomialRoots-Utils.cc¶
↰ Return to documentation for file (PolynomialRoots-Utils.cc)
/*--------------------------------------------------------------------------*\
| |
| Copyright (C) 2017 |
| |
| , __ , __ |
| /|/ \ /|/ \ |
| | __/ _ ,_ | __/ _ ,_ |
| | \|/ / | | | | \|/ / | | | |
| |(__/|__/ |_/ \_/|/|(__/|__/ |_/ \_/|/ |
| /| /| |
| \| \| |
| |
| Enrico Bertolazzi |
| Dipartimento di Ingegneria Industriale |
| Universita` degli Studi di Trento |
| email: enrico.bertolazzi@unitn.it |
| |
\*--------------------------------------------------------------------------*/
#ifdef __GNUC__
#pragma GCC diagnostic ignored "-Wunused-function"
#endif
#ifdef __clang__
#pragma clang diagnostic ignored "-Wglobal-constructors"
#pragma clang diagnostic ignored "-Wvla-extension"
#pragma clang diagnostic ignored "-Wvla"
#pragma clang diagnostic ignored "-Wunused-function"
#pragma clang diagnostic ignored "-Wdocumentation-unknown-command"
#endif
#include "PolynomialRoots.hh"
#include "PolynomialRoots-Utils.hh"
#ifndef DOXYGEN_SHOULD_SKIP_THIS
namespace PolynomialRoots {
// static valueType const machepsi = std::numeric_limits<valueType>::epsilon();
using std::pair;
using std::abs;
using std::pow;
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// stable computation of polinomial
//
// op[0] * x^n + .... + op[n-1]*x + op[n]
//
// (op[0] + op[1]/x.... + op[n]/x^n)*x^n
//
valueType
evalPoly(
valueType const op[],
indexType Degree,
valueType x
) {
bool reverse = std::abs(x) > 1;
if ( reverse ) {
valueType res(op[Degree]);
valueType xn(1);
for ( indexType i = 1; i <= Degree; ++i ) {
res = res/x + op[Degree-i];
xn *= x;
}
res *= xn;
return res;
} else {
valueType res(op[0]);
for ( indexType i = 1; i <= Degree; ++i ) res = res*x + op[i];
return res;
}
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// stable computation of Newton step
//
// op[0] * x^n + .... + op[n-1]*x + op[n]
//
// (op[0] + op[1]/x.... + op[n]/x^n)*x^n
//
bool
NewtonStep(
valueType const op[],
indexType Degree,
valueType & x
) {
valueType p, dp;
bool reverse = std::abs(x) > 1;
if ( reverse ) {
p = op[Degree];
valueType xn(1);
for ( indexType i = 1; i <= Degree; ++i ) {
p = p/x + op[Degree-i];
xn *= x;
}
p *= xn;
dp = op[Degree];
xn = 1;
for ( indexType i = 2; i <= Degree; ++i ) {
dp = dp/x + (Degree-i)*op[Degree-i];
xn *= x;
}
dp *= xn;
} else {
dp = Degree*op[0];
p = op[0];
for ( indexType i = 1; i < Degree; ++i ) {
p = p*x + op[i];
dp = dp*x + (Degree-i)*op[i];
}
p = p*x + op[Degree];
}
x -= p/dp;
return true;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// stable computation of polinomial and its derivative
//
// op[0] * x^n + .... + op[n-1]*x + op[n]
//
// (op[0] + op[1]/x.... + op[n]/x^n)*x^n
//
void
evalPolyDPoly(
valueType const op[],
indexType Degree,
valueType x,
valueType & p,
valueType & dp
) {
bool reverse = std::abs(x) > 1;
if ( reverse ) {
valueType const * pc = op+Degree;
valueType rD = Degree;
p = *pc--;
dp = rD*p;
valueType xn(1);
for ( indexType i = 1; i < Degree; ++i ) {
p = p/x + *pc;
dp = dp/x + rD*(*pc--);
xn *= x;
}
dp *= xn;
p = p/x + *pc;
xn *= x;
p *= xn;
} else {
p = op[0]*x+op[1];
for ( indexType i = 2; i <= Degree; ++i ) {
p = p*x + op[i];
dp = dp*x + i*op[i];
}
}
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
//
// op[0] * x^n + .... + op[n-1]*x + op[n]
std::complex<valueType>
evalPolyC(
valueType const op[],
indexType Degree,
std::complex<valueType> const & x
) {
bool reverse = std::abs(x) > 1;
if ( reverse ) {
std::complex<valueType> res(op[Degree]);
std::complex<valueType> xn(1,0);
for ( indexType i = 1; i <= Degree; ++i ) {
res = res/x + op[Degree-i];
xn *= x;
}
res *= xn;
return res;
} else {
std::complex<valueType> res(op[0]);
for ( indexType i = 1; i <= Degree; ++i ) res = res*x + op[i];
return res;
}
}
//============================================================================
/*
.. scale roots by pair.second, after scaling the polinomial has the form
..
.. p[0] + p[1]*x + ... + p[n]*x^n
..
.. pair.first is the index such that p[pair.first] = +-1
..
*/
static
pair<indexType,valueType>
scalePolynomial(
indexType n, // degree
valueType const p[],
valueType ps[]
) {
indexType i_max = n;
valueType an = p[n];
valueType scale = -1;
indexType i = n;
while ( --i >= 0 ) {
ps[i] = p[i]/an;
valueType scale1 = std::pow( std::abs(ps[i]), 1.0/(n-i) );
if ( scale1 > scale ) { scale = scale1; i_max = i; }
}
// scale coeffs
pair<indexType,valueType> res(i_max,scale);
valueType s = scale;
for ( i = 1; i <= n; ++i, s *= scale ) ps[n-i] /= s;
ps[n] = 1;
return res;
}
//============================================================================
/*
.. divide a(x) by (x-r) so that a(x) = (x-r) * b(x)
*/
static
void
deflatePolynomial(
indexType n, // degree
valueType const a[],
valueType r,
valueType b[]
) {
// crossover index for forward/backward deflation
// G. Peters and J. H. Wilkinson.
// Practical problems arising in the solution of polynomial equations.
// J. Inst. Math. Appl. 8 (1971), 16–35.
indexType i_cross = 0;
valueType v_cross = std::abs(a[0]);
valueType r1 = r;
for ( indexType i = 1; i < n; ++i, r1 *= r ) {
valueType v_cross1 = std::abs(a[i]*r1);
if ( v_cross1 > v_cross )
{ v_cross = v_cross1; i_cross = i; }
}
b[n-1] = 1;
if ( isZero(a[n]-1) ) {
if ( i_cross > 0 ) {
b[0] = -a[0] / r;
for ( indexType j = 1; j < i_cross; ++j )
b[j] = (a[j]-b[j-1]) / r;
}
for ( indexType j = n-2; j >= i_cross; --j )
b[j] = a[j+1]+r*b[j+1];
} else {
valueType an = a[n];
if ( i_cross > 0 ) {
b[0] = -(a[0]/an) / r;
for ( indexType j = 1; j < i_cross; ++j )
b[j] = (a[j]/an-b[j-1]) / r;
}
for ( indexType j = n-2; j >= i_cross; --j )
b[j] = a[j+1]/an+r*b[j+1];
}
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// stable computation of polinomial
// p0 + p1*x + p2*x^2 + ... + pn*x^n
//
// (p0/x^n + p1/x^(n-1) + p2/(x^(n-2) + ... + pn)*x^n
//
#if 0
// UNUSED
valueType
CompHorner(
valueType const p[],
indexType Degree,
valueType x
) {
valueType xabs = std::abs(x);
bool reverse = xabs > 1;
//bool reverse = false;
if ( reverse ) x = valueType(1)/x;
indexType ii0 = reverse ? Degree : 0;
valueType res(p[ii0]);
valueType c = 0;
for ( indexType i = 1; i <= Degree; ++i ) {
indexType ii = reverse ? Degree-i : i;
valueType tmp, pi, sigma;
//TwoProduct( res, x, tmp, pi );
//TwoSum( tmp, p[ii], res, sigma );
res = res * x + p[ii];
//c = c * x + (pi+sigma);
}
//res += c;
if ( reverse ) res *= std::pow(x,Degree);
return res;
}
#endif
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
//
// x^3 + A x^2 + B x + C
static
inline
void
scaleCubicMonicPolynomial(
valueType A,
valueType B,
valueType C,
valueType & AS,
valueType & BS,
valueType & CS,
indexType & i_case,
valueType & scale
) {
valueType a = std::abs(A);
valueType b = std::sqrt(abs(B));
valueType c = std::cbrt(abs(C));
if ( a < b ) {
if ( b < c ) i_case = 0; // a < b < c --> c MAX
else i_case = 1; // a < b and c <= b --> b MAX
} else {
if ( a < c ) i_case = 0; // b <= a < c --> c MAX
else i_case = 2; // b <= a and c <= a --> a MAX
}
switch ( i_case ) {
case 0:
scale = c;
AS = A/c;
BS = (B/c)/c;
CS = C > 0 ? 1 : -1;
break;
case 1:
scale = b;
AS = A/b;
BS = B > 0 ? 1 : -1;
CS = ((C/b)/b)/b;
break;
case 2:
scale = a;
AS = A > 0 ? 1 : -1;
BS = (B/a)/a;
CS = ((C/a)/a)/a;
break;
}
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
//
// a3*x^3 + a2*x^2 + a1*x + a0 = (x-r)*(a3*x^2+b1*x+b0)
static
void
deflateCubicPolynomial(
valueType a3,
valueType a2,
valueType a1,
valueType a0,
valueType r,
valueType & b1,
valueType & b0
) {
indexType i_cross = 0;
valueType r2 = r*r;
valueType v_cross = std::abs(a0);
valueType v_cross1 = std::abs(a1*r);
if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 1; }
v_cross1 = std::abs(a2*r2);
if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 2; }
v_cross1 = std::abs(a3*r*r2);
if ( v_cross1 > v_cross ) i_cross = 3;
switch ( i_cross ) {
case 0: b1 = a2+a3*r; b0 = a1+r*b1; break;
case 1: b1 = a2+a3*r; b0 = -a0/r; break;
case 2:
case 3: b0 = -a0/r; b1 = (b0-a1)/r; break;
}
}
}
#endif
// EOF: PolynomialRoots-Utils.cc
