我有一个类似于这里的多项式类：Polynomial.java。除了我还没有弄清楚如何找到多项式的零点。我有Jenkins-Traub Algorithm的负责人，但我不知道如何在 Java 中实现它。我确实设法找到了 FORTRAN 版本，但我没有使用 FORTRAN 的经验。是否有演示该算法的开源项目？也许有人可以在这里用伪代码写出来？
问问题
1039 次
3 回答
2
我已将 Jenkins–Traub 算法的 Fortran 源代码移植到 Java（实系数版本）。它是我的Java DSP 集合的一部分。
于 2013-08-25T21:38:57.667 回答
0
尝试数值方法。他们有一个实现。
于 2012-02-29T07:20:40.697 回答
0
在我求解多项式的过程中，我达到了 5 点：你不能合理地求解 5 次或更高阶的方程。因此，您必须转向数值分析（或其他近似方法），因为您正在离开绝对领域，并且您正在进入“让我们进入这里，也许有炸弹”的情况。
甚至 MathWorld 也将其描述为带有明显嘲笑的“复杂方法”，而其他人则称其为黑盒方法中最好的（这意味着您不知道它是如何工作的，但您认为它是有效的）。
我将在这里发布我的版本，该版本已经过测试并且与小度数的精确求解器（例如四次）相同。
（由// http://hvks.com/Numerical/Downloads/cpoly_ver_1.cpp提供，从 Fortran 翻译成 C）
public class Poly {
double DBL_EPSILON=2.2204460492503131e-16, DBL_RADIX=2, DBL_MAX=1E+37, DBL_MIN=1E-37;
double sr, si, are, mre, eta, infin;
double[] pvri=new double[2], tri=new double[2];
int nn;
double[] pr, pi, hr, hi, qpr, qpi, qhr, qhi, shr, shi;
public int poly(double[] opr, double[] opi, int degree, double[] zeror, double[] zeroi) {
int cnt1, cnt2, idnn2, i, conv=0;
double xx, yy, cosr, sinr, smalno, base, xxx, bnd;
double[] zri=new double[2];
base = DBL_RADIX;
eta = DBL_EPSILON;
infin = DBL_MAX;
smalno = DBL_MIN;
are = eta;
mre = 2.0 * Math.sqrt( 2.0 ) * eta;
xx = 0.70710678;
yy = -xx;
cosr = -0.060756474;
sinr = -0.99756405;
nn = degree;
// Algorithm fails if the leading coefficient is zero
if( opr[ 0 ] == 0 && opi[ 0 ] == 0 )
return -1;
// Allocate arrays
pr = new double [ degree+1 ];
pi = new double [ degree+1 ];
hr = new double [ degree+1 ];
hi = new double [ degree+1 ];
qpr= new double [ degree+1 ];
qpi= new double [ degree+1 ];
qhr= new double [ degree+1 ];
qhi= new double [ degree+1 ];
shr= new double [ degree+1 ];
shi= new double [ degree+1 ];
// Remove the zeros at the origin if any
while( opr[ nn ] == 0 && opi[ nn ] == 0 ) {
idnn2 = degree - nn;
zeror[ idnn2 ] = 0;
zeroi[ idnn2 ] = 0;
nn--;
}
// Make a copy of the coefficients
for( i = 0; i <= nn; i++ ) {
pr[ i ] = opr[ i ];
pi[ i ] = opi[ i ];
shr[ i ] = cmod( pr[ i ], pi[ i ] );
}
// Scale the polynomial
bnd = scale( nn, shr, eta, infin, smalno, base );
if( bnd != 1 ) {
for( i = 0; i <= nn; i++ ) {
pr[ i ] *= bnd;
pi[ i ] *= bnd;
}
}
do {
if( nn <= 1 ) {
double[] cri=cdivid( -pr[ 1 ], -pi[ 1 ], new double[]{pr[ 0 ], pi[ 0 ]});
zeror[ degree-1 ]=cri[0]; zeroi[ degree-1 ]=cri[1];
break;
}
// Calculate bnd, alower bound on the modulus of the zeros
for( i = 0; i<= nn; i++ )
shr[ i ] = cmod( pr[ i ], pi[ i ] );
bnd=cauchy( nn, shr, shi);
// Outer loop to control 2 Major passes with different sequences of shifts
for( cnt1 = 1; cnt1 <= 2; cnt1++ ) {
// First stage calculation , no shift
noshft( 5 );
// Inner loop to select a shift
for( cnt2 = 1; cnt2 <= 9; cnt2++ ) {
// Shift is chosen with modulus bnd and amplitude rotated by 94 degree from the previous shif
xxx = cosr * xx - sinr * yy;
yy = sinr * xx + cosr * yy;
xx = xxx;
sr = bnd * xx;
si = bnd * yy;
// Second stage calculation, fixed shift
conv=fxshft( 10 * cnt2, zri);
if( conv==1 ) {
// The second stage jumps directly to the third stage ieration
// If successful the zero is stored and the polynomial deflated
idnn2 = degree - nn;
zeror[ idnn2 ] = zri[0];
zeroi[ idnn2 ] = zri[1];
nn--;
for( i = 0; i <= nn; i++ ) {
pr[ i ] = qpr[ i ];
pi[ i ] = qpi[ i ];
}
break;
}
// If the iteration is unsuccessful another shift is chosen
}
// if 9 shifts fail, the outer loop is repeated with another sequence of shifts
}
}
while(true);
// The zerofinder has failed on two major passes
// return empty handed with the number of roots found (less than the original degree)
degree -= nn;
return degree;
}
// COMPUTES THE DERIVATIVE POLYNOMIAL AS THE INITIAL H
// POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS.
//
public void noshft(int l1 ) {
int i, j, jj, n, nm1;
double xni, t1, t2;
n = nn;
nm1 = n - 1;
for( i = 0; i < n; i++ ) {
xni = nn - i;
hr[ i ] = xni * pr[ i ] / n;
hi[ i ] = xni * pi[ i ] / n;
}
for( jj = 1; jj <= l1; jj++ ) {
if( cmod( hr[ n - 1 ], hi[ n - 1 ] ) > eta * 10 * cmod( pr[ n - 1 ], pi[ n - 1 ] ) ) {
double[] tri2=cdivid( -pr[ nn ], -pi[ nn ], new double[]{hr[ n - 1 ], hi[ n - 1 ]});
tri[0]=tri2[0]; tri[1]=tri2[1];
for( i = 0; i < nm1; i++ ) {
j = nn - i - 1;
t1 = hr[ j - 1 ];
t2 = hi[ j - 1 ];
hr[ j ] = tri[0] * t1 - tri[1] * t2 + pr[ j ];
hi[ j ] = tri[0] * t2 + tri[1] * t1 + pi[ j ];
}
hr[ 0 ] = pr[ 0 ];
hi[ 0 ] = pi[ 0 ];
}
else {
// If the constant term is essentially zero, shift H coefficients
for( i = 0; i < nm1; i++ ) {
j = nn - i - 1;
hr[ j ] = hr[ j - 1 ];
hi[ j ] = hi[ j - 1 ];
}
hr[ 0 ] = 0;
hi[ 0 ] = 0;
}
}
}
// COMPUTES L2 FIXED-SHIFT H POLYNOMIALS AND TESTS FOR CONVERGENCE.
// INITIATES A VARIABLE-SHIFT ITERATION AND RETURNS WITH THE
// APPROXIMATE ZERO IF SUCCESSFUL.
// L2 - LIMIT OF FIXED SHIFT STEPS
// ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE.
// CONV - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION
//
public int fxshft(int l2, double[] zri) {
int i, j, n, conv=0;
int test, pasd;
double otr=0, oti=0, svsr, svsi;
n = nn;
polyev( nn, sr, si, pr, pi, qpr, qpi, pvri);
test = 1;
pasd = 0;
// Calculate first T = -P(S)/H(S)
int bool=calct();
// Main loop for second stage
for( j = 1; j <= l2; j++ ) {
otr = tri[0];
oti = tri[1];
// Compute the next H Polynomial and new t
nexth( bool );
bool=calct();
zri[0] = sr + tri[0];
zri[1] = si + tri[1];
// Test for convergence unless stage 3 has failed once or this
// is the last H Polynomial
if( !( bool==1 || test==0 || j == 12 ) ) {
if( cmod( tri[0] - otr, tri[1] - oti ) < 0.5 * cmod( zri[0], zri[1] ) ) {
if( pasd==1 ) {
// The weak convergence test has been passwed twice, start the third stage
// Iteration, after saving the current H polynomial and shift
for( i = 0; i < n; i++ ) {
shr[ i ] = hr[ i ];
shi[ i ] = hi[ i ];
}
svsr = sr;
svsi = si;
conv=vrshft( 10, zri);
if(conv==1) return conv;
//The iteration failed to converge. Turn off testing and restore h,s,pv and T
test = 0;
for( i = 0; i < n; i++ ) {
hr[ i ] = shr[ i ];
hi[ i ] = shi[ i ];
}
sr = svsr;
si = svsi;
polyev( nn, sr, si, pr, pi, qpr, qpi, pvri);
bool=calct();
continue;
}
pasd = 1;
}
else
pasd = 0;
}
}
// Attempt an iteration with final H polynomial from second stage
conv=vrshft( 10, zri);
return conv;
}
// CARRIES OUT THE THIRD STAGE ITERATION.
// L3 - LIMIT OF STEPS IN STAGE 3.
// ZR,ZI - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE
// ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE ON EXIT.
// CONV - .TRUE. IF ITERATION CONVERGES
//
public int vrshft(int l3, double[] zri) {
int b, bool=0, conv=0;
int i, j;
double mp=0, ms, omp=0, relstp=0, r1, r2, tp;
conv = 0;
b = 0;
sr = zri[0];
si = zri[1];
// Main loop for stage three
for( i = 1; i <= l3; i++ ) {
// Evaluate P at S and test for convergence
polyev( nn, sr, si, pr, pi, qpr, qpi, pvri);
mp = cmod( pvri[0], pvri[1] );
ms = cmod( sr, si );
if( mp <= 20 * errev( nn, qpr, qpi, ms, mp, are, mre ) ) {
// Polynomial value is smaller in value than a bound onthe error
// in evaluationg P, terminate the ietartion
conv = 1;
zri[0] = sr;
zri[1] = si;
return conv;
}
if( i != 1 ) {
if( !( b==1 || mp < omp || relstp >= 0.05 ) ) {
// Iteration has stalled. Probably a cluster of zeros. Do 5 fixed
// shift steps into the cluster to force one zero to dominate
tp = relstp;
b = 1;
if( relstp < eta ) tp = eta;
r1 = Math.sqrt( tp );
r2 = sr * ( 1 + r1 ) - si * r1;
si = sr * r1 + si * ( 1 + r1 );
sr = r2;
polyev( nn, sr, si, pr, pi, qpr, qpi, pvri);
for( j = 1; j <= 5; j++ ) {
bool=calct();
nexth( bool );
}
omp = infin;
bool=calct();
nexth( bool );
bool=calct();
if( bool==0 ) {
relstp = cmod( tri[0], tri[1] ) / cmod( sr, si );
sr += tri[0];
si += tri[1];
}
continue;
}
// Exit if polynomial value increase significantly
if( mp *0.1 > omp ) return conv;
}
omp = mp;
}
return conv;
}
// COMPUTES T = -P(S)/H(S).
// BOOL - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO.
public int calct() {
int n;
double[] hvri=new double[2];
n = nn;
// evaluate h(s)
polyev( n - 1, sr, si, hr, hi, qhr, qhi, hvri);
int bool = cmod( hvri[0], hvri[1] ) <= are * 10 * cmod( hr[ n - 1 ], hi[ n - 1 ] ) ? 1 : 0;
if( bool==0 ) {
double[] tri2=cdivid( -pvri[0], -pvri[1], new double[]{hvri[0], hvri[1]});
tri[0]=tri2[0]; tri[1]=tri2[1];
return bool;
}
tri[0] = 0;
tri[1] = 0;
return bool;
}
// CALCULATES THE NEXT SHIFTED H POLYNOMIAL.
// BOOL - LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO
//
public void nexth(int bool ) {
int j, n;
double t1, t2;
n = nn;
if( bool==0 ) {
for( j = 1; j < n; j++ ) {
t1 = qhr[ j - 1 ];
t2 = qhi[ j - 1 ];
hr[ j ] = tri[0] * t1 - tri[1] * t2 + qpr[ j ];
hi[ j ] = tri[0] * t2 + tri[1] * t1 + qpi[ j ];
}
hr[ 0 ] = qpr[ 0 ];
hi[ 0 ] = qpi[ 0 ];
return;
}
// If h[s] is zero replace H with qh
for( j = 1; j < n; j++ ) {
hr[ j ] = qhr[ j - 1 ];
hi[ j ] = qhi[ j - 1 ];
}
hr[ 0 ] = 0;
hi[ 0 ] = 0;
}
// EVALUATES A POLYNOMIAL P AT S BY THE HORNER RECURRENCE
// PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV.
//
public void polyev(int nn, double sr, double si, double pr[], double pi[], double qr[], double qi[], double[] pvri) {
int i;
double t;
qr[ 0 ] = pr[ 0 ];
qi[ 0 ] = pi[ 0 ];
pvri[0] = qr[ 0 ];
pvri[1] = qi[ 0 ];
for( i = 1; i <= nn; i++ ) {
t = ( pvri[0] ) * sr - ( pvri[1] ) * si + pr[ i ];
pvri[1] = ( pvri[0] ) * si + ( pvri[1] ) * sr + pi[ i ];
pvri[0] = t;
qr[ i ] = pvri[0];
qi[ i ] = pvri[1];
}
}
// BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER RECURRENCE.
// QR,QI - THE PARTIAL SUMS
// MS -MODULUS OF THE POINT
// MP -MODULUS OF POLYNOMIAL VALUE
// ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION
//
public double errev(int nn, double qr[], double qi[], double ms, double mp, double are, double mre ) {
int i;
double e;
e = cmod( qr[ 0 ], qi[ 0 ] ) * mre / ( are + mre );
for( i = 0; i <= nn; i++ )
e = e * ms + cmod( qr[ i ], qi[ i ] );
return e * ( are + mre ) - mp * mre;
}
// CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A
// POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS.
//
public double cauchy(int nn, double[] pt, double[] q) {
int i, n;
double x, xm, f, dx, df;
pt[ nn ] = -pt[ nn ];
// Compute upper estimate bound
n = nn;
x = Math.exp( Math.log( -pt[ nn ] ) - Math.log( pt[ 0 ] ) ) / n;
if( pt[ n - 1 ] != 0 ) {
// Newton step at the origin is better, use it
xm = -pt[ nn ] / pt[ n - 1 ];
if( xm < x ) x = xm;
}
// Chop the interval (0,x) until f < 0
while(true) {
xm = x * 0.1;
f = pt[ 0 ];
for( i = 1; i <= nn; i++ )
f = f * xm + pt[ i ];
if( f <= 0 )
break;
x = xm;
}
dx = x;
// Do Newton iteration until x converges to two decimal places
while( Math.abs( dx / x ) > 0.005 ) {
q[ 0 ] = pt[ 0 ];
for( i = 1; i <= nn; i++ )
q[ i ] = q[ i - 1 ] * x + pt[ i ];
f = q[ nn ];
df = q[ 0 ];
for( i = 1; i < n; i++ )
df = df * x + q[ i ];
dx = f / df;
x -= dx;
}
return x;
}
// RETURNS A SCALE FACTOR TO MULTIPLY THE COEFFICIENTS OF THE POLYNOMIAL.
// THE SCALING IS DONE TO AVOID OVERFLOW AND TO AVOID UNDETECTED UNDERFLOW
// INTERFERING WITH THE CONVERGENCE CRITERION. THE FACTOR IS A POWER OF THE
// BASE.
// PT - MODULUS OF COEFFICIENTS OF P
// ETA, INFIN, SMALNO, BASE - CONSTANTS DESCRIBING THE FLOATING POINT ARITHMETIC.
//
public double scale(int nn, double pt[], double eta, double infin, double smalno, double base ) {
int i, l;
double hi, lo, max, min, x, sc;
double fn_val;
// Find largest and smallest moduli of coefficients
hi = Math.sqrt( infin );
lo = smalno / eta;
max = 0;
min = infin;
for( i = 0; i <= nn; i++ ) {
x = pt[ i ];
if( x > max ) max = x;
if( x != 0 && x < min ) min = x;
}
// Scale only if there are very large or very small components
fn_val = 1;
if( min >= lo && max <= hi ) return fn_val;
x = lo / min;
if( x <= 1 )
sc = 1 / ( Math.sqrt( max )* Math.sqrt( min ) );
else {
sc = x;
if( infin / sc > max ) sc = 1;
}
l = (int)( Math.log( sc ) / Math.log(base ) + 0.5 );
fn_val = Math.pow( base , l );
return fn_val;
}
// COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW.
//
public double[] cdivid(double ar, double ai, double[] bri) {
double r, d, t, infin;
double[] cri=new double[2];
if( bri[0] == 0 && bri[1] == 0 ) {
// Division by zero, c = infinity
double[] mc=mcon();
t=mc[1]; infin=mc[2]; t=mc[3]; t=mc[0];
cri[0] = infin;
cri[1] = infin;
return cri;
}
if( Math.abs( bri[0] ) < Math.abs( bri[1] ) ) {
r = bri[0]/ bri[1];
d = bri[1] + r * bri[0];
cri[0] = ( ar * r + ai ) / d;
cri[1] = ( ai * r - ar ) / d;
return cri;
}
r = bri[1] / bri[0];
d = bri[0] + r * bri[1];
cri[0] = ( ar + ai * r ) / d;
cri[1] = ( ai - ar * r ) / d;
return cri;
}
// MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW.
//
public double cmod( double r, double i ) {
double ar, ai;
ar = Math.abs( r );
ai = Math.abs( i );
if( ar < ai )
return ai * Math.sqrt( 1.0 + Math.pow( ( ar / ai ), 2.0 ) );
if( ar > ai )
return ar * Math.sqrt( 1.0 + Math.pow( ( ai / ar ), 2.0 ) );
return ar * Math.sqrt( 2.0 );
}
// MCON PROVIDES MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE PROGRAM.
// THE USER MAY EITHER SET THEM DIRECTLY OR USE THE STATEMENTS BELOW TO
// COMPUTE THEM. THE MEANING OF THE FOUR CONSTANTS ARE -
// ETA THE MAXIMUM RELATIVE REPRESENTATION ERROR WHICH CAN BE DESCRIBED
// AS THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
// 1.0_dp + ETA > 1.0.
// INFINY THE LARGEST FLOATING-POINT NUMBER
// SMALNO THE SMALLEST POSITIVE FLOATING-POINT NUMBER
// BASE THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED
//
public double[] mcon() {
return new double[]{DBL_RADIX, Double.MIN_VALUE, Double.MAX_VALUE, Double.MIN_VALUE};
}
}
于 2019-11-15T13:49:52.190 回答