📄 rpoly.c
字号:
#include "v3p_netlib.h"
#undef abs
#undef min
#undef max
#include <math.h>
#define abs(x) ((x) >= 0 ? (x) : -(x))
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))
static void calcsc_(integer *type,
v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg);
static void fxshfr_(integer *l2, integer *nz,
v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg);
static void newest_(integer *type, doublereal *uu, doublereal *vv,
v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg);
static void nextk_(integer *type,
v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg);
static void quadit_(doublereal *uu, doublereal *vv, integer *nz,
v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg);
static void realit_(doublereal *sss, integer *nz, integer *iflag,
v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg);
static void quadsd_(integer *nn, doublereal *u, doublereal *v,
doublereal *p, doublereal *q, doublereal *a,
doublereal *b);
static void quad_(doublereal *a, doublereal *b1, doublereal *c,
doublereal *sr, doublereal *si, doublereal *lr,
doublereal *li);
#define global_1 (*v3p_netlib_rpoly_global_arg)
/* Table of constant values */
static doublereal c_b41 = 1.;
#ifdef _MSC_VER
// This needs to be before the start of the function that contains the offending code
# pragma warning ( disable : 4756)
#endif
/* ====================================================================== */
/* NIST Guide to Available Math Software. */
/* Fullsource for module 493 from package TOMS. */
/* Retrieved from NETLIB on Wed Jul 3 11:47:53 1996. */
/* ====================================================================== */
/* Subroutine */ void rpoly_(
doublereal* op, integer* degree,
doublereal* zeror, doublereal* zeroi,
logical* fail,
v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg
)
{
/* Builtin functions */
double log(doublereal), pow_di(doublereal *, integer *), exp(doublereal);
/* System generated locals */
integer i__1;
/* Local variables */
doublereal base;
doublereal temp[101];
real cosr, sinr;
integer i, j, l;
doublereal t;
real x, infin;
logical zerok;
doublereal aa, bb, cc;
real df, ff;
integer jj;
real sc, lo, dx, pt[101], xm;
integer nz;
doublereal factor;
real xx, yy, smalno;
integer nm1;
real bnd, min_, max_;
integer cnt;
real xxx;
/* FINDS THE ZEROS OF A REAL POLYNOMIAL */
/* OP - DOUBLE PRECISION VECTOR OF COEFFICIENTS IN */
/* ORDER OF DECREASING POWERS. */
/* DEGREE - INTEGER DEGREE OF POLYNOMIAL. */
/* ZEROR, ZEROI - OUTPUT DOUBLE PRECISION VECTORS OF */
/* REAL AND IMAGINARY PARTS OF THE */
/* ZEROS. */
/* FAIL - OUTPUT LOGICAL PARAMETER, TRUE ONLY IF */
/* LEADING COEFFICIENT IS ZERO OR IF RPOLY */
/* HAS FOUND FEWER THAN DEGREE ZEROS. */
/* IN THE LATTER CASE DEGREE IS RESET TO */
/* THE NUMBER OF ZEROS FOUND. */
/* TO CHANGE THE SIZE OF POLYNOMIALS WHICH CAN BE */
/* SOLVED, RESET THE DIMENSIONS OF THE ARRAYS IN THE */
/* COMMON AREA AND IN THE FOLLOWING DECLARATIONS. */
/* THE SUBROUTINE USES SINGLE PRECISION CALCULATIONS */
/* FOR SCALING, BOUNDS AND ERROR CALCULATIONS. ALL */
/* CALCULATIONS FOR THE ITERATIONS ARE DONE IN DOUBLE */
/* PRECISION. */
/* THE FOLLOWING STATEMENTS SET MACHINE CONSTANTS USED */
/* IN VARIOUS PARTS OF THE PROGRAM. 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.D0+ETA IS GREATER THAN 1. */
/* INFIN THE LARGEST FLOATING-POINT NUMBER. */
/* SMALNO THE SMALLEST POSITIVE FLOATING-POINT NUMBER */
/* IF THE EXPONENT RANGE DIFFERS IN SINGLE AND */
/* DOUBLE PRECISION THEN SMALNO AND INFIN */
/* SHOULD INDICATE THE SMALLER RANGE. */
/* BASE THE BASE OF THE FLOATING-POINT NUMBER */
/* SYSTEM USED. */
/* THE VALUES BELOW CORRESPOND TO THE BURROUGHS B6700 */
/* changed for sparc, but these seem better -- awf */
base = 2.0;
global_1.eta = 2.23e-16f;
infin = 3.40282346638528860e+38f;
smalno = 1e-33f;
/* ARE AND MRE REFER TO THE UNIT ERROR IN + AND * */
/* RESPECTIVELY. THEY ARE ASSUMED TO BE THE SAME AS */
/* ETA. */
global_1.are = global_1.eta;
global_1.mre = global_1.eta;
lo = smalno / global_1.eta;
/* INITIALIZATION OF CONSTANTS FOR SHIFT ROTATION */
xx = 0.70710678f;
yy = -xx;
cosr = -0.069756474f;
sinr = 0.99756405f;
*fail = FALSE_;
global_1.n = *degree;
global_1.nn = global_1.n + 1;
/* ALGORITHM FAILS IF THE LEADING COEFFICIENT IS ZERO. */
if (op[0] == 0.) {
*fail = TRUE_;
*degree = 0;
return;
}
/* REMOVE THE ZEROS AT THE ORIGIN IF ANY */
while (op[global_1.nn - 1] == 0.) {
j = *degree - global_1.n;
zeror[j] = 0.;
zeroi[j] = 0.;
--global_1.nn;
--global_1.n;
}
/* MAKE A COPY OF THE COEFFICIENTS */
for (i = 0; i < global_1.nn; ++i) {
global_1.p[i] = op[i];
}
/* START THE ALGORITHM FOR ONE ZERO */
L40:
if (global_1.n > 2) {
goto L60;
}
if (global_1.n < 1) {
return;
}
/* CALCULATE THE FINAL ZERO OR PAIR OF ZEROS */
if (global_1.n != 2) {
zeror[*degree-1] = -global_1.p[1] / global_1.p[0];
zeroi[*degree-1] = 0.;
return;
}
quad_(global_1.p, &global_1.p[1], &global_1.p[2],
&zeror[*degree-2], &zeroi[*degree-2], &zeror[*degree-1], &zeroi[*degree-1]);
return;
/* FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS. */
L60:
max_ = 0.0f;
min_ = infin;
for (i = 0; i < global_1.nn; ++i) {
x = (real)abs(global_1.p[i]);
if (x > max_) {
max_ = x;
}
if (x != 0.0f && x < min_) {
min_ = x;
}
}
/* SCALE IF THERE ARE LARGE OR VERY SMALL COEFFICIENTS */
/* COMPUTES 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 */
sc = lo / min_;
if (sc > 1.0f) {
goto L80;
}
if (max_ < 10.0f) {
goto L110;
}
if (sc == 0.0f) {
sc = smalno;
}
goto L90;
L80:
if (infin / sc < max_) {
goto L110;
}
L90:
l = (int)(log((doublereal)sc) / log(base) + 0.5);
factor = base;
factor = pow_di(&factor, &l);
if (factor == 1.) {
goto L110;
}
for (i = 0; i < global_1.nn; ++i) {
global_1.p[i] *= factor;
}
/* COMPUTE LOWER BOUND ON MODULI OF ZEROS. */
L110:
for (i = 0; i < global_1.nn; ++i) {
pt[i] = (real)abs(global_1.p[i]);
}
pt[global_1.nn - 1] = -pt[global_1.nn - 1];
/* COMPUTE UPPER ESTIMATE OF BOUND */
x = (real)exp((log(-pt[global_1.nn - 1]) - log(pt[0])) / global_1.n);
if (pt[global_1.n - 1] == 0.0f) {
goto L130;
}
/* IF NEWTON STEP AT THE ORIGIN IS BETTER, USE IT. */
xm = -pt[global_1.nn - 1] / pt[global_1.n - 1];
if (xm < x) {
x = xm;
}
/* CHOP THE INTERVAL (0,X) UNTIL FF .LE. 0 */
L130:
xm = x * 0.1f;
ff = pt[0];
for (i = 1; i < global_1.nn; ++i) {
ff = ff * xm + pt[i];
}
if (ff > 0.0f) {
x = xm;
goto L130;
}
dx = x;
/* DO NEWTON ITERATION UNTIL X CONVERGES TO TWO */
/* DECIMAL PLACES */
while (abs(dx/x) > 0.005f) {
ff = pt[0];
df = ff;
for (i = 1; i < global_1.n; ++i) {
ff = ff * x + pt[i];
df = df * x + ff;
}
ff = ff * x + pt[global_1.nn - 1];
dx = ff / df;
x -= dx;
}
bnd = x;
/* COMPUTE THE DERIVATIVE AS THE INITIAL K POLYNOMIAL */
/* AND DO 5 STEPS WITH NO SHIFT */
nm1 = global_1.n - 1;
for (i = 1; i < global_1.n; ++i) {
global_1.k[i] = (global_1.nn - i - 1) * global_1.p[i] / global_1.n;
}
global_1.k[0] = global_1.p[0];
aa = global_1.p[global_1.nn - 1];
bb = global_1.p[global_1.n - 1];
zerok = global_1.k[global_1.n - 1] == 0.;
for (jj = 1; jj <= 5; ++jj) {
cc = global_1.k[global_1.n - 1];
if (zerok) { /* USE UNSCALED FORM OF RECURRENCE */
for (i = 0; i < nm1; ++i) {
j = global_1.nn - i - 2;
global_1.k[j] = global_1.k[j - 1];
}
global_1.k[0] = 0.;
zerok = global_1.k[global_1.n - 1] == 0.;
}
else { /* USE SCALED FORM OF RECURRENCE */
t = -aa / cc;
for (i = 0; i < nm1; ++i) {
j = global_1.nn - i - 2;
global_1.k[j] = t * global_1.k[j - 1] + global_1.p[j];
}
global_1.k[0] = global_1.p[0];
zerok = abs(global_1.k[global_1.n - 1]) <= abs(bb) * global_1.eta * 10.0;
}
}
/* SAVE K FOR RESTARTS WITH NEW SHIFTS */
for (i = 0; i < global_1.n; ++i) {
temp[i] = global_1.k[i];
}
/* LOOP TO SELECT THE QUADRATIC CORRESPONDING TO EACH */
/* NEW SHIFT */
for (cnt = 1; cnt <= 20; ++cnt) {
/* QUADRATIC CORRESPONDS TO A DOUBLE SHIFT TO A */
/* NON-REAL POINT AND ITS COMPLEX CONJUGATE. THE POINT */
/* HAS MODULUS BND AND AMPLITUDE ROTATED BY 94 DEGREES */
/* FROM THE PREVIOUS SHIFT */
xxx = cosr * xx - sinr * yy;
yy = sinr * xx + cosr * yy;
xx = xxx;
global_1.sr = bnd * xx;
global_1.si = bnd * yy;
global_1.u = global_1.sr * -2.;
global_1.v = bnd;
/* SECOND STAGE CALCULATION, FIXED QUADRATIC */
i__1 = cnt * 20;
fxshfr_(&i__1, &nz, v3p_netlib_rpoly_global_arg);
if (nz == 0) {
goto L260;
}
/* THE SECOND STAGE JUMPS DIRECTLY TO ONE OF THE THIRD */
/* STAGE ITERATIONS AND RETURNS HERE IF SUCCESSFUL. */
/* DEFLATE THE POLYNOMIAL, STORE THE ZERO OR ZEROS AND */
/* RETURN TO THE MAIN ALGORITHM. */
j = *degree - global_1.n;
zeror[j] = global_1.szr;
zeroi[j] = global_1.szi;
global_1.nn -= nz;
global_1.n = global_1.nn - 1;
for (i = 0; i < global_1.nn; ++i) {
global_1.p[i] = global_1.qp[i];
}
if (nz == 1) {
goto L40;
}
zeror[j + 1] = global_1.lzr;
zeroi[j + 1] = global_1.lzi;
goto L40;
/* IF THE ITERATION IS UNSUCCESSFUL ANOTHER QUADRATIC */
/* IS CHOSEN AFTER RESTORING K */
L260:
for (i = 0; i < global_1.n; ++i) {
global_1.k[i] = temp[i];
}
}
/* RETURN WITH FAILURE IF NO CONVERGENCE WITH 20 */
/* SHIFTS */
*fail = TRUE_;
*degree -= global_1.n;
} /* rpoly_ */
/* Subroutine */
static void fxshfr_(integer *l2, integer *nz,
v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg)
{
/* Local variables */
integer type;
logical stry, vtry;
integer i, j, iflag;
doublereal s;
real betas, betav;
logical spass;
logical vpass;
doublereal ui, vi;
real ts, tv, vv;
real ots=0, otv=0, tss;
doublereal ss, oss, ovv, svu, svv;
real tvv;
/* COMPUTES UP TO L2 FIXED SHIFT K-POLYNOMIALS, */
/* TESTING FOR CONVERGENCE IN THE LINEAR OR QUADRATIC */
/* CASE. INITIATES ONE OF THE VARIABLE SHIFT */
/* ITERATIONS AND RETURNS WITH THE NUMBER OF ZEROS */
/* FOUND. */
/* L2 - LIMIT OF FIXED SHIFT STEPS */
/* NZ - NUMBER OF ZEROS FOUND */
*nz = 0;
betav = .25f;
betas = .25f;
oss = global_1.sr;
ovv = global_1.v;
/* EVALUATE POLYNOMIAL BY SYNTHETIC DIVISION */
quadsd_(&global_1.nn, &global_1.u, &global_1.v, global_1.p, global_1.qp, &global_1.a, &global_1.b);
calcsc_(&type, v3p_netlib_rpoly_global_arg);
for (j = 1; j <= *l2; ++j) {
/* CALCULATE NEXT K POLYNOMIAL AND ESTIMATE V */
nextk_(&type, v3p_netlib_rpoly_global_arg);
calcsc_(&type, v3p_netlib_rpoly_global_arg);
newest_(&type, &ui, &vi, v3p_netlib_rpoly_global_arg);
vv = (real)vi;
/* ESTIMATE S */
ss = 0.0;
if (global_1.k[global_1.n - 1] != 0.) {
ss = -global_1.p[global_1.nn - 1] / global_1.k[global_1.n - 1];
}
tv = 1.0f;
ts = 1.0f;
if (j == 1 || type == 3) {
goto L70;
}
/* COMPUTE RELATIVE MEASURES OF CONVERGENCE OF S AND V */
/* SEQUENCES */
if (vv != 0.0f) {
tv = (real)abs((vv - ovv) / vv);
}
if (ss != 0.0) {
ts = (real)abs((ss - oss) / ss);
}
/* IF DECREASING, MULTIPLY TWO MOST RECENT */
/* CONVERGENCE MEASURES */
tvv = 1.0f;
if (tv < otv) {
tvv = tv * otv;
}
tss = 1.0f;
if (ts < ots) {
tss = ts * ots;
}
/* COMPARE WITH CONVERGENCE CRITERIA */
vpass = tvv < betav;
spass = tss < betas;
if (! (spass || vpass)) {
goto L70;
}
/* AT LEAST ONE SEQUENCE HAS PASSED THE CONVERGENCE */
/* TEST. STORE VARIABLES BEFORE ITERATING */
svu = global_1.u;
svv = global_1.v;
for (i = 1; i <= global_1.n; ++i) {
global_1.svk[i - 1] = global_1.k[i - 1];
}
s = ss;
/* CHOOSE ITERATION ACCORDING TO THE FASTEST */
/* CONVERGING SEQUENCE */
vtry = FALSE_;
stry = FALSE_;
if (spass && (! vpass || tss < tvv)) {
goto L40;
}
L20:
quadit_(&ui, &vi, nz, v3p_netlib_rpoly_global_arg);
if (*nz > 0) {
return;
}
/* QUADRATIC ITERATION HAS FAILED. FLAG THAT IT HAS */
/* BEEN TRIED AND DECREASE THE CONVERGENCE CRITERION. */
vtry = TRUE_;
betav *= 0.25f;
/* TRY LINEAR ITERATION IF IT HAS NOT BEEN TRIED AND */
/* THE S SEQUENCE IS CONVERGING */
if (stry || ! spass) {
goto L50;
}
for (i = 1; i <= global_1.n; ++i) {
global_1.k[i - 1] = global_1.svk[i - 1];
}
L40:
realit_(&s, nz, &iflag, v3p_netlib_rpoly_global_arg);
if (*nz > 0) {
return;
}
/* LINEAR ITERATION HAS FAILED. FLAG THAT IT HAS BEEN */
/* TRIED AND DECREASE THE CONVERGENCE CRITERION */
stry = TRUE_;
betas *= 0.25f;
/* IF LINEAR ITERATION SIGNALS AN ALMOST DOUBLE REAL */
/* ZERO ATTEMPT QUADRATIC ITERATION */
if (iflag != 0) {
ui = -(s + s);
vi = s * s;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -