我有一个具有复数系数的第五多项式。我有一个用C编写的工作程序，我想在其中包含一个寻根算法。例如等式：
x ^ 5 + i x ^ 4 +（2 + i）x ^ 3 +（5 + i）x ^ 2 +（2 + 4i）x + 6 i = 0。
我一直在使用GSL库。例如：
#include <stdio.h>
#include <gsl/gsl_poly.h>
int
main (void)
{
int i;
/* coefficients of P(x) = -1 + x^5 */
double a[6] = { -1, 0, 0, 0, 0, 1 };
double z[10];
gsl_poly_complex_workspace * w
= gsl_poly_complex_workspace_alloc (6);
gsl_poly_complex_solve (a, 6, w, z);
gsl_poly_complex_workspace_free (w);
for (i = 0; i < 5; i++)
{
printf ("z%d = %+.18f %+.18f\n",
i, z[2*i], z[2*i+1]);
}
return 0;
}
但是，它只能用于实系数。我如何才能通过它进行复杂系数的计算。甚至可以用C来做到这一点。
谢谢。
[这是我所做的。 https://wp.csiro.au/alanmiller/toms.html
FORTRAN f77和f90中有一个代码https://mathblog.com/polynomial-root-finding-with-the-jenkins-traub-algorithm/
能够使用f2c库将f77代码转换为c。 cpoly.c的代码如下
/* cpoly.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Common Block Declarations */
struct {
doublereal pr[50], pi[50], hr[50], hi[50], qpr[50], qpi[50], qhr[50], qhi[
50], shr[50], shi[50], sr, si, tr, ti, pvr, pvi, are, mre, eta,
infin;
integer nn;
} global_;
#define global_1 global_
/* Table of constant values */
static integer c__5 = 5;
static integer c__10 = 10;
/* Subroutine */ int cpoly_(doublereal *opr, doublereal *opi, integer *degree,
doublereal *zeror, doublereal *zeroi, logical *fail)
{
/* System generated locals */
integer i__1;
doublereal d__1, d__2;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static integer i__;
static doublereal zi, zr, xx, yy, bnd, xxx;
static integer cnt1, cnt2;
static doublereal base;
extern doublereal cmod_(doublereal *, doublereal *);
extern /* Subroutine */ int mcon_(doublereal *, doublereal *, doublereal *
, doublereal *);
static logical conv;
static doublereal cosr, sinr;
static integer idnn2;
extern doublereal scale_(integer *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *);
extern /* Subroutine */ int cdivid_(doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *);
extern doublereal cauchy_(integer *, doublereal *, doublereal *);
static doublereal smalno;
extern /* Subroutine */ int noshft_(integer *), fxshft_(integer *,
doublereal *, doublereal *, logical *);
/* Added changes from Remark on Algorithm 419 by David H. Withers */
/* CACM (March 1974) Vol 17 No 3 p. 157 */
/* FINDS THE ZEROS OF A COMPLEX POLYNOMIAL. */
/* OPR, OPI - DOUBLE PRECISION VECTORS OF REAL AND */
/* IMAGINARY PARTS OF THE 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 CPOLY */
/* HAS FOUND FEWER THAN DEGREE ZEROS. */
/* THE PROGRAM HAS BEEN WRITTEN TO REDUCE THE CHANCE OF OVERFLOW */
/* OCCURRING. IF IT DOES OCCUR, THERE IS STILL A POSSIBILITY THAT */
/* THE ZEROFINDER WILL WORK PROVIDED THE OVERFLOWED QUANTITY IS */
/* REPLACED BY A LARGE NUMBER. */
/* COMMON AREA */
/* TO CHANGE THE SIZE OF POLYNOMIALS WHICH CAN BE SOLVED, REPLACE */
/* THE DIMENSION OF THE ARRAYS IN THE COMMON AREA. */
/* INITIALIZATION OF CONSTANTS */
/* Parameter adjustments */
--zeroi;
--zeror;
--opi;
--opr;
/* Function Body */
mcon_(&global_1.eta, &global_1.infin, &smalno, &base);
global_1.are = global_1.eta;
global_1.mre = sqrt(2.) * 2. * global_1.eta;
xx = .70710678f;
yy = -xx;
cosr = -.069756474f;
sinr = .99756405f;
*fail = FALSE_;
global_1.nn = *degree + 1;
/* ALGORITHM FAILS IF THE LEADING COEFFICIENT IS ZERO. */
if (opr[1] != 0. || opi[1] != 0.) {
goto L10;
}
*fail = TRUE_;
return 0;
/* REMOVE THE ZEROS AT THE ORIGIN IF ANY. */
L10:
if (opr[global_1.nn] != 0. || opi[global_1.nn] != 0.) {
goto L20;
}
idnn2 = *degree - global_1.nn + 2;
zeror[idnn2] = 0.;
zeroi[idnn2] = 0.;
--global_1.nn;
goto L10;
/* MAKE A COPY OF THE COEFFICIENTS. */
L20:
i__1 = global_1.nn;
for (i__ = 1; i__ <= i__1; ++i__) {
global_1.pr[i__ - 1] = opr[i__];
global_1.pi[i__ - 1] = opi[i__];
global_1.shr[i__ - 1] = cmod_(&global_1.pr[i__ - 1], &global_1.pi[i__
- 1]);
/* L30: */
}
/* SCALE THE POLYNOMIAL. */
bnd = scale_(&global_1.nn, global_1.shr, &global_1.eta, &global_1.infin, &
smalno, &base);
if (bnd == 1.) {
goto L40;
}
i__1 = global_1.nn;
for (i__ = 1; i__ <= i__1; ++i__) {
global_1.pr[i__ - 1] = bnd * global_1.pr[i__ - 1];
global_1.pi[i__ - 1] = bnd * global_1.pi[i__ - 1];
/* L35: */
}
/* START THE ALGORITHM FOR ONE ZERO . */
L40:
if (global_1.nn > 2) {
goto L50;
}
/* CALCULATE THE FINAL ZERO AND RETURN. */
d__1 = -global_1.pr[1];
d__2 = -global_1.pi[1];
cdivid_(&d__1, &d__2, global_1.pr, global_1.pi, &zeror[*degree], &zeroi[*
degree]);
return 0;
/* CALCULATE BND, A LOWER BOUND ON THE MODULUS OF THE ZEROS. */
L50:
i__1 = global_1.nn;
for (i__ = 1; i__ <= i__1; ++i__) {
global_1.shr[i__ - 1] = cmod_(&global_1.pr[i__ - 1], &global_1.pi[i__
- 1]);
/* L60: */
}
bnd = cauchy_(&global_1.nn, global_1.shr, global_1.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_(&c__5);
/* INNER LOOP TO SELECT A SHIFT. */
for (cnt2 = 1; cnt2 <= 9; ++cnt2) {
/* SHIFT IS CHOSEN WITH 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;
/* SECOND STAGE CALCULATION, FIXED SHIFT. */
i__1 = cnt2 * 10;
fxshft_(&i__1, &zr, &zi, &conv);
if (! conv) {
goto L80;
}
/* THE SECOND STAGE JUMPS DIRECTLY TO THE THIRD STAGE ITERATION. */
/* IF SUCCESSFUL THE ZERO IS STORED AND THE POLYNOMIAL DEFLATED. */
idnn2 = *degree - global_1.nn + 2;
zeror[idnn2] = zr;
zeroi[idnn2] = zi;
--global_1.nn;
i__1 = global_1.nn;
for (i__ = 1; i__ <= i__1; ++i__) {
global_1.pr[i__ - 1] = global_1.qpr[i__ - 1];
global_1.pi[i__ - 1] = global_1.qpi[i__ - 1];
/* L70: */
}
goto L40;
L80:
/* IF THE ITERATION IS UNSUCCESSFUL ANOTHER SHIFT IS CHOSEN. */
/* L90: */
;
}
/* IF 9 SHIFTS FAIL, THE OUTER LOOP IS REPEATED WITH ANOTHER */
/* SEQUENCE OF SHIFTS. */
/* L100: */
}
/* THE ZEROFINDER HAS FAILED ON TWO MAJOR PASSES. */
/* RETURN EMPTY HANDED. */
*fail = TRUE_;
return 0;
} /* cpoly_ */
/* Subroutine */ int noshft_(integer *l1)
{
/* System generated locals */
integer i__1, i__2;
doublereal d__1, d__2;
/* Local variables */
static integer i__, j, n;
static doublereal t1, t2;
static integer jj, nm1;
static doublereal xni;
extern doublereal cmod_(doublereal *, doublereal *);
extern /* Subroutine */ int cdivid_(doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *);
/* COMPUTES THE DERIVATIVE POLYNOMIAL AS THE INITIAL H */
/* POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS. */
/* COMMON AREA */
n = global_1.nn - 1;
nm1 = n - 1;
i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
xni = (doublereal) (global_1.nn - i__);
global_1.hr[i__ - 1] = xni * global_1.pr[i__ - 1] / (real) n;
global_1.hi[i__ - 1] = xni * global_1.pi[i__ - 1] / (real) n;
/* L10: */
}
i__1 = *l1;
for (jj = 1; jj <= i__1; ++jj) {
if (cmod_(&global_1.hr[n - 1], &global_1.hi[n - 1]) <= global_1.eta *
10. * cmod_(&global_1.pr[n - 1], &global_1.pi[n - 1])) {
goto L30;
}
d__1 = -global_1.pr[global_1.nn - 1];
d__2 = -global_1.pi[global_1.nn - 1];
cdivid_(&d__1, &d__2, &global_1.hr[n - 1], &global_1.hi[n - 1], &
global_1.tr, &global_1.ti);
i__2 = nm1;
for (i__ = 1; i__ <= i__2; ++i__) {
j = global_1.nn - i__;
t1 = global_1.hr[j - 2];
t2 = global_1.hi[j - 2];
global_1.hr[j - 1] = global_1.tr * t1 - global_1.ti * t2 +
global_1.pr[j - 1];
global_1.hi[j - 1] = global_1.tr * t2 + global_1.ti * t1 +
global_1.pi[j - 1];
/* L20: */
}
global_1.hr[0] = global_1.pr[0];
global_1.hi[0] = global_1.pi[0];
goto L50;
/* IF THE CONSTANT TERM IS ESSENTIALLY ZERO, SHIFT H COEFFICIENTS. */
L30:
i__2 = nm1;
for (i__ = 1; i__ <= i__2; ++i__) {
j = global_1.nn - i__;
global_1.hr[j - 1] = global_1.hr[j - 2];
global_1.hi[j - 1] = global_1.hi[j - 2];
/* L40: */
}
global_1.hr[0] = 0.;
global_1.hi[0] = 0.;
L50:
;
}
return 0;
} /* noshft_ */
/* Subroutine */ int fxshft_(integer *l2, doublereal *zr, doublereal *zi,
logical *conv)
{
/* System generated locals */
integer i__1, i__2;
doublereal d__1, d__2;
/* Local variables */
static integer i__, j, n;
static doublereal oti, otr;
extern doublereal cmod_(doublereal *, doublereal *);
static logical pasd, bool, test;
static doublereal svsi, svsr;
extern /* Subroutine */ int calct_(logical *), nexth_(logical *), vrshft_(
integer *, doublereal *, doublereal *, logical *), polyev_(
integer *, doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *);
/* 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 */
/* COMMON AREA */
n = global_1.nn - 1;
/* EVALUATE P AT S. */
polyev_(&global_1.nn, &global_1.sr, &global_1.si, global_1.pr,
global_1.pi, global_1.qpr, global_1.qpi, &global_1.pvr, &
global_1.pvi);
test = TRUE_;
pasd = FALSE_;
/* CALCULATE FIRST T = -P(S)/H(S). */
calct_(&bool);
/* MAIN LOOP FOR ONE SECOND STAGE STEP. */
i__1 = *l2;
for (j = 1; j <= i__1; ++j) {
otr = global_1.tr;
oti = global_1.ti;
/* COMPUTE NEXT H POLYNOMIAL AND NEW T. */
nexth_(&bool);
calct_(&bool);
*zr = global_1.sr + global_1.tr;
*zi = global_1.si + global_1.ti;
/* TEST FOR CONVERGENCE UNLESS STAGE 3 HAS FAILED ONCE OR THIS */
/* IS THE LAST H POLYNOMIAL . */
if (bool || ! test || j == *l2) {
goto L50;
}
d__1 = global_1.tr - otr;
d__2 = global_1.ti - oti;
if (cmod_(&d__1, &d__2) >= cmod_(zr, zi) * .5) {
goto L40;
}
if (! pasd) {
goto L30;
}
/* THE WEAK CONVERGENCE TEST HAS BEEN PASSED TWICE, START THE */
/* THIRD STAGE ITERATION, AFTER SAVING THE CURRENT H POLYNOMIAL */
/* AND SHIFT. */
i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
global_1.shr[i__ - 1] = global_1.hr[i__ - 1];
global_1.shi[i__ - 1] = global_1.hi[i__ - 1];
/* L10: */
}
svsr = global_1.sr;
svsi = global_1.si;
vrshft_(&c__10, zr, zi, conv);
if (*conv) {
return 0;
}
/* THE ITERATION FAILED TO CONVERGE. TURN OFF TESTING AND RESTORE */
/* H,S,PV AND T. */
test = FALSE_;
i__2 = n;
for (i__ = 1; i__ <= i__2; ++i__) {
global_1.hr[i__ - 1] = global_1.shr[i__ - 1];
global_1.hi[i__ - 1] = global_1.shi[i__ - 1];
/* L20: */
}
global_1.sr = svsr;
global_1.si = svsi;
polyev_(&global_1.nn, &global_1.sr, &global_1.si, global_1.pr,
global_1.pi, global_1.qpr, global_1.qpi, &global_1.pvr, &
global_1.pvi);
calct_(&bool);
goto L50;
L30:
pasd = TRUE_;
goto L50;
L40:
pasd = FALSE_;
L50:
;
}
/* ATTEMPT AN ITERATION WITH FINAL H POLYNOMIAL FROM SECOND STAGE. */
vrshft_(&c__10, zr, zi, conv);
return 0;
} /* fxshft_ */
/* Subroutine */ int vrshft_(integer *l3, doublereal *zr, doublereal *zi,
logical *conv)
{
/* System generated locals */
integer i__1;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static logical b;
static integer i__, j;
static doublereal r1, r2, mp, ms, tp, omp;
extern doublereal cmod_(doublereal *, doublereal *);
static logical bool;
extern /* Subroutine */ int calct_(logical *);
extern doublereal errev_(integer *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *);
extern /* Subroutine */ int nexth_(logical *);
static doublereal relstp;
extern /* Subroutine */ int polyev_(integer *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *);
/* 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 */
/* COMMON AREA */
*conv = FALSE_;
b = FALSE_;
global_1.sr = *zr;
global_1.si = *zi;
/* MAIN LOOP FOR STAGE THREE */
i__1 = *l3;
for (i__ = 1; i__ <= i__1; ++i__) {
/* EVALUATE P AT S AND TEST FOR CONVERGENCE. */
polyev_(&global_1.nn, &global_1.sr, &global_1.si, global_1.pr,
global_1.pi, global_1.qpr, global_1.qpi, &global_1.pvr, &
global_1.pvi);
mp = cmod_(&global_1.pvr, &global_1.pvi);
ms = cmod_(&global_1.sr, &global_1.si);
if (mp > errev_(&global_1.nn, global_1.qpr, global_1.qpi, &ms, &mp, &
global_1.are, &global_1.mre) * 20.) {
goto L10;
}
/* POLYNOMIAL VALUE IS SMALLER IN VALUE THAN A BOUND ON THE ERROR */
/* IN EVALUATING P, TERMINATE THE ITERATION. */
*conv = TRUE_;
*zr = global_1.sr;
*zi = global_1.si;
return 0;
L10:
if (i__ == 1) {
goto L40;
}
if (b || mp < omp || relstp >= .05) {
goto L30;
}
/* 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 = TRUE_;
if (relstp < global_1.eta) {
tp = global_1.eta;
}
r1 = sqrt(tp);
r2 = global_1.sr * (r1 + 1.) - global_1.si * r1;
global_1.si = global_1.sr * r1 + global_1.si * (r1 + 1.);
global_1.sr = r2;
polyev_(&global_1.nn, &global_1.sr, &global_1.si, global_1.pr,
global_1.pi, global_1.qpr, global_1.qpi, &global_1.pvr, &
global_1.pvi);
for (j = 1; j <= 5; ++j) {
calct_(&bool);
nexth_(&bool);
/* L20: */
}
omp = global_1.infin;
goto L50;
/* EXIT IF POLYNOMIAL VALUE INCREASES SIGNIFICANTLY. */
L30:
if (mp * .1 > omp) {
return 0;
}
L40:
omp = mp;
/* CALCULATE NEXT ITERATE. */
L50:
calct_(&bool);
nexth_(&bool);
calct_(&bool);
if (bool) {
goto L60;
}
relstp = cmod_(&global_1.tr, &global_1.ti) / cmod_(&global_1.sr, &
global_1.si);
global_1.sr += global_1.tr;
global_1.si += global_1.ti;
L60:
;
}
return 0;
} /* vrshft_ */
/* Subroutine */ int calct_(logical *bool)
{
/* System generated locals */
doublereal d__1, d__2;
/* Local variables */
static integer n;
static doublereal hvi, hvr;
extern doublereal cmod_(doublereal *, doublereal *);
extern /* Subroutine */ int cdivid_(doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *), polyev_(
integer *, doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *);
/* COMPUTES T = -P(S)/H(S). */
/* BOOL - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO. */
/* COMMON AREA */
n = global_1.nn - 1;
/* EVALUATE H(S). */
polyev_(&n, &global_1.sr, &global_1.si, global_1.hr, global_1.hi,
global_1.qhr, global_1.qhi, &hvr, &hvi);
*bool = cmod_(&hvr, &hvi) <= global_1.are * 10. * cmod_(&global_1.hr[n -
1], &global_1.hi[n - 1]);
if (*bool) {
goto L10;
}
d__1 = -global_1.pvr;
d__2 = -global_1.pvi;
cdivid_(&d__1, &d__2, &hvr, &hvi, &global_1.tr, &global_1.ti);
return 0;
L10:
global_1.tr = 0.;
global_1.ti = 0.;
return 0;
} /* calct_ */
/* Subroutine */ int nexth_(logical *bool)
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer j, n;
static doublereal t1, t2;
static integer nm1;
/* CALCULATES THE NEXT SHIFTED H POLYNOMIAL. */
/* BOOL - LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO */
/* COMMON AREA */
n = global_1.nn - 1;
nm1 = n - 1;
if (*bool) {
goto L20;
}
i__1 = n;
for (j = 2; j <= i__1; ++j) {
t1 = global_1.qhr[j - 2];
t2 = global_1.qhi[j - 2];
global_1.hr[j - 1] = global_1.tr * t1 - global_1.ti * t2 +
global_1.qpr[j - 1];
global_1.hi[j - 1] = global_1.tr * t2 + global_1.ti * t1 +
global_1.qpi[j - 1];
/* L10: */
}
global_1.hr[0] = global_1.qpr[0];
global_1.hi[0] = global_1.qpi[0];
return 0;
/* IF H(S) IS ZERO REPLACE H WITH QH. */
L20:
i__1 = n;
for (j = 2; j <= i__1; ++j) {
global_1.hr[j - 1] = global_1.qhr[j - 2];
global_1.hi[j - 1] = global_1.qhi[j - 2];
/* L30: */
}
global_1.hr[0] = 0.;
global_1.hi[0] = 0.;
return 0;
} /* nexth_ */
/* Subroutine */ int polyev_(integer *nn, doublereal *sr, doublereal *si,
doublereal *pr, doublereal *pi, doublereal *qr, doublereal *qi,
doublereal *pvr, doublereal *pvi)
{
/* System generated locals */
integer i__1;
/* Local variables */
static integer i__;
static doublereal t;
/* EVALUATES A POLYNOMIAL P AT S BY THE HORNER RECURRENCE */
/* PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV. */
/* Parameter adjustments */
--qi;
--qr;
--pi;
--pr;
/* Function Body */
qr[1] = pr[1];
qi[1] = pi[1];
*pvr = qr[1];
*pvi = qi[1];
i__1 = *nn;
for (i__ = 2; i__ <= i__1; ++i__) {
t = *pvr * *sr - *pvi * *si + pr[i__];
*pvi = *pvr * *si + *pvi * *sr + pi[i__];
*pvr = t;
qr[i__] = *pvr;
qi[i__] = *pvi;
/* L10: */
}
return 0;
} /* polyev_ */
doublereal errev_(integer *nn, doublereal *qr, doublereal *qi, doublereal *ms,
doublereal *mp, doublereal *are, doublereal *mre)
{
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Local variables */
static doublereal e;
static integer i__;
extern doublereal cmod_(doublereal *, doublereal *);
/* 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 */
/* Parameter adjustments */
--qi;
--qr;
/* Function Body */
e = cmod_(&qr[1], &qi[1]) * *mre / (*are + *mre);
i__1 = *nn;
for (i__ = 1; i__ <= i__1; ++i__) {
e = e * *ms + cmod_(&qr[i__], &qi[i__]);
/* L10: */
}
ret_val = e * (*are + *mre) - *mp * *mre;
return ret_val;
} /* errev_ */
doublereal cauchy_(integer *nn, doublereal *pt, doublereal *q)
{
/* System generated locals */
integer i__1;
doublereal ret_val, d__1;
/* Builtin functions */
double log(doublereal), exp(doublereal);
/* Local variables */
static doublereal f;
static integer i__, n;
static doublereal x, df, dx, xm;
/* CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A */
/* POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS. */
/* Parameter adjustments */
--q;
--pt;
/* Function Body */
pt[*nn] = -pt[*nn];
/* COMPUTE UPPER ESTIMATE OF BOUND. */
n = *nn - 1;
x = exp((log(-pt[*nn]) - log(pt[1])) / (real) n);
if (pt[n] == 0.) {
goto L20;
}
/* IF NEWTON STEP AT THE ORIGIN IS BETTER, USE IT. */
xm = -pt[*nn] / pt[n];
if (xm < x) {
x = xm;
}
/* CHOP THE INTERVAL (0,X) UNITL F LE 0. */
L20:
xm = x * .1;
f = pt[1];
i__1 = *nn;
for (i__ = 2; i__ <= i__1; ++i__) {
f = f * xm + pt[i__];
/* L30: */
}
if (f <= 0.) {
goto L40;
}
x = xm;
goto L20;
L40:
dx = x;
/* DO NEWTON ITERATION UNTIL X CONVERGES TO TWO DECIMAL PLACES. */
L50:
if ((d__1 = dx / x, abs(d__1)) <= .005) {
goto L70;
}
q[1] = pt[1];
i__1 = *nn;
for (i__ = 2; i__ <= i__1; ++i__) {
q[i__] = q[i__ - 1] * x + pt[i__];
/* L60: */
}
f = q[*nn];
df = q[1];
i__1 = n;
for (i__ = 2; i__ <= i__1; ++i__) {
df = df * x + q[i__];
/* L65: */
}
dx = f / df;
x -= dx;
goto L50;
L70:
ret_val = x;
return ret_val;
} /* cauchy_ */
doublereal scale_(integer *nn, doublereal *pt, doublereal *eta, doublereal *
infin, doublereal *smalno, doublereal *base)
{
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Builtin functions */
double sqrt(doublereal), log(doublereal), pow_di(doublereal *, integer *);
/* Local variables */
static integer i__, l;
static doublereal x, hi, sc, lo, min__, max__;
/* 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. */
/* FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS. */
/* Parameter adjustments */
--pt;
/* Function Body */
hi = sqrt(*infin);
lo = *smalno / *eta;
max__ = 0.;
min__ = *infin;
i__1 = *nn;
for (i__ = 1; i__ <= i__1; ++i__) {
x = pt[i__];
if (x > max__) {
max__ = x;
}
if (x != 0. && x < min__) {
min__ = x;
}
/* L10: */
}
/* SCALE ONLY IF THERE ARE VERY LARGE OR VERY SMALL COMPONENTS. */
ret_val = 1.;
if (min__ >= lo && max__ <= hi) {
return ret_val;
}
x = lo / min__;
if (x > 1.) {
goto L20;
}
sc = 1. / (sqrt(max__) * sqrt(min__));
goto L30;
L20:
sc = x;
if (*infin / sc > max__) {
sc = 1.;
}
L30:
l = (integer) (log(sc) / log(*base) + .5f);
ret_val = pow_di(base, &l);
return ret_val;
} /* scale_ */
/* Subroutine */ int cdivid_(doublereal *ar, doublereal *ai, doublereal *br,
doublereal *bi, doublereal *cr, doublereal *ci)
{
static doublereal d__, r__, t;
extern /* Subroutine */ int mcon_(doublereal *, doublereal *, doublereal *
, doublereal *);
static doublereal infin;
/* COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW. */
if (*br != 0. || *bi != 0.) {
goto L10;
}
/* DIVISION BY ZERO, C = INFINITY. */
mcon_(&t, &infin, &t, &t);
*cr = infin;
*ci = infin;
return 0;
L10:
if (abs(*br) >= abs(*bi)) {
goto L20;
}
r__ = *br / *bi;
d__ = *bi + r__ * *br;
*cr = (*ar * r__ + *ai) / d__;
*ci = (*ai * r__ - *ar) / d__;
return 0;
L20:
r__ = *bi / *br;
d__ = *br + r__ * *bi;
*cr = (*ar + *ai * r__) / d__;
*ci = (*ai - *ar * r__) / d__;
return 0;
} /* cdivid_ */
doublereal cmod_(doublereal *r__, doublereal *i__)
{
/* System generated locals */
doublereal ret_val, d__1;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static doublereal ai, ar;
/* MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW. */
ar = abs(*r__);
ai = abs(*i__);
if (ar >= ai) {
goto L10;
}
/* Computing 2nd power */
d__1 = ar / ai;
ret_val = ai * sqrt(d__1 * d__1 + 1.);
return ret_val;
L10:
if (ar <= ai) {
goto L20;
}
/* Computing 2nd power */
d__1 = ai / ar;
ret_val = ar * sqrt(d__1 * d__1 + 1.);
return ret_val;
L20:
ret_val = ar * sqrt(2.);
return ret_val;
} /* cmod_ */
/* Subroutine */ int mcon_(doublereal *eta, doublereal *infiny, doublereal *
smalno, doublereal *base)
{
/* System generated locals */
integer i__1, i__2;
doublereal d__1;
/* Builtin functions */
double pow_di(doublereal *, integer *);
/* Local variables */
static integer m, n, t;
/* 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.0D0 + ETA IS */
/* GREATER THAN 1.0D0. */
/* INFINY THE LARGEST FLOATING-POINT NUMBER */
/* SMALNO THE SMALLEST POSITIVE FLOATING-POINT NUMBER */
/* BASE THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED */
/* LET T BE THE NUMBER OF BASE-DIGITS IN EACH FLOATING-POINT */
/* NUMBER(DOUBLE PRECISION). THEN ETA IS EITHER .5*B**(1-T) */
/* OR B**(1-T) DEPENDING ON WHETHER ROUNDING OR TRUNCATION */
/* IS USED. */
/* LET M BE THE LARGEST EXPONENT AND N THE SMALLEST EXPONENT */
/* IN THE NUMBER SYSTEM. THEN INFINY IS (1-BASE**(-T))*BASE**M */
/* AND SMALNO IS BASE**N. */
/* THE VALUES FOR BASE,T,M,N BELOW CORRESPOND TO THE IBM/360. */
*base = 16.;
t = 14;
m = 63;
n = -65;
i__1 = 1 - t;
*eta = pow_di(base, &i__1);
i__1 = -t;
i__2 = m - 1;
*infiny = *base * (1. - pow_di(base, &i__1)) * pow_di(base, &i__2);
i__1 = n + 3;
/* Computing 3rd power */
d__1 = *base;
*smalno = pow_di(base, &i__1) / (d__1 * (d__1 * d__1));
return 0;
} /* mcon_ */
/*
int main(){
doublereal opr[11] = {1, -55, 1320, -18150, 157773, -902055, 3416930, -8409500, 12753576, -10628640, 3628800};
doublereal opi[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
integer degree = 10;
doublereal zr[11];
doublereal zi[11];
logical fail;
cpoly_(opr, opi, °ree, zr, zi, &fail);
}
*/
下面是测试代码
#include<stdio.h>
#include "f2c.h"
extern int cpoly_(doublereal *opr, doublereal *opi, integer *degree, doublereal *zeror, doublereal *zeroi, logical *fail);
int main(){
doublereal opr[11] = {1, -55, 1320, -18150, 157773, -902055, 3416930, -8409500, 12753576, -10628640, 3628800};
doublereal opi[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
integer degree = 10;
doublereal zr[11];
doublereal zi[11];
logical fail;
int res = cpoly_(opr, opi, °ree, zr, zi, &fail);
printf("%d %d\n", res, fail);
for(int i=0; i<11; i++){
printf("%e %e\n", zr[i], zi[i]);
}
}
编译和链接可以通过以下方式完成
gcc -c cpoly.c
gcc -c cpoly_test.c
gcc cpoly_test.o cpoly.o -lm -lf2c -o test
但是，需要安装f2c库。这是链接：
https://askubuntu.com/questions/1077151/how-to-run-a-csh-script-to-install-f2c