rpoly.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

            goto L20;
        }
/* RESTORE VARIABLES */
L50:
        global_1.u = svu;
        global_1.v = svv;
        for (i = 1; i <= global_1.n; ++i) {
            global_1.k[i - 1] = global_1.svk[i - 1];
        }
/* TRY QUADRATIC ITERATION IF IT HAS NOT BEEN TRIED */
/* AND THE V SEQUENCE IS CONVERGING */
        if (vpass && ! vtry) {
            goto L20;
        }
/* RECOMPUTE QP AND SCALAR VALUES TO CONTINUE THE */
/* SECOND STAGE */
        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);
L70:
        ovv = vv;
        oss = ss;
        otv = tv;
        ots = ts;
    }
} /* fxshfr_ */

/* Subroutine */
static void quadit_(doublereal *uu, doublereal *vv, integer *nz,
                    v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg)
{
    /* Local variables */
    integer type, i, j;
    doublereal t;
    logical tried;
    real ee;
    doublereal ui, vi;
    real mp, zm;
    real relstp=0, omp=0;

/* VARIABLE-SHIFT K-POLYNOMIAL ITERATION FOR A */
/* QUADRATIC FACTOR CONVERGES ONLY IF THE ZEROS ARE */
/* EQUIMODULAR OR NEARLY SO. */
/* UU,VV - COEFFICIENTS OF STARTING QUADRATIC */
/* NZ - NUMBER OF ZERO FOUND */
    *nz = 0;
    tried = FALSE_;
    global_1.u = *uu;
    global_1.v = *vv;
    j = 0;
/* MAIN LOOP */
L10:
    quad_(&c_b41, &global_1.u, &global_1.v, &global_1.szr, &global_1.szi, &global_1.lzr, &global_1.lzi);
/* RETURN IF ROOTS OF THE QUADRATIC ARE REAL AND NOT */
/* CLOSE TO MULTIPLE OR NEARLY EQUAL AND  OF OPPOSITE */
/* SIGN */
    if (abs(abs(global_1.szr) - abs(global_1.lzr)) > abs(global_1.lzr) * .01) {
        return;
    }
/* EVALUATE POLYNOMIAL BY QUADRATIC SYNTHETIC DIVISION */
    quadsd_(&global_1.nn, &global_1.u, &global_1.v, global_1.p, global_1.qp, &global_1.a, &global_1.b);
    mp = (real)abs(global_1.a - global_1.szr * global_1.b)
       + (real)abs(global_1.szi * global_1.b);
/* COMPUTE A RIGOROUS  BOUND ON THE ROUNDING ERROR IN */
/* EVALUTING P */
    zm = (real)sqrt(abs(global_1.v));
    ee = (real)abs(global_1.qp[0]) * 2.0f;
    t = -global_1.szr * global_1.b;
    for (i = 2; i <= global_1.n; ++i) {
        ee = ee * zm + (real)abs(global_1.qp[i - 1]);
    }
    ee = ee * zm + (real)abs(global_1.a + t);
    ee = (real)((global_1.mre * 5.0 + global_1.are * 4.0) * ee
       - (global_1.mre * 5.0 + global_1.are * 2.0) * (abs(global_1.a + t) + abs(global_1.b) * zm)
       + global_1.are * 2.0 * abs(t));
/* ITERATION HAS CONVERGED SUFFICIENTLY IF THE */
/* POLYNOMIAL VALUE IS LESS THAN 20 TIMES THIS BOUND */
    if (mp <= ee * 20.0f) {
        *nz = 2;
        return;
    }
/* STOP ITERATION AFTER 20 STEPS */
    if (++j > 20) {
        return;
    }
    if (j < 2) {
        goto L50;
    }
    if (relstp > 0.01f || mp < omp || tried) {
        goto L50;
    }
/* A CLUSTER APPEARS TO BE STALLING THE CONVERGENCE. */
/* FIVE FIXED SHIFT STEPS ARE TAKEN WITH A U,V CLOSE */
/* TO THE CLUSTER */
    if (relstp < global_1.eta) {
        relstp = global_1.eta;
    }
    relstp = (float)sqrt(relstp);
    global_1.u -= global_1.u * relstp;
    global_1.v += global_1.v * relstp;
    quadsd_(&global_1.nn, &global_1.u, &global_1.v, global_1.p, global_1.qp, &global_1.a, &global_1.b);
    for (i = 1; i <= 5; ++i) {
        calcsc_(&type, v3p_netlib_rpoly_global_arg);
        nextk_(&type, v3p_netlib_rpoly_global_arg);
    }
    tried = TRUE_;
    j = 0;
L50:
    omp = mp;
/* CALCULATE NEXT K POLYNOMIAL AND NEW U AND V */
    calcsc_(&type, v3p_netlib_rpoly_global_arg);
    nextk_(&type, v3p_netlib_rpoly_global_arg);
    calcsc_(&type, v3p_netlib_rpoly_global_arg);
    newest_(&type, &ui, &vi, v3p_netlib_rpoly_global_arg);
/* IF VI IS ZERO THE ITERATION IS NOT CONVERGING */
    if (vi == 0.) {
        return;
    }
    relstp = (real)abs((vi - global_1.v) / vi);
    global_1.u = ui;
    global_1.v = vi;
    goto L10;
} /* quadit_ */

/* Subroutine */
static void realit_(doublereal *sss, integer *nz, integer *iflag,
                    v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg)
{
    /* Local variables */
    integer i, j;
    doublereal s, t=0;
    real ee, mp, ms;
    doublereal kv, pv;
    real omp=0;

/* VARIABLE-SHIFT H POLYNOMIAL ITERATION FOR A REAL */
/* ZERO. */
/* SSS   - STARTING ITERATE */
/* NZ    - NUMBER OF ZERO FOUND */
/* IFLAG - FLAG TO INDICATE A PAIR OF ZEROS NEAR REAL */
/*         AXIS. */
    *nz = 0;
    s = *sss;
    *iflag = 0;
    j = 0;
/* MAIN LOOP */
L10:
    pv = global_1.p[0];
/* EVALUATE P AT S */
    global_1.qp[0] = pv;
    for (i = 2; i <= global_1.nn; ++i) {
        pv = pv * s + global_1.p[i - 1];
        global_1.qp[i - 1] = pv;
    }
    mp = (real)abs(pv);
/* COMPUTE A RIGOROUS BOUND ON THE ERROR IN EVALUATING */
/* P */
    ms = (real)abs(s);
    ee = (real)(global_1.mre / (global_1.are + global_1.mre) * abs(global_1.qp[0]));
    for (i = 2; i <= global_1.nn; ++i) {
        ee = ee * ms + (real)abs(global_1.qp[i - 1]);
    }
/* ITERATION HAS CONVERGED SUFFICIENTLY IF THE */
/* POLYNOMIAL VALUE IS LESS THAN 20 TIMES THIS BOUND */
    if (mp <= ((global_1.are + global_1.mre) * ee - global_1.mre * mp) * 20.0f) {
        *nz = 1;
        global_1.szr = s;
        global_1.szi = 0.;
        return;
    }
/* STOP ITERATION AFTER 10 STEPS */
    if (++j > 10) {
        return;
    }
    if (j < 2) {
        goto L50;
    }
    if (abs(t) > abs(s - t) * 0.001 || mp <= omp) {
        goto L50;
    }
/* A CLUSTER OF ZEROS NEAR THE REAL AXIS HAS BEEN */
/* ENCOUNTERED RETURN WITH IFLAG SET TO INITIATE A */
/* QUADRATIC ITERATION */
    *iflag = 1;
    *sss = s;
    return;
/* RETURN IF THE POLYNOMIAL VALUE HAS INCREASED */
/* SIGNIFICANTLY */
L50:
    omp = mp;
/* COMPUTE T, THE NEXT POLYNOMIAL, AND THE NEW ITERATE */
    kv = global_1.k[0];
    global_1.qk[0] = kv;
    for (i = 2; i <= global_1.n; ++i) {
        kv = kv * s + global_1.k[i - 1];
        global_1.qk[i - 1] = kv;
    }
    if (abs(kv) <= abs(global_1.k[global_1.n - 1]) * 10.0 * global_1.eta) {
        goto L80;
    }
/* USE THE SCALED FORM OF THE RECURRENCE IF THE VALUE */
/* OF K AT S IS NONZERO */
    t = -pv / kv;
    global_1.k[0] = global_1.qp[0];
    for (i = 2; i <= global_1.n; ++i) {
        global_1.k[i - 1] = t * global_1.qk[i - 2] + global_1.qp[i - 1];
    }
    goto L100;
/* USE UNSCALED FORM */
L80:
    global_1.k[0] = 0.;
    for (i = 2; i <= global_1.n; ++i) {
        global_1.k[i - 1] = global_1.qk[i - 2];
    }
L100:
    kv = global_1.k[0];
    for (i = 2; i <= global_1.n; ++i) {
        kv = kv * s + global_1.k[i - 1];
    }
    t = 0.;
    if (abs(kv) > abs(global_1.k[global_1.n - 1]) * 10.0 * global_1.eta) {
        t = -pv / kv;
    }
    s += t;
    goto L10;
} /* realit_ */

