Go to the documentation of this file.
11 int roots(
float* a,
int n,
float* wr,
float* wi) {
12 float sq, b2, c, disc;
30 wr[
m - 2] =
fabs(b2) + sq;
31 if (b2 < 0.0) wr[
m - 2] = -wr[
m - 2];
35 wr[
m - 1] = c / wr[
m - 2];
54 void deflate(
float* a,
int n,
float* b,
float* quad,
float* err) {
63 for (
i = 2;
i <= n;
i++) {
64 b[
i] = a[
i] - r * b[
i - 1] - s * b[
i - 2];
79 void find_quad(
float* a,
int n,
float* b,
float* quad,
float* err,
int* iter) {
80 float *c, dn, dr, ds, drn, dsn, eps, r, s;
92 if (((*iter) % 200) == 0) { eps *= 10.0; }
96 for (
i = 2;
i <= n;
i++) {
97 b[
i] = a[
i] - r * b[
i - 1] - s * b[
i - 2];
98 c[
i] = b[
i] - r * c[
i - 1] - s * c[
i - 2];
100 dn = c[n - 1] * c[n - 3] - c[n - 2] * c[n - 2];
101 drn = b[n] * c[n - 3] - b[n - 1] * c[n - 2];
102 dsn = b[n - 1] * c[n - 1] - b[n] * c[n - 2];
104 if (
fabs(dn) < 1e-10) {
115 }
while ((
fabs(dr) +
fabs(ds)) > eps);
130 for (
i = 1;
i < n;
i++) {
131 b[
i] = a[
i] * ((float) (n -
i)) / coef;
161 float *c, *
x, rs[2], tst;
163 if (
fabs(b[
m]) < 1e-16)
m--;
171 c =
new float[
m + 1];
172 x =
new float[n + 1];
178 tst =
fabs(rs[0] - quad[0]) +
fabs(rs[1] - quad[1]);
184 if (((*iter > 5) && (tst < 1e-4)) || ((*iter > 20) && (tst < 1e-1))) {
186 recurse(a, n, c,
m - 1, rs, err, iter);
199 float *b, *z, err, tmp;
202 if ((tmp = a[0]) != 1.0) {
204 for (
i = 1;
i <= n;
i++) {
217 b =
new float[n + 1];
218 z =
new float[n + 1];
220 for (
i = 0;
i <= n;
i++) {
226 quad[0] = 3.14159e-1;
227 quad[1] = 2.78127e-1;
230 for (
i = 0;
i < 5;
i++) {
232 if ((err > 1e-7) || (iter >
maxiter)) {
234 recurse(z,
m, b,
m - 1, quad, &err, &iter);
237 if (err < 0.001)
break;
242 printf(
"Error! Convergence failure in quadratic x^2 + r*x + s.\n");
245 }
while (err > 0.01);
249 for (
i = 0;
i <=
m;
i++) {
263 float* quad =
new float[2];
264 float*
x =
new float[n];
267 quad[0] = 2.71828e-1;
268 quad[1] = 3.14159e-1;
272 numr =
roots(
x, n, wr, wi);
void diff_poly(float *a, int n, float *b)
friend F32vec4 sqrt(const F32vec4 &a)
void find_quad(float *a, int n, float *b, float *quad, float *err, int *iter)
void polynomComplexRoots(float *wr, float *wi, int n, float *a, int &numr)
void deflate(float *a, int n, float *b, float *quad, float *err)
friend F32vec4 fabs(const F32vec4 &a)
void recurse(float *a, int n, float *b, int m, float *quad, float *err, int *iter)
int roots(float *a, int n, float *wr, float *wi)
void get_quads(float *a, int n, float *quad, float *x)