/* bf_dka -- A solver for the reverse Bessel polynomials using the Durand-Kerner-Aberth method.
* Rev. 1.21 (extended precision, minimal edition) (May 31, 2024) (c) 2024, Takayuki HOSODA
* The calculation results are valid for 14-digits upto 12-th order, 12-digits upto 16-th order
* and 11-digis at 17-th order. C99 (e.g. gcc 4.2 or later) is required for compilation.
* http://www.finetune.co.jp/~lyuka/technote/herb/bessel.html
* DISCLAIMER : This program is free software. It is provided 'as is', without any warranty,
* express or implied, under the extent permitted by applicable law. You are free
* to redistribute and/or modify it.
* This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
*/
#include <float.h>
#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef __FreeBSD__
#include <floatingpoint.h>
#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923
#endif
#ifndef LDBL_EPSILON
#define LDBL_EPSILON 1.0842021724855044340E-19L
#endif
#define DECELERATE ((LDBL_EPSILON) * (long double)(1LL << 32))
#define THRESHOLD ((LDBL_EPSILON) * (long double)(1LL << 27))
#define MAXN 17
#define MAX_ITTR 100
#define F_NORMALIZE 1
#define E_HELP 1
#define E_OPTION 2
#define E_DOMAIN 3
typedef struct fltr {
int order;
} filter;
/* Though omitted for simplicity, the function prototypes should better be written here. */
int flags = 0;
int order;
char * myname = "bf_dka_mini";
char * descr = "A root-finding program for Bessel polynomials using the DKA method.";
char * version = "Rev. 1.21 (extended precision, minimal edition) (May 31, 2024)";
char * copyright = "(c) 2024, Takayuki HOSODA";
char * disclaim = "is free software with ABSOLUTELY NO WARRANTY.";
long double coef[MAXN + 1]; /* global work variables to be used by ostl() */
/* complex polynomial */
long double complex poly(long double complex z, long double *c, int n) {
long double complex p = 0.0;
for (int i = 0; i <= n; i++) p = z * p + c[i];
return p;
}
/* normalize coefficients of the polynomial */
int normalizePoly(long double *a, long double h, int n) {
long double hi = h;
for (int i = n; i >= 0; i--, hi *= h) a[i] *= hi;
return n;
}
/* the function to be solved by ostl() */
long double magfc(long double x) {
extern long double coef[];
extern int order;
return cabsl(coef[order] / poly(x * I, coef, order)) - sqrtl(0.5L);
}
/* ostl - Solve the equation f(x) = 0 to within desired error value for x.
* Based on ost Rev.2.3 (Mar. 3, 2010) (c) 2010 Takayuki HOSODA
* http://www.finetune.co.jp/~lyuka/technote/fract/ost.html
*/
int ostl(long double x1, long double x2, long double (*func)(long double), long double eps, int ittr, long double *x) {
long double a, b, c, d, e, f, t;
int i;
if (func == NULL) return 10;
if (x == NULL) return 9;
if (0 >= ittr) { return 8; }
a = x1; b = x2;
b = (a == b) ? a == 0 ? 0.01 : a * 1.01 : b;
c = (a + b) * 0.5;
if (0.0 == (d = (*func)(a))) { *x = a; return 0; }
if (0.0 == (e = (*func)(b))) { *x = b; return 0; }
if (0.0 == (f = (*func)(c))) { *x = c; return 0; }
if (0.0 == (f - e) * (f - d)) {
c = (a + 3 * b) * 0.25;
if (0.0 == (f = (*func)(c))) { *x = c; return 0; }
if (0.0 == (f - e) * (f - d)) { *x = c; return 7; }
}
for (i = 0; i < ittr; i++) {
if (0.0 == f) { *x = c; return 0; }
if (fabsl(f) <= eps) { *x = c; return 1; }
if (0.0 == (c - a)) { *x = c; return 2; }
if (0 == (f - e)) { *x = c; return 3; }
if ((fabsl(f) > fabsl(e)) && (fabsl(e) <= sqrtl(eps))) { *x = b; return 4; }
/* Ostrowski's method :
t &=& \frac{x_\mathrm{n+2} - x_\mathrm{n+1}}{x_\mathrm{n+2} - x_\mathrm{n}}
\frac{f(x_\mathrm{n+2}) - f(x_\mathrm{n})}{f(x_\mathrm{n+2}) - f(x_\mathrm{n+1})}
\frac{f(x_\mathrm{n+1})}{f(x_\mathrm{n})}\\
x_\mathrm{n+3} &=& \frac{x_\mathrm{n+1} - tx_\mathrm{n}}{1 - t} */
t = e * (f - d) * (c - b);
t /= d * (f - e) * (c - a);
if (t == 1.0) { *x = b; return 6; }
t = (t * a - b) / (t - 1.0);
a = b; b = c; c = t;
d = e; e = f; f = (*func)(c);
}
*x = c; return 5;
}
/* EhrilichAberth method :
\widehat{x_i} &=& x_i - \frac{{p_n(x_i)}}{{{p'_n(x_i)} - {p_n(x_i)} \displaystyle\sum_{\substack{j = 0\\j \neq i}}^ {n-1} \frac{1}{x_i - x_j}}} */
long double complex diffEhrilichAberth(long double complex *x, long double *a, int k, int n) {
long double complex s = 0.0;
long double complex p = 0.0;
long double complex d = 0.0;
long double complex z = x[k];
for (int i = 0; i < n; i++) {
p = z * p + a[i];
d = z * d + a[i] * (n - i); /* d = p'_n(z_i) / p_n(z_i) */
}
p = a[n] + z * p;
for (int i = 0; i < n; i++)
if (i != k) s += 1.0 / (z - x[i]); /* s = \sum_{j = 0, j \neq i}^{n - 1}(1 / (z_i - z_j) */
z = (p == 0.0) ? 0.0 : p / (d - p * s); /* workarround for complex division of zero by infinity */
return z;
}
/* maximum */
long double maxfabsl(long double *a, int n) {
long double max = 0;
for (int i = 0; i <= n; i++) if (max < fabsl(a[i])) max = fabsl(a[i]);
return max;
}
/* magnitude maximum */
long double maxcabsl(long double complex *c, int n) {
long double max = 0;
for (int i = 0; i <= n; i++) if (max < cabsl(c[i])) max = cabsl(c[i]);
return max;
}
/* set Aberth's initial values for c[] to x[] and returns its radious r0 */
/* valid only for the reverse Bessel polynomials (upto at least n = 17) */
long double setAberth(long double complex *x, long double *a, int n) {
long double complex g;
long double r0;
g = - a[1] / n; /* g = -\frac{a_1}{n} */
r0 = fabsl(g) + pow(maxfabsl(a, n), 1.0 / n); /* r_0 = |g| + \sqrt[\displaystyle n]{\max_{0\le k \le n}|a_k|} */
for (int k = 0; k < n; k++)
x[k] = g + r0 * cexp(I * (2 * M_PI * k + M_PI_2 ) / n);
/* x_k = g + r_0\exp\left\{\left(\frac{2\pi k}{n}+\frac{\pi}{2n}}\right)\mathrm{i}\right\}\;|\;0\le k<n,\;\mathrm{i}^2=-1 */
return r0;
}
/* integer reduce */
long long reduce(long long *a, long long *b) {
long long t, p, q;
p = *a;
q = *b;
do {
if (p <= q) {
t = p; p = q; q = t;
}
p %= q;
} while (p);
*a /= q;
*b /= q;
return q;
}
/* factorial */
long long fact(long long ld) {
long long m;
if (ld <= 0) return 1LL;
m = ld--;
while(ld > 1LL)
m *= ld--;
return m;
}
/* partial factorial */
long long pfact(long long j, long long k) {
long long m;
if (j <= 0) return 1LL;
m = j--;
while(j > k)
m *= j--;
return m;
}
/* Bessel polynomial
\displaystyle{a_{i}={\frac {(2n-i)!}{2^{n-i}i!(n-i)!}}\quad i=0..n} */
long long rBessCoef(long double *a, int n) {
long long coef = 1LL;
long long p, q, r, b;
for (int i = n; i >= 0; i--) {
if ((2 * n - i) >= 22) {
q = pfact((2 * n - i), 21); p = fact(i);
r = pfact(21, (n - i)); b = 1LL << (n - i);
reduce(&q, &p); reduce(&r, &b);
reduce(&q, &b); reduce(&r, &p);
coef = (q * r) / (p * b);
a[n - i] = (long double)coef;
} else {
a[n - i] = (long double)(coef = (pfact((2LL * n - i), (n - i)) / (fact(i) * (1LL << (n - i)))));
}
}
return coef;
}
/* the comparator function of qsort() */
int cCompar(const void * c1, const void * c2) {
if (cimagl(*(long double complex *)c1) > cimagl(*(long double complex *)c2)) {
return -1;
} else if (cimagl(*(long double complex *)c1) < cimagl(*(long double complex *)c2)) {
return 1;
} else {
return 0;
}
}
/* polynomial solver by Durand-Aberth-Kerner (Ehrilich-Aberth) method */
long double solvePoly(long double complex *x, long double *a, int n) {
long double ep;
long double er = DBL_MAX;
long double decelerate;
long double complex xn[n];
long double complex dx[n];
for (int k = 1; k < MAX_ITTR; k++) {
decelerate = er <= DECELERATE ? 0.25 : 1.0; /* Decelerate near convergence */
ep = er;
for (int i = 0; i < n; i++)
xn[i] = x[i] - (dx[i] = decelerate * diffEhrilichAberth(x, a, i, n)); /* Ehrilich-Aberth method */
er = maxcabsl(dx, n - 1);
for (int i = 0; i < n; i++) x[i] = xn[i]; /* copy array */
if (er < (DBL_EPSILON * (long double)(2 * n)) || ((ep <= er) && (er < THRESHOLD))) break; /* convergence check */
}
return er;
}
/* print the roots of the reverse Bessel polynomials solved, i.e. poles of the Bessel filters */
void printResults(long double complex *x, int n) {
qsort(x, n, sizeof(long double complex), cCompar); /* sort x[] in descending order of imaginary part */
for (int i = 0; i < n; i++) { /* print results with Q and W */
if (cimagl(x[i]) > DBL_EPSILON) {
printf("# %sPole%d[%d] = %15.11Lf +-%15.11Lf * I (Q = %15.11Lf, w = %15.11Lf)\n",
flags & F_NORMALIZE ? "n" : "", n, i, creall(x[i]), cimagl(x[i]),\
cabsl(x[i]) / (-2 * creall(x[i])), cabsl(x[i]));
} else if (cabsl(cimagl(x[i])) <= DBL_EPSILON) {
printf("# %sPole%d[%d] = %15.11Lf (Q = %15.11Lf, w = %15.11Lf)\n",
flags & F_NORMALIZE ? "n" : "", n, i, creall(x[i]),\
cabsl(x[i]) / (-2 * creall(x[i])), cabsl(x[i]));
}
}
}
/* command-line parser */
int parse(int argc, char *argv[], filter *bf) {
int rtnv = 0;
bf->order = 0;
for (int i = 1; i < argc; i++) {
if (*argv[i] == '-') {
switch (*++argv[i]) {
case 'n':
if (!strncmp(argv[i], "normalize", strlen("normalize"))) {
flags |= F_NORMALIZE;
} else {
rtnv = E_OPTION;
}
break;
case 'o':
if (!strncmp(argv[i], "order=", strlen("order="))) {
sscanf(&argv[i][strlen("order=")], "%d", &bf->order);
if (bf->order <= 0) {
rtnv = E_DOMAIN;
} else if (bf->order > MAXN) {
rtnv = E_DOMAIN;
}
} else {
rtnv = E_OPTION;
}
break;
default:
rtnv = E_OPTION;
break;
}
}
}
if (bf->order == 0) {
bf->order = MAXN;
}
return (rtnv);
}
int main(int argc, char *argv[]) {
extern int order;
extern long double coef[MAXN + 1];
long double a[MAXN + 1]; /* coefficients */
long double complex x[MAXN]; /* approximation values of roots */
long double fc; /* cut-off frequency */
filter bf;
int r;
int n;
#ifdef __FreeBSD__
fpsetprec(FP_PE);
fpsetround(FP_RN);
#endif
if (argc == 1) {
printf("%s -- %s\n %s %s\n%s %s\n Usage : %s [-normalize] -order=n | n = 0..17\n",
myname, descr, version, copyright, myname, disclaim, myname);
return E_HELP;
}
if ((argc == 1) || (r = parse(argc, argv, &bf))) {
printf("%s -- %s\n %s %s\n Usage : %s [-normalize] -order=n | n = 0..17\n", myname, descr, version, copyright, myname);
return r;
} else {
order = n = bf.order; /* set order */
rBessCoef(a, n); /* set reverse Bessel coefficients to a[0..n] */
for (int i = 0; i <= n; i++) coef[i] = a[i]; /* copy array to global work variables to be used by ostl() */
r = ostl(1.0L, 5.0L, magfc, LDBL_EPSILON, MAX_ITTR, &fc); /* solve for cut-off frequency */
if (flags & F_NORMALIZE) normalizePoly(a, fc, n); /* normalize the polynomial coefficients */
printf("# fc = %.11Lg\n", flags & F_NORMALIZE ? 1 : fc); /* print cut-off frequency */
for (int i = n; i >= 0; i--) a[i] /= a[0]; /* make a[] monic */
setAberth(x, a, n); /* set initial values */
solvePoly(x, a, n); /* solve polynomial by DKA method */
printResults(x, n); /* print results */
}
return 0;
}
English edition is here.
![[Mail]](/~lyuka/images/mail.gif)
© 2000 Takayuki HOSODA.
![g &=& -\frac{a_1}{n}\\
r_0 &=& |g| + \sqrt[\displaystyle n]{\max_{0\le k \le n}|a_k|} \\
x_k &=& g + r_0\exp\left\{\mathrm{i}\left(\frac{2\pi k}{n}+\frac{\pi}{2n}}\right)\right\}\quad |\;0\le k<n,\;\mathrm{i}^2=-1\\
{x_k} &\leftarrow& x_k - \frac{{p_n(x_k)}}{{{p'_n(x_k)} - {p_n(x_k)} \displaystyle\sum_{\substack{j = 0\\j \neq k}}^ {n-1} \frac{1}{x_k - x_j}}}](eqn-Durand-Kerner-Ehrlich-Aberth-Bessel-midnightblue.png)