/* Subroutine */
static void calcsc_(integer *type,
                    v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg)
{
/* THIS ROUTINE CALCULATES SCALAR QUANTITIES USED TO */
/* COMPUTE THE NEXT K POLYNOMIAL AND NEW ESTIMATES OF */
/* THE QUADRATIC COEFFICIENTS. */
/* TYPE - INTEGER VARIABLE SET HERE INDICATING HOW THE */
/* CALCULATIONS ARE NORMALIZED TO AVOID OVERFLOW */
/* SYNTHETIC DIVISION OF K BY THE QUADRATIC 1,U,V */
    quadsd_(&global_1.n, &global_1.u, &global_1.v, global_1.k, global_1.qk, &global_1.c, &global_1.d);
    if (abs(global_1.c) > abs(global_1.k[global_1.n - 1]) * 100.0 * global_1.eta) {
        goto L10;
    }
    if (abs(global_1.d) > abs(global_1.k[global_1.n - 2]) * 100.0 * global_1.eta) {
        goto L10;
    }
    *type = 3;
/* TYPE=3 INDICATES THE QUADRATIC IS ALMOST A FACTOR */
/* OF K */
    return;
L10:
    if (abs(global_1.d) < abs(global_1.c)) {
        goto L20;
    }
    *type = 2;
/* TYPE=2 INDICATES THAT ALL FORMULAS ARE DIVIDED BY D */
    global_1.e = global_1.a / global_1.d;
    global_1.f = global_1.c / global_1.d;
    global_1.g = global_1.u * global_1.b;
    global_1.h = global_1.v * global_1.b;
    global_1.a3 = (global_1.a + global_1.g) * global_1.e + global_1.h * (global_1.b / global_1.d);
    global_1.a1 = global_1.b * global_1.f - global_1.a;
    global_1.a7 = (global_1.f + global_1.u) * global_1.a + global_1.h;
    return;
L20:
    *type = 1;
/* TYPE=1 INDICATES THAT ALL FORMULAS ARE DIVIDED BY C */
    global_1.e = global_1.a / global_1.c;
    global_1.f = global_1.d / global_1.c;
    global_1.g = global_1.u * global_1.e;
    global_1.h = global_1.v * global_1.b;
    global_1.a3 = global_1.a * global_1.e + (global_1.h / global_1.c + global_1.g) * global_1.b;
    global_1.a1 = global_1.b - global_1.a * (global_1.d / global_1.c);
    global_1.a7 = global_1.a + global_1.g * global_1.d + global_1.h * global_1.f;
    return;
} /* calcsc_ */

/* Subroutine */
static void nextk_(integer *type,
                    v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg)
{
    /* Local variables */
    doublereal temp;
    integer i;

/* COMPUTES THE NEXT K POLYNOMIALS USING SCALARS */
/* COMPUTED IN CALCSC */
    if (*type == 3) {
        goto L40;
    }
    temp = global_1.a;
    if (*type == 1) {
        temp = global_1.b;
    }
    if (abs(global_1.a1) > abs(temp) * global_1.eta * 10.0) {
        goto L20;
    }
/* IF A1 IS NEARLY ZERO THEN USE A SPECIAL FORM OF THE */
/* RECURRENCE */
    global_1.k[0] = 0.;
    global_1.k[1] = -global_1.a7 * global_1.qp[0];
    for (i = 3; i <= global_1.n; ++i) {
        global_1.k[i - 1] = global_1.a3 * global_1.qk[i - 3] - global_1.a7 * global_1.qp[i - 2];
    }
    return;
/* USE SCALED FORM OF THE RECURRENCE */
L20:
    global_1.a7 /= global_1.a1;
    global_1.a3 /= global_1.a1;
    global_1.k[0] = global_1.qp[0];
    global_1.k[1] = global_1.qp[1] - global_1.a7 * global_1.qp[0];
    for (i = 3; i <= global_1.n; ++i) {
        global_1.k[i - 1] = global_1.a3 * global_1.qk[i - 3] - global_1.a7 * global_1.qp[i - 2] + global_1.qp[i - 1];
    }
    return;
/* USE UNSCALED FORM OF THE RECURRENCE IF TYPE IS 3 */
L40:
    global_1.k[0] = 0.;
    global_1.k[1] = 0.;
    for (i = 3; i <= global_1.n; ++i) {
        global_1.k[i - 1] = global_1.qk[i - 3];
    }
} /* nextk_ */

/* Subroutine */
static void newest_(integer *type, doublereal *uu, doublereal *vv,
                    v3p_netlib_rpoly_global_t* v3p_netlib_rpoly_global_arg)
{
    doublereal temp, a4, a5, b1, b2, c1, c2, c3, c4;

/* COMPUTE NEW ESTIMATES OF THE QUADRATIC COEFFICIENTS */
/* USING THE SCALARS COMPUTED IN CALCSC. */
/* USE FORMULAS APPROPRIATE TO SETTING OF TYPE. */
    if (*type == 3) {
        goto L30;
    }
    if (*type == 2) {
        goto L10;
    }
    a4 = global_1.a + global_1.u * global_1.b + global_1.h * global_1.f;
    a5 = global_1.c + (global_1.u + global_1.v * global_1.f) * global_1.d;
    goto L20;
L10:
    a4 = (global_1.a + global_1.g) * global_1.f + global_1.h;
    a5 = (global_1.f + global_1.u) * global_1.c + global_1.v * global_1.d;
/* EVALUATE NEW QUADRATIC COEFFICIENTS. */
L20:
    b1 = -global_1.k[global_1.n - 1] / global_1.p[global_1.nn - 1];
    b2 = -(global_1.k[global_1.n - 2] + b1 * global_1.p[global_1.n - 1]) / global_1.p[global_1.nn - 1];
    c1 = global_1.v * b2 * global_1.a1;
    c2 = b1 * global_1.a7;
    c3 = b1 * b1 * global_1.a3;
    c4 = c1 - c2 - c3;
    temp = a5 + b1 * a4 - c4;
    if (temp == 0.) {
        goto L30;
    }
    *uu = global_1.u - (global_1.u * (c3 + c2) + global_1.v * (b1 * global_1.a1 + b2 * global_1.a7)) / temp;
    *vv = global_1.v * (c4 / temp + 1.0);
    return;
/* IF TYPE=3 THE QUADRATIC IS ZEROED */
L30:
    *uu = 0.;
    *vv = 0.;
    return;
} /* newest_ */

/* Subroutine */
static void quadsd_(
  integer *nn, doublereal *u, doublereal *v,
  doublereal *p, doublereal *q, doublereal *a, doublereal *b)
{
    /* Local variables */
    doublereal c;
    integer i;

/* DIVIDES P BY THE QUADRATIC  1,U,V  PLACING THE */
/* QUOTIENT IN Q AND THE REMAINDER IN A,B */

    *b = p[0];
    q[0] = *b;
    *a = p[1] - *u * *b;
    q[1] = *a;
    for (i = 2; i < *nn; ++i) {
        c = p[i] - *u * *a - *v * *b;
        q[i] = c;
        *b = *a;
        *a = c;
    }
    return;
} /* quadsd_ */

/* Subroutine */
static void quad_(
  doublereal *a, doublereal *b1, doublereal *c,
  doublereal *sr, doublereal *si, doublereal *lr, doublereal *li
  )
{
    /* Local variables */
    doublereal b, d, e;

/* CALCULATE THE ZEROS OF THE QUADRATIC A*Z**2+B1*Z+C. */
/* THE QUADRATIC FORMULA, MODIFIED TO AVOID */
/* OVERFLOW, IS USED TO FIND THE LARGER ZERO IF THE */
/* ZEROS ARE REAL AND BOTH ZEROS ARE COMPLEX. */
/* THE SMALLER REAL ZERO IS FOUND DIRECTLY FROM THE */
/* PRODUCT OF THE ZEROS C/A. */
    if (*a != 0.) {
        goto L20;
    }
    *sr = 0.;
    if (*b1 != 0.) {
        *sr = -(*c) / *b1;
    }
    *lr = 0.;
L10:
    *si = 0.;
    *li = 0.;
    return;
L20:
    if (*c == 0.) {
        *sr = 0.;
        *lr = -(*b1) / *a;
        goto L10;
    }
/* COMPUTE DISCRIMINANT AVOIDING OVERFLOW */
    b = *b1 / 2.;
    if (abs(b) >= abs(*c)) {
        e = 1. - *a / b * (*c / b);
        d = sqrt(abs(e)) * abs(b);
    }
    else {
        e = *a;
        if (*c < 0.) {
            e = -(*a);
        }
        e = b * (b / abs(*c)) - e;
        d = sqrt(abs(e)) * sqrt(abs(*c));
    }
    if (e < 0.) {
        goto L60;
    }
/* REAL ZEROS */
    if (b >= 0.) {
        d = -d;
    }
    *lr = (-b + d) / *a;
    *sr = 0.;
    if (*lr != 0.) {
        *sr = *c / *lr / *a;
    }
    goto L10;
/* COMPLEX CONJUGATE ZEROS */
L60:
    *sr = -b / *a;
    *lr = *sr;
    *si = abs(d / *a);
    *li = -(*si);
} /* quad_ */