1 /*
2 * This file is part of MPSolve 3.2.1
3 *
4 * Copyright (C) 2001-2020, Dipartimento di Matematica "L. Tonelli", Pisa.
5 * License: http://www.gnu.org/licenses/gpl.html GPL version 3 or higher
6 *
7 * Authors:
8 * Dario Andrea Bini <bini@dm.unipi.it>
9 * Giuseppe Fiorentino <fiorent@dm.unipi.it>
10 * Leonardo Robol <leonardo.robol@unipi.it>
11 */
12
13
14 #include <math.h>
15 #include <float.h>
16 #include <limits.h>
17 #include <assert.h>
18 #include <mps/mps.h>
19
20 #define MPS_STARTING_SIGMA (0.66 * (PI / s->n))
21 #define pi2 6.283184
22
23 /**
24 * @brief Compute the greatest common divisor of <code>a</code>
25 * and <code>b</code>.
26 *
27 * @param a first integer
28 * @param b second integer
29 * @return <code>GCD(a,b)</code>
30 */
31 static int
mps_gcd(int a,int b)32 mps_gcd (int a, int b)
33 {
34 int temp;
35
36 do
37 {
38 temp = b;
39 b = a % b;
40 a = temp;
41 } while (b != 0);
42 return a;
43 }
44
45 /**
46 * @brief Find the sigma that maximize distance between
47 * starting approximation in the last annulus and the one
48 * in the new annulus.
49 *
50 * This function also set <code>s->last_sigma</code> to the value
51 * that is computed.
52 *
53 * @param s the mps_context struct pointer.
54 * @param last_sigma the last value of sigma.
55 * @param cluster_item The element of the <code>mps_clusterization</code> of which
56 * we are computing the starting points, or NULL if we are computing the starting
57 * points for all the approximations.
58 * @param n the number of roots in the cluster.
59 * @return the shift advised for the starting approximation in the annulus.
60 */
61 static double
mps_maximize_distance(mps_context * s,double last_sigma,mps_cluster_item * cluster_item,int n)62 mps_maximize_distance (mps_context * s, double last_sigma,
63 mps_cluster_item * cluster_item, int n)
64 {
65 double delta_sigma;
66
67 /* Find number of roots in the last cluster */
68 if (!cluster_item || cluster_item->prev == NULL)
69 return s->last_sigma;
70 int old_clust_n = (cluster_item->prev->cluster->n);
71
72 /* Compute right shifting angle for the new approximations, i.e.
73 * pi / [m,n] where [m,n] is the least common multiply of m and n.
74 * This is done by computing the gcd and then dividing m*n by it. */
75 delta_sigma = PI * (old_clust_n * mps_gcd (old_clust_n, n)) / (4 * n);
76
77 /* Return shifted value, archiving it for the next pass */
78 s->last_sigma = last_sigma + delta_sigma;
79 return s->last_sigma;
80 }
81
82 /**
83 * @brief Compute radii of the circles where the initial approximation
84 * will be disposed by mps_fstart()
85 *
86 * @param s mps_context* stuct pointer.
87 * @param n number of roots in the cluster.
88 * @param cluster_item The element of the <code>mps_clusterization</code> of which
89 * we are computing the starting points, or NULL if we are computing the starting
90 * points for all the approximations.
91 * @param clust_rad radius of the cluster.
92 * @param g new gravity center where the polynomial has been shifted.
93 * @param eps out epsilon.
94 * @param fap[] Array with the moduli of the coefficients.
95 *
96 * @return A mps_starting_configuration containing the necessary information
97 * to dispose the initial approximations on the circles.
98 *
99 * @see mps_fstart()
100 */
101 MPS_PRIVATE mps_starting_configuration
mps_fcompute_starting_radii(mps_context * s,int n,mps_cluster_item * cluster_item,double clust_rad,double g,rdpe_t eps,double fap[])102 mps_fcompute_starting_radii (mps_context * s, int n, mps_cluster_item * cluster_item,
103 double clust_rad, double g, rdpe_t eps,
104 double fap[])
105 {
106 MPS_DEBUG_THIS_CALL (s);
107
108 const double big = DBL_MAX, small = DBL_MIN;
109 const double xbig = log (big), xsmall = log (small);
110
111 int i, j, k, nzeros, iold, ni, offset;
112 double temp, r;
113
114 mps_starting_configuration c = {
115 .n_radii = 0,
116 .partitioning = NULL,
117 .fradii = NULL,
118 .dradii = NULL
119 };
120
121 ni = 0;
122 nzeros = 0;
123 r = 0.0;
124
125 /**********************************************
126 check for possible null entries in the trailing
127 coefficients only in the case where the polynomial
128 has been shifted in g, replace null coefficients
129 with small numbers according to the working precision
130 and to the number of null coefficients
131 **********************************************/
132 if (g != 0.0)
133 {
134 for (i = 0; i <= n; i++)
135 if (fap[i] != 0.0)
136 {
137 ni = i;
138 break;
139 }
140 if (ni == 0)
141 temp = 2 * xsmall;
142 else
143 temp = log (fap[ni]) + ni * (log (DBL_EPSILON) + log (g * ni * 10.0));
144 }
145 else
146 temp = 2 * xsmall;
147
148 for (i = 0; i <= n; i++)
149 {
150 if (fap[i] != 0.0)
151 s->fap2[i] = log (fap[i]);
152 else
153 s->fap2[i] = temp;
154 }
155
156 /* Compute convex hull */
157 int * h = mps_fconvex (s, n, s->fap2);
158
159 j = 0;
160 for (i = 1; i <= n; i++)
161 if (h[i])
162 j++;
163
164 c.fradii = mps_newv (double, j + 2);
165 c.partitioning = mps_newv (int, j + 2);
166
167 /* compute the radii of the circles containing starting approximations */
168 c.n_radii = 0;
169 c.partitioning[0] = 0;
170
171 for (i = 1; i <= n; i++)
172 if (h[i])
173 {
174 iold = c.partitioning[c.n_radii];
175 nzeros = i - iold;
176 temp = (s->fap2[iold] - s->fap2[i]) / nzeros;
177
178 /* if the radius is too small to be represented as double, set it
179 * to the minimum representable double */
180 if (temp < xsmall) /* if (temp < MAX(xsmall, -xbig)) DARIO Giugno 23 */
181 r = DBL_MIN; /* r = small; */
182
183 /* if the radius is too big to be represented as double, set it
184 * to the maximum representable double */
185 if (temp > xbig)
186 r = DBL_MAX; /* big; DARIO Giugno 23 */
187
188 /* if the radius is representable as double, compute it */
189 if ((temp <= xbig) && (temp > xsmall))
190 /* if ((temp <= xbig) && (temp > MAX(-xbig, xsmall))) DARIO Giugno 23 */
191 r = exp (temp);
192
193 /* if the radius is greater than the radius of the cluster
194 * set the radius equal to the radius of the cluster */
195 if (clust_rad != 0 && r > clust_rad)
196 r = clust_rad;
197
198 c.fradii[c.n_radii] = r;
199 c.partitioning[++c.n_radii] = i;
200 }
201
202 /* Close partitioning */
203 c.partitioning[c.n_radii] = n;
204
205 /* Compact radius that are too near */
206 for (i = 0; i < c.n_radii; i++)
207 {
208 /* Scan next radii to see if they are near the
209 * i-th that we are considering now */
210 for (j = i + 1; j < c.n_radii; j++)
211 {
212 /* Get an estimate of the distance between approximations that would
213 * be obtained collapsing the circles. It's a little understimated
214 * but this is due to the fact that changing approximation is only
215 * a fallback, and not the preferred action to perform. */
216 if (fabs ((c.fradii[j] - c.fradii[i])) >
217 MIN (c.fradii[j],
218 c.fradii[i]) * PI / (c.partitioning[j + 1] -
219 c.partitioning[i]))
220 {
221 break;
222 }
223 }
224
225 /* This is the number of circles that are near */
226 j--;
227 offset = j - i;
228
229 /* If there is nothing to compact, do not compact it */
230 if (offset == 0)
231 {
232 continue;
233 }
234
235 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
236 MPS_DEBUG (s, "Compacting circles from %d to %d", i, j);
237
238 /* We shall now compact circles between i and j, so
239 * we start computing the mean of the radius */
240 c.fradii[i] *= (c.partitioning[i + 1] - c.partitioning[i]);
241 for (k = i + 1; k <= j; k++)
242 {
243 c.fradii[i] +=
244 c.fradii[k] * (c.partitioning[k + 1] - c.partitioning[k]);
245 }
246
247 c.fradii[i] /= c.partitioning[j + 1] - c.partitioning[i];
248
249 /* Move other circles backward */
250 for (k = j + 1; k < c.n_radii; k++)
251 {
252 c.fradii[k - offset] = c.fradii[k];
253 c.partitioning[k - offset] = c.partitioning[k];
254 }
255
256 /* Set new c.n_radii and new partitioning */
257 c.n_radii = c.n_radii - offset;
258 c.partitioning[c.n_radii] = n;
259 }
260
261 free (h);
262
263 return c;
264 }
265
266 /**
267 * @brief Compute new starting approximations to the roots
268 * of the polynomial \f$p(x)\f$ having coefficients of modulus apoly.
269 *
270 * Computations is done by
271 * means of the Rouche'-based criterion of Bini (Numer. Algo. 1996).
272 * The program can compute all the approximations
273 * (if \f$n\f$ is the degree of \f$p(x)\f$) or it may compute the
274 * approximations of the cluster in the <code>cluster_item</code>.
275 * The status vector is changed into <code>'o'</code> for the components
276 * that belong to a cluster with relative radius less than <code>eps</code>.
277 * The status vector is changed into <code>'x'</code> for the components that
278 * cannot be represented as double.
279 *
280 * @param s The <code>mps_context</code> associated with the current computation.
281 * @param n number of roots in the cluster.
282 * @param cluster_item The element of the <code>mps_clusterization</code> of which
283 * we are computing the starting points, or NULL if we are computing the starting
284 * points for all the approximations.
285 * @param clust_rad radius of cluster.
286 * @param g gravity center of the cluster.
287 * @param eps a double that represent the maximum value
288 * of relative radius (with respect to <code>g</code>) of
289 * roots whose status must be set to <code>o</code>.
290 * @param fap array of moduli of the coefficients as double.
291 *
292 * @see status
293 */
294 MPS_PRIVATE void
mps_fstart(mps_context * s,int n,mps_cluster_item * cluster_item,double clust_rad,double g,rdpe_t eps,double fap[])295 mps_fstart (mps_context * s, int n, mps_cluster_item * cluster_item,
296 double clust_rad, double g, rdpe_t eps, double fap[])
297 {
298 MPS_DEBUG_THIS_CALL (s);
299 int i, j, jj, l, nzeros = 0;
300 double sigma, th, ang, r = 0;
301 mps_root * it_root = NULL;
302 rdpe_t tmp;
303
304 mps_cluster * cluster = NULL;
305 mps_root * root = NULL;
306
307 if (cluster_item)
308 cluster = cluster_item->cluster;
309
310 if (s->random_seed)
311 sigma = drand ();
312 else
313 {
314 /* If this is the first cluster select sigma = 0. In the other
315 * case try to maximize starting points distance. */
316 if ((cluster == NULL) || (cluster_item->prev == NULL))
317 {
318 sigma = s->last_sigma = MPS_STARTING_SIGMA;
319 }
320 else
321 {
322 sigma = mps_maximize_distance (s, s->last_sigma, cluster_item, n);
323 }
324 }
325
326 th = pi2 / n;
327
328 /* In the general case apply the Rouche-based criterion */
329 mps_starting_configuration c =
330 mps_fcompute_starting_radii (s, n, cluster_item, clust_rad, g, eps, fap);
331
332 if (g != 0.0 && cluster != NULL)
333 it_root = cluster->first;
334
335 for (i = 0; i < c.n_radii; i++)
336 {
337 nzeros = c.partitioning[i + 1] - c.partitioning[i];
338 ang = pi2 / nzeros;
339 r = c.fradii[i];
340
341 for (j = c.partitioning[i]; j < c.partitioning[i + 1]; j++)
342 {
343 if (g != 0.0 && it_root)
344 l = it_root->k;
345 else
346 l = j;
347
348 jj = j - c.partitioning[i];
349
350 /* if the radius reaches extreme values then set the component
351 * of status, corresponding to approximation which fall out the
352 * representable range, to 'x' (out) */
353 if ((r == DBL_MIN) || (r == DBL_MAX))
354 /* if ((r == small) || (r == big)) DARIO Giugno 23 */
355 s->root[l]->status = MPS_ROOT_STATUS_NOT_FLOAT;
356 cplx_set_d (s->root[l]->fvalue,
357 r * cos (ang * jj + th * c.partitioning[i + 1] +
358 sigma),
359 r * sin (ang * jj + th * c.partitioning[i + 1] +
360 sigma));
361 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
362 {
363 MPS_DEBUG_CPLX (s, s->root[l]->fvalue, "s->froot[%d]", l);
364 }
365
366 if (it_root)
367 it_root = it_root->next;
368 }
369
370
371 /* If the new radius of the cluster is relatively smaller, then
372 * set the status component equal to 'o' (output) */
373 if (g != 0.0 && cluster)
374 {
375 rdpe_mul_d (tmp, eps, g);
376 if (r * nzeros <= rdpe_get_d (tmp))
377 {
378 for (root = cluster->first; root != NULL; root = root->next)
379 {
380 l = root->k;
381 s->root[l]->status = MPS_ROOT_STATUS_APPROXIMATED_IN_CLUSTER;
382 s->root[l]->frad = r * nzeros;
383 }
384 }
385 }
386 }
387
388 mps_starting_configuration_clear (s, &c);
389 }
390
391 /**
392 * @brief Compute radii of the circles where the initial approximation
393 * will be disposed by mps_dstart()
394 *
395 * @param s mps_context* stuct pointer.
396 * @param n number of roots in the cluster.
397 * @param cluster_item The element of the <code>mps_clusterization</code> of which
398 * we are computing the starting points, or NULL if we are computing the starting
399 * points for all the approximations.
400 * @param clust_rad radius of the cluster.
401 * @param g new gravity center where the polynomial has been shifted.
402 * @param eps out epsilon.
403 * @param dap[] Array with the moduli of the coefficients.
404 *
405 * @see mps_dstart()
406 */
407 MPS_PRIVATE mps_starting_configuration
mps_dcompute_starting_radii(mps_context * s,int n,mps_cluster_item * cluster_item,rdpe_t clust_rad,rdpe_t g,rdpe_t eps,rdpe_t dap[])408 mps_dcompute_starting_radii (mps_context * s, int n,
409 mps_cluster_item * cluster_item,
410 rdpe_t clust_rad, rdpe_t g, rdpe_t eps,
411 rdpe_t dap[])
412 {
413 /* Compute big and small values */
414 double xbig, xsmall;
415
416 int i, j, k, nzeros, ni, offset, iold;
417 double temp;
418 rdpe_t r, tmp;
419
420 xbig = LONG_MAX;
421 xsmall = LONG_MIN;
422
423 mps_starting_configuration c = {
424 .n_radii = 0,
425 .partitioning = NULL,
426 .fradii = NULL,
427 .dradii = NULL
428 };
429
430 ni = 0;
431 nzeros = 0;
432 rdpe_set (r, rdpe_zero);
433
434 if (rdpe_ne (g, rdpe_zero))
435 {
436 for (i = 0; i <= n; i++)
437 if (rdpe_ne (dap[i], rdpe_zero))
438 {
439 ni = i;
440 break;
441 }
442 if (ni == 0)
443 temp = -2.0 * (LONG_MAX * LOG2);
444 else
445 {
446 /* temp = log(dap[ni])+ni*(log(DBL_EPSILON)+log(g*ni*10.d0)) */
447 temp = ni * 10.0;
448 rdpe_mul_d (tmp, g, temp);
449 temp = rdpe_log (tmp);
450 temp += log (DBL_EPSILON);
451 temp *= ni;
452 temp += rdpe_log (dap[ni]);
453 }
454 }
455 else
456 temp = -2.0 * (LONG_MAX * LOG2);
457 for (i = 0; i <= n; i++)
458 if (rdpe_ne (dap[i], rdpe_zero))
459 s->fap2[i] = rdpe_log (dap[i]);
460 else
461 s->fap2[i] = temp;
462
463 /* compute the convex hull */
464 int * h = mps_fconvex (s, n, s->fap2);
465
466 /* Count the number of vertexes in the convex hull. */
467 j = 0;
468 for (i = 1; i <= n; i++)
469 if (h[i])
470 j++;
471
472 c.dradii = mps_newv (rdpe_t, j + 2);
473 c.partitioning = mps_newv (int, j + 2);
474
475 /* compute the radii of the circles containing starting approximations */
476 c.n_radii = 0;
477 c.partitioning[0] = 0;
478 for (i = 1; i <= n; i++)
479 if (h[i])
480 {
481 iold = c.partitioning[c.n_radii];
482 nzeros = i - iold;
483 temp = (s->fap2[iold] - s->fap2[i]) / nzeros;
484
485 /* if the radius is too small to be represented as double, set it
486 * to the minimum representable double */
487 if (temp < xsmall)
488 rdpe_set (r, RDPE_MIN); /* r = small; */
489
490 /* if the radius is too big to be represented as double, set it
491 * to the maximum representable double */
492 if (temp > xbig)
493 rdpe_set (r, RDPE_MAX);
494
495 /* if the radius is representable as double, compute it */
496 if ((temp < xbig) && (temp > xsmall))
497 rdpe_set_d (r, temp);
498 rdpe_exp_eq (r);
499
500 /* if the radius is greater than the radius of the cluster
501 * set the radius equal to the radius of the cluster */
502 if (rdpe_ne (clust_rad, rdpe_zero) && rdpe_gt (r, clust_rad))
503 rdpe_set (r, clust_rad);
504
505 rdpe_set (c.dradii[c.n_radii], r);
506 c.partitioning[++c.n_radii] = i;
507 }
508
509 /* Close partitioning */
510 c.partitioning[c.n_radii] = n;
511
512 /* Compact radius that are too near */
513 for (i = 0; i < c.n_radii; i++)
514 {
515 /* Scan next radii to see if they are near the
516 * i-th that we are considering now */
517 for (j = i + 1; j < c.n_radii; j++)
518 {
519 /* Get an estimate of the distance between approximations that would
520 * be obtained collapsing the circles. It's a little understimated
521 * but this is due to the fact that changing approximation is only
522 * a fallback, and not the preferred action to perform. */
523 rdpe_sub (tmp, c.dradii[j], c.dradii[i]);
524 rdpe_abs_eq (tmp);
525 if (rdpe_lt (c.dradii[i], c.dradii[j]))
526 rdpe_div_eq (tmp, c.dradii[i]);
527 else
528 rdpe_div_eq (tmp, c.dradii[j]);
529 rdpe_div_eq_d (tmp, PI);
530 rdpe_mul_eq_d (tmp, c.partitioning[j + 1] - c.partitioning[i]);
531 if (rdpe_gt (tmp, rdpe_one))
532 {
533 break;
534 }
535 }
536
537 /* This is the number of circles that are near */
538 j--;
539 offset = j - i;
540
541
542 /* If there is nothing to compact, do not compact it */
543 if (offset == 0)
544 {
545 continue;
546 }
547
548 MPS_DEBUG (s, "Compacting circles from %d to %d", i, j);
549
550 /* We shall now compact circles between i and j, so
551 * we start computing the mean of the radius */
552 rdpe_mul_eq_d (c.dradii[i],
553 c.partitioning[i + 1] - c.partitioning[i]);
554 for (k = i + 1; k <= j; k++)
555 {
556 rdpe_mul_d (tmp, c.dradii[k],
557 c.partitioning[k + 1] - c.partitioning[k]);
558 rdpe_add_eq (c.dradii[i], tmp);
559 }
560
561 rdpe_div_eq_d (c.dradii[i],
562 c.partitioning[j + 1] - c.partitioning[i]);
563
564 /* Move other circles backward */
565 for (k = j + 1; k < c.n_radii; k++)
566 {
567 rdpe_set (c.dradii[k - offset], c.dradii[k]);
568 c.partitioning[k - offset] = c.partitioning[k];
569 }
570
571 /* Set new c.n_radii and new partitioning */
572 c.n_radii = c.n_radii - offset;
573 c.partitioning[c.n_radii] = n;
574 }
575
576 free (h);
577 return c;
578 }
579
580 /**
581 * @brief Compute new starting approximations to the roots of the
582 * polynomial \f$p(x)\f$ having coefficients of modulus apoly, by
583 * means of the Rouche'-based criterion of Bini (Numer. Algo. 1996).
584 *
585 * The program can compute all the approximations
586 * (if \f$n\f$ is the degree of \f$p(x)\f$) or it may compute the
587 * approximations of the cluster of the <code>cluster_item</code>.
588 * The status vector is changed into <code>'o'</code> for the components
589 * that belong to a cluster with relative radius less than <code>eps</code>.
590 * The status vector is changed into <code>'f'</code> for the components
591 * that cannot be represented as <code>>dpe</code>.
592 *
593 * @param s mps_context struct pointer.
594 * @param n number of root in the cluster to consider
595 * @param cluster_item The item of the mps_clusterization containing the cluster
596 * that we are analyzing.
597 * @param clust_rad radius of the cluster.
598 * @param g new center in which the the polynomial will be shifted.
599 * @param eps maximum radius considered small enough to be negligible.
600 * @param dap[] moduli of the coefficients as <code>dpe</code> numbers.
601 */
602 MPS_PRIVATE void
mps_dstart(mps_context * s,int n,mps_cluster_item * cluster_item,rdpe_t clust_rad,rdpe_t g,rdpe_t eps,rdpe_t dap[])603 mps_dstart (mps_context * s, int n, mps_cluster_item * cluster_item,
604 rdpe_t clust_rad, rdpe_t g, rdpe_t eps, rdpe_t dap[])
605 {
606 int l = 0, i, j, jj, nzeros = 0;
607 rdpe_t r, tmp, tmp1;
608 double sigma, th, ang;
609 mps_boolean flag = false;
610
611 mps_cluster * cluster = NULL;
612 mps_root * root = NULL;
613
614 if (cluster_item)
615 cluster = cluster_item->cluster;
616
617 if (s->random_seed)
618 sigma = drand ();
619 else
620 {
621 /* If this is the first cluster select sigma = 0. In the other
622 * case try to maximize starting points distance. */
623 if (!cluster_item || cluster_item->prev == NULL)
624 {
625 sigma = s->last_sigma = MPS_STARTING_SIGMA;
626 }
627 else
628 {
629 sigma = mps_maximize_distance (s, s->last_sigma, cluster_item, n);
630 }
631 }
632
633 /* In the case of user-defined polynomial choose as starting
634 * approximations equispaced points in the unit circle. */
635 if (!MPS_IS_MONOMIAL_POLY (s->active_poly))
636 {
637 ang = pi2 / n;
638 for (i = 0; i < n; i++)
639 {
640 cdpe_set_d (s->root[i]->dvalue, cos (ang * i + sigma),
641 sin (ang * i + sigma));
642 }
643 return;
644 }
645
646
647 /* check if it is the case dpe_after_float, in this case set flag=true */
648 for (i = 0; i < n; i++)
649 {
650 flag = (s->root[i]->status == MPS_ROOT_STATUS_NOT_FLOAT);
651 if (flag)
652 break;
653 }
654
655 /* Compute starting radii with the Rouche based criterion */
656 mps_starting_configuration c =
657 mps_dcompute_starting_radii (s, n, cluster_item, clust_rad, g, eps, dap);
658
659 th = pi2 / n;
660
661 for (i = 0; i < c.n_radii; i++)
662 {
663 nzeros = c.partitioning[i + 1] - c.partitioning[i];
664 ang = pi2 / nzeros;
665 rdpe_set (r, c.dradii[i]);
666
667 if (cluster_item)
668 root = cluster->first;
669
670 for (j = c.partitioning[i]; j < c.partitioning[i + 1]; j++)
671 {
672 if (cluster_item)
673 {
674 l = root->k;
675 root = root->next;
676 }
677 else
678 l = j;
679
680
681 jj = j - c.partitioning[i];
682
683 /* If dpe_after_float (i.e., flag is true) recompute the starting
684 * values of only the approximations falling out of the range */
685 if (flag)
686 {
687 if (s->root[l]->status == MPS_ROOT_STATUS_NOT_FLOAT)
688 {
689 cdpe_set_d (s->root[l]->dvalue,
690 cos (ang * jj + th * c.partitioning[i + 1] +
691 sigma),
692 sin (ang * jj + th * c.partitioning[i + 1] +
693 sigma));
694 cdpe_mul_eq_e (s->root[l]->dvalue, r);
695 s->root[l]->status = MPS_ROOT_STATUS_CLUSTERED;
696 /*#G 27/4/98 if (rdpe_eq(r, big) || rdpe_eq(r, small)) */
697 if (rdpe_eq (r, RDPE_MIN) || rdpe_eq (r, RDPE_MAX))
698 s->root[l]->status = MPS_ROOT_STATUS_NOT_DPE;
699 }
700 }
701 else
702 {
703 /* else compute all the initial approximations */
704 cdpe_set_d (s->root[l]->dvalue,
705 cos (ang * jj + th * c.partitioning[i + 1] +
706 sigma),
707 sin (ang * jj + th * c.partitioning[i + 1] +
708 sigma));
709 cdpe_mul_eq_e (s->root[l]->dvalue, r);
710 /*#G 27/4/98 if (rdpe_eq(r, big) || rdpe_eq(r, small)) */
711 if (rdpe_eq (r, RDPE_MIN) || rdpe_eq (r, RDPE_MAX))
712 s->root[l]->status = MPS_ROOT_STATUS_NOT_DPE;
713 }
714
715 if (s->debug_level & MPS_DEBUG_APPROXIMATIONS)
716 {
717 MPS_DEBUG_CDPE (s, s->root[l]->dvalue, "s->droot[%d]", l);
718 }
719 }
720
721 /* If the new radius of the cluster is relatively small, then
722 * set the status component equal to 'o' (output) */
723 if (cluster_item)
724 {
725 rdpe_mul (tmp, g, eps);
726 rdpe_mul_d (tmp1, r, (double)nzeros);
727 if (rdpe_lt (tmp1, tmp))
728 {
729 for (root = cluster->first; root != NULL; root = root->next)
730 {
731 l = root->k;
732 s->root[l]->status = MPS_ROOT_STATUS_APPROXIMATED_IN_CLUSTER;
733 rdpe_set (s->root[l]->drad, tmp1);
734 }
735 }
736 }
737 }
738
739 mps_starting_configuration_clear (s, &c);
740 }
741
742 /**
743 * @brief Compute radii of the circles where the initial approximation
744 * will be disposed by mps_mstart()
745 *
746 * @param s mps_context* stuct pointer.
747 * @param n number of roots in the cluster.
748 * @param cluster_item The element of the <code>mps_clusterization</code> of which
749 * we are computing the starting points, or NULL if we are computing the starting
750 * points for all the approximations.
751 * @param clust_rad radius of the cluster.
752 * @param g new gravity center where the polynomial has been shifted.
753 * @param dap[] Array with the moduli of the coefficients.
754 *
755 * @see mps_mstart()
756 */
757 MPS_PRIVATE mps_starting_configuration
mps_mcompute_starting_radii(mps_context * s,int n,mps_cluster_item * cluster_item,rdpe_t clust_rad,rdpe_t g,rdpe_t dap[])758 mps_mcompute_starting_radii (mps_context * s, int n, mps_cluster_item * cluster_item,
759 rdpe_t clust_rad, rdpe_t g, rdpe_t dap[])
760 {
761 int i, offset, iold, nzeros, j, k;
762 rdpe_t big, small, tmp;
763 double xbig, xsmall, temp;
764
765 xsmall = rdpe_log (RDPE_MIN);
766 xbig = rdpe_log (RDPE_MAX);
767 rdpe_set (small, RDPE_MIN);
768 rdpe_set (big, RDPE_MAX);
769
770 mps_starting_configuration c = {
771 .n_radii = 0,
772 .partitioning = NULL,
773 .fradii = NULL,
774 .dradii = NULL
775 };
776
777 if (rdpe_eq (dap[0], rdpe_zero))
778 s->fap2[0] = -s->mpwp * LOG2;
779
780 /* check for possible null entries in the trailing coefficients */
781 for (i = 0; i <= n; i++)
782 if (rdpe_ne (dap[i], rdpe_zero))
783 s->fap2[i] = rdpe_log (dap[i]);
784 else
785 s->fap2[i] = s->fap2[0];
786
787 /* compute the convex hull */
788 int * h = mps_fconvex (s, n, s->fap2);
789
790 /* Count the number of vertexes in the convex hull. */
791 j = 0;
792 for (i = 1; i <= n; i++)
793 if (h[i])
794 j++;
795
796 c.dradii = mps_newv (rdpe_t, j + 2);
797 c.partitioning = mps_newv (int, j + 2);
798
799 /* Scan all the vertices of the convex hull */
800 c.partitioning[0] = 0;
801 c.n_radii = 0;
802 for (i = 1; i <= n; i++)
803 {
804 if (h[i])
805 {
806 iold = c.partitioning[c.n_radii];
807 nzeros = i - iold;
808 temp = (s->fap2[iold] - s->fap2[i]) / nzeros;
809 /* if the radius is too small or too big to be represented as dpe,
810 * output a warning message */
811 if (temp < xsmall)
812 {
813 rdpe_set (c.dradii[c.n_radii], small);
814 MPS_DEBUG (s, "Warning: Some zeros are too small to be\n"
815 "represented as cdpe, they are replaced by\n"
816 "small numbers and the status is set to 'F'.");
817 }
818 if (temp > xbig)
819 {
820 rdpe_set (c.dradii[c.n_radii], big);
821 MPS_DEBUG (s, "Warning: Some zeros are too big to be\n"
822 "represented as cdpe, they are replaced by\n"
823 "big numbers and the status is set to 'F'.");
824 }
825
826 /* if the radius is representable as dpe, compute it */
827 if (temp <= xbig && temp >= xsmall)
828 {
829 rdpe_set_d (c.dradii[c.n_radii], temp);
830 rdpe_exp_eq (c.dradii[c.n_radii]);
831 }
832
833 /* if the radius is greater than the radius of the cluster
834 * set the radius equal to the radius of the cluster */
835 if (rdpe_gt (c.dradii[c.n_radii], clust_rad))
836 rdpe_set (c.dradii[c.n_radii], clust_rad);
837
838 MPS_DEBUG_RDPE (s, c.dradii[c.n_radii], "c.dradii[%d]", c.n_radii);
839
840 /* Close partitioning and start a new one */
841 c.partitioning[++c.n_radii] = i;
842 }
843 }
844
845 /* Set last point of the partitioning */
846 c.partitioning[c.n_radii] = n;
847
848 /* Compact radius that are too near, but not if we are shifting. */
849 if (rdpe_ne (g, rdpe_zero))
850 return c;
851
852 for (i = 0; i < c.n_radii; i++)
853 {
854 /* Scan next radii to see if they are near the
855 * i-th that we are considering now */
856 for (j = i + 1; j < c.n_radii; j++)
857 {
858 /* Get an estimate of the distance between approximations that would
859 * be obtained collapsing the circles. It's a little understimated
860 * but this is due to the fact that changing approximation is only
861 * a fallback, and not the preferred action to perform. */
862 rdpe_sub (tmp, c.dradii[j], c.dradii[i]);
863 rdpe_abs_eq (tmp);
864 if (rdpe_lt (c.dradii[i], c.dradii[j]))
865 rdpe_div_eq (tmp, c.dradii[i]);
866 else
867 rdpe_div_eq (tmp, c.dradii[j]);
868 rdpe_div_eq_d (tmp, PI);
869 rdpe_mul_eq_d (tmp, c.partitioning[j + 1] - c.partitioning[i]);
870 if (rdpe_gt (tmp, rdpe_one))
871 {
872 break;
873 }
874 }
875
876 /* This is the number of circles that are near */
877 j--;
878 offset = j - i;
879
880 /* If there is nothing to compact, do not compact it */
881 if (offset == 0)
882 {
883 continue;
884 }
885
886 MPS_DEBUG (s, "Compacting circles from %d to %d", i, j);
887
888 /* We shall now compact circles between i and j, so
889 * we start computing the mean of the radius */
890 rdpe_mul_eq_d (c.dradii[i],
891 c.partitioning[i + 1] - c.partitioning[i]);
892 for (k = i + 1; k <= j; k++)
893 {
894 rdpe_mul_d (tmp, c.dradii[j],
895 c.partitioning[k + 1] - c.partitioning[k]);
896 rdpe_add_eq (c.dradii[i], tmp);
897 }
898
899 rdpe_div_eq_d (c.dradii[i],
900 c.partitioning[j + 1] - c.partitioning[i]);
901
902 /* Move other circles backward */
903 for (k = j + 1; k < c.n_radii; k++)
904 {
905 rdpe_set (c.dradii[k - offset], c.dradii[k]);
906 c.partitioning[k - offset] = c.partitioning[k];
907 }
908
909 /* Set new c.n_radii and new partitioning */
910 c.n_radii = c.n_radii - offset;
911 c.partitioning[c.n_radii] = n;
912 }
913
914 free (h);
915 return c;
916 }
917
918 /**
919 * @brief Multiprecision version of mps_fstart()
920 * @see mps_fstart()
921 */
922 MPS_PRIVATE void
mps_mstart(mps_context * s,int n,mps_cluster_item * cluster_item,rdpe_t clust_rad,rdpe_t g,rdpe_t dap[],mpc_t gg)923 mps_mstart (mps_context * s, int n, mps_cluster_item * cluster_item,
924 rdpe_t clust_rad,
925 rdpe_t g, rdpe_t dap[], mpc_t gg)
926 {
927 mps_cluster *cluster = NULL;
928 int i, j, jj, iold, l, nzeros;
929 double sigma, ang, th;
930 rdpe_t big, small, rtmp1, rtmp2;
931 cdpe_t ctmp;
932 mpc_t mtmp;
933 mps_boolean need_recomputing = true;
934 mps_root * root = NULL;
935
936 if (cluster_item)
937 cluster = cluster_item->cluster;
938
939 mpc_init2 (mtmp, s->mpwp);
940
941 rdpe_set (small, RDPE_MIN);
942 rdpe_set (big, RDPE_MAX);
943
944 if (s->random_seed)
945 sigma = drand ();
946 else
947 {
948 /* If this is the first cluster select sigma = 0. In the other
949 * case try to maximize starting points distance. */
950 if (!cluster_item || cluster_item->prev == NULL)
951 {
952 sigma = s->last_sigma = MPS_STARTING_SIGMA;
953 }
954 else
955 {
956 sigma = mps_maximize_distance (s, s->last_sigma, cluster_item, n);
957 }
958 }
959
960 nzeros = 0;
961
962 mps_starting_configuration c = MPS_STARTING_CONFIGURATION_INIT;
963
964 /* Continue to cycle while clusters are _real clusters_,
965 * because until now we have detached roots that were not
966 * certainly outside of the cluster. */
967 while (need_recomputing)
968 {
969 /* In the general case apply the Rouche-based criterion */
970 mps_starting_configuration_clear (s, &c);
971 c = mps_mcompute_starting_radii (s, n, cluster_item, clust_rad, g, dap);
972
973 /* We need to check that the points that we have kept out of
974 * the cluster are really out of the clusters. */
975 need_recomputing = false;
976
977 /* Find the maximum radii */
978 rdpe_set (rtmp1, rdpe_zero);
979 for (i = 0; i < c.n_radii; i++)
980 {
981 if (rdpe_lt (rtmp1, c.dradii[i]))
982 {
983 rdpe_set (rtmp1, c.dradii[i]);
984 }
985 }
986
987 mps_cluster_item * c_item;
988 for (c_item = s->clusterization->first; c_item != NULL; c_item = c_item->next)
989 {
990 if (c_item->detached && c_item->detached->cluster == cluster &&
991 cluster != NULL)
992 {
993 mps_root * root = c_item->cluster->first;
994 i = root->k;
995
996 /* Check if the root touches the cluster */
997 mpc_sub (mtmp, s->root[i]->mvalue, gg);
998 mpc_get_cdpe (ctmp, mtmp);
999 cdpe_mod (rtmp2, ctmp);
1000
1001 /* i is the cluster detached from i_clust, so the unique
1002 * root in it is the first. */
1003 rdpe_sub_eq (rtmp2, s->root[i]->drad);
1004 rdpe_sub_eq (rtmp2, rtmp1);
1005
1006 /* If they are too near we need to recompact them */
1007 if (rdpe_lt (rtmp2, rdpe_zero))
1008 {
1009 mps_cluster_item * next = c_item->next;
1010
1011 if (s->debug_level & MPS_DEBUG_CLUSTER)
1012 MPS_DEBUG (s, "Recompacting cluster with root %d", i);
1013
1014 need_recomputing = true;
1015 mps_clusterization_remove_cluster (s, s->clusterization, c_item);
1016 mps_cluster_insert_root (s, cluster, i);
1017
1018 c_item = next;
1019 }
1020 }
1021 }
1022 }
1023
1024 th = pi2 / n;
1025
1026 /* Set initial approximations accordingly to the computed
1027 * circles */
1028 root = (cluster) ? cluster->first : s->clusterization->first->cluster->first;
1029 for (i = 0; i < c.n_radii; i++)
1030 {
1031 mps_root * starting_root = root;
1032
1033 nzeros = c.partitioning[i + 1] - c.partitioning[i];
1034 ang = pi2 / nzeros;
1035 iold = c.partitioning[i];
1036
1037 /* Compute the initial approximations */
1038 for (j = iold; j < c.partitioning[i + 1]; j++)
1039 {
1040 jj = j - iold;
1041
1042 /* Take index relative to the cluster
1043 * that we are analyzing. */
1044 /* l = s->clust[s->punt[i_clust] + j]; */
1045 l = root->k;
1046
1047 cdpe_set_d (ctmp,
1048 cos (ang * jj + th * c.partitioning[i + 1] + sigma),
1049 sin (ang * jj + th * c.partitioning[i + 1] + sigma));
1050 cdpe_mul_eq_e (ctmp, c.dradii[i]);
1051 cdpe_set (s->root[l]->dvalue, ctmp);
1052
1053 if (rdpe_eq (c.dradii[i], big) || rdpe_eq (c.dradii[i], small))
1054 {
1055 s->root[l]->status = MPS_ROOT_STATUS_NOT_DPE;
1056 }
1057
1058 root = root->next;
1059 }
1060
1061
1062 /* If the new radius of the cluster is relatively small, then
1063 * set the status component equal to 'o' (output)
1064 * and set the corresponding radius */
1065 rdpe_set (rtmp1, c.dradii[i]);
1066 rdpe_mul_eq_d (rtmp1, (double)nzeros);
1067 if (rdpe_lt (rtmp1, clust_rad) && s->algorithm == MPS_ALGORITHM_SECULAR_GA)
1068 rdpe_set (rtmp1, clust_rad);
1069 rdpe_set (rtmp2, g);
1070 rdpe_mul_eq (rtmp2, s->eps_out);
1071 MPS_DEBUG (s, "Performing relatively small check");
1072 if (rdpe_le (rtmp1, rtmp2))
1073 {
1074 mps_root * root2 = starting_root;
1075 for (j = c.partitioning[i]; j < c.partitioning[i + 1]; j++)
1076 {
1077 l = root2->k;
1078 s->root[l]->status = MPS_ROOT_STATUS_APPROXIMATED_IN_CLUSTER;
1079 rdpe_mul_d (s->root[l]->drad, rtmp1, (double)nzeros);
1080 root2 = root2->next;
1081 }
1082 }
1083 rdpe_set (clust_rad, c.dradii[i]);
1084 }
1085
1086 mpc_clear (mtmp);
1087 mps_starting_configuration_clear (s, &c);
1088 }
1089
1090 /**
1091 * @brief This function scans the existing clusters and selects the
1092 * ones where shift in the gravity center must be done.
1093 * Then computes the gravity center g, performs the shift of
1094 * the variable and compute new starting approximations in the
1095 * cluster (floating point version).
1096 *
1097 * @param s A pointer to the current mps_context.
1098 *
1099 * Shift in g is perfomed if the approximation is included in the search set
1100 * or its inclusion status has not been determined yet.
1101 *
1102 * To compute g, first compute the weighted mean (super center sc)
1103 * of the approximations in the cluster, where the weight are the
1104 * radii, then compute the radius (super radius sr) of the disk
1105 * centered in the super center containing all the disks of the cluster.
1106 * Apply few steps of Newton's iteration to the (m-1)-st derivative
1107 * of the polynomial starting from the super center and obtain the
1108 * point g where to shift the variable.
1109 * If g is outside the super disk of center sc and radius sr
1110 * output a warning message.
1111 */
1112 MPS_PRIVATE void
mps_frestart(mps_context * s)1113 mps_frestart (mps_context * s)
1114 {
1115 int k, j, l;
1116 double sr, sum, rtmp, rtmp1;
1117 cplx_t sc, corr, ctmp;
1118 mps_boolean tst;
1119 mps_monomial_poly *p = MPS_MONOMIAL_POLY (s->active_poly);
1120 mps_approximation * g = NULL;
1121
1122 s->operation = MPS_OPERATION_SHIFT;
1123
1124 /* Variables for cluster analysis */
1125 mps_cluster * cluster;
1126 mps_cluster_item * c_item;
1127 mps_root * root;
1128
1129 /* For user's polynomials skip the restart stage (not yet implemented) */
1130 if (!MPS_IS_MONOMIAL_POLY (s->active_poly))
1131 return;
1132
1133 /* scan the existing clusters and select the ones where shift in
1134 * the gravity center must be done. tst=true means do not perform shift */
1135 for (c_item = s->clusterization->first; c_item != NULL; c_item = c_item->next)
1136 { /* loop1: */
1137 cluster = c_item->cluster;
1138
1139 /* Continue if the root is isolated */
1140 if (cluster->n == 1)
1141 continue;
1142
1143 tst = true;
1144 for (root = cluster->first; root != NULL; root = root->next)
1145 { /* looptst : */
1146 l = root->k;
1147 if (!s->root[l]->again)
1148 goto loop1;
1149 if (s->output_config->goal == MPS_OUTPUT_GOAL_COUNT)
1150 {
1151 if (s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1152 s->root[l]->inclusion == MPS_ROOT_INCLUSION_UNKNOWN)
1153 {
1154 tst = false;
1155 break;
1156 }
1157 }
1158 else if ((s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1159 s->root[l]->inclusion == MPS_ROOT_INCLUSION_UNKNOWN)
1160 || (s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1161 s->root[l]->inclusion == MPS_ROOT_INCLUSION_IN))
1162 {
1163 tst = false;
1164 break;
1165 }
1166 } /* for */
1167 if (tst)
1168 goto loop1;
1169
1170 /* Compute super center sc and super radius sr */
1171 mps_fsrad (s, cluster, sc, &sr);
1172
1173 /* Check the relative width of the cluster
1174 * If it is greater than 1 do not shift
1175 * and set status(:1)='c' that means
1176 * keep iterating Aberth's step. */
1177 if (sr > cplx_mod (sc))
1178 {
1179 mps_root * r;
1180 for (r = cluster->first; r != NULL; r = r->next)
1181 s->root[r->k]->status = MPS_ROOT_STATUS_CLUSTERED;
1182 MPS_DEBUG (s, "Cluster rel. large: skip to the next component");
1183 goto loop1;
1184 }
1185
1186 /* Now check the Newton isolation of the cluster */
1187 mps_cluster_item * c_item2;
1188 mps_cluster * cluster2;
1189 int m;
1190 for (c_item2 = s->clusterization->first; c_item2 != NULL; c_item2 = c_item2->next)
1191 if (c_item2 != c_item)
1192 {
1193 cluster2 = c_item2->cluster;
1194 for (root = cluster2->first; root != NULL; root = root->next)
1195 {
1196 m = root->k;
1197 cplx_sub (ctmp, sc, s->root[m]->fvalue);
1198 rtmp = cplx_mod (ctmp);
1199 rtmp1 = (sr + s->root[m]->frad) * 5 * s->n;
1200 if (rtmp < rtmp1)
1201 {
1202 for (root = cluster->first; root != NULL; root = root->next)
1203 s->root[root->k]->status = MPS_ROOT_STATUS_CLUSTERED;
1204 MPS_DEBUG_WITH_INFO (s, "Cluster not Newton isolated: skip "
1205 "to the next component.");
1206 goto loop1;
1207 }
1208 }
1209 }
1210 /* Compute the coefficients of the derivative of p(x) having order
1211 * equal to the multiplicity of the cluster -1. */
1212 sum = 0.0;
1213 for (j = 0; j <= MPS_POLYNOMIAL (p)->degree; j++)
1214 {
1215 sum += cplx_mod (p->fpc[j]);
1216 cplx_set (p->fppc[j], p->fpc[j]);
1217 }
1218 for (j = 1; j < cluster->n; j++)
1219 {
1220 for (k = 0; k <= MPS_POLYNOMIAL (p)->degree - j; k++)
1221 cplx_mul_d (p->fppc[k], p->fppc[k + 1], (double)(k + 1));
1222 }
1223 for (j = 0; j < MPS_POLYNOMIAL (p)->degree - cluster->n + 2; j++)
1224 s->fap1[j] = cplx_mod (p->fppc[j]);
1225
1226 /* Apply at most max_newt_it steps of Newton's iterations
1227 * to the above derivative starting from the super center
1228 * of the cluster. */
1229
1230 g = mps_approximation_new (s);
1231 cplx_set (g->fvalue, sc);
1232
1233 for (j = 0; j < s->max_newt_it; j++)
1234 {
1235 mps_fnewton (s, MPS_POLYNOMIAL (p), g, corr);
1236 cplx_sub_eq (g->fvalue, corr);
1237 if (!g->again)
1238 {
1239 break;
1240 }
1241 }
1242
1243 if (j == s->max_newt_it)
1244 {
1245 MPS_DEBUG (s, "Exceeded maximum Newton iterations in frestart");
1246 mps_approximation_free (s, g);
1247 return;
1248 }
1249 cplx_sub (ctmp, sc, g->fvalue);
1250 if (cplx_mod (ctmp) > sr)
1251 {
1252 MPS_DEBUG (s, "The gravity center falls outside the cluster");
1253 mps_approximation_free (s, g);
1254 return;
1255 }
1256 /* Compute the coefficients of the shifted polynomial p(x+g)
1257 * and compute new starting approximations
1258 * First check if shift may cause overflow, in this case skip
1259 * the shift stage */
1260
1261 if (s->n * log (cplx_mod (g->fvalue)) + log (sum) > log (DBL_MAX))
1262 goto loop1;
1263 MPS_DEBUG_CALL (s, "mps_fshift");
1264 mps_fshift (s, cluster->n, c_item, sr, g->fvalue, s->eps_out);
1265 rtmp = cplx_mod (g->fvalue);
1266 rtmp *= DBL_EPSILON * 2;
1267 for (root = cluster->first; root != NULL; root = root->next)
1268 {
1269 l = root->k;
1270 /* Choose as new incl. radius 2*multiplicity*(radius of the circle) */
1271 s->root[l]->frad = 2 * cluster->n * cplx_mod (s->root[l]->fvalue);
1272 cplx_add_eq (s->root[l]->fvalue, g->fvalue);
1273 if (s->root[l]->frad < rtmp) /* DARIO* aggiunto 1/5/97 */
1274 s->root[l]->frad = rtmp;
1275 }
1276
1277 loop1:
1278 if (g)
1279 {
1280 mps_approximation_free (s, g);
1281 g = NULL;
1282 }
1283 ;
1284 }
1285 }
1286
1287 /**
1288 * @brief This function scans the existing clusters and selects the
1289 * ones where shift in the gravity center must be done.
1290 * Then computes the gravity center g, performs the shift of
1291 * the variable and compute new starting approximations in the
1292 * cluster (DPE version).
1293 *
1294 * @param s A pointer to the current mps_context.
1295 *
1296 * Shift in g is perfomed if the approximation is included in the search set
1297 * or its inclusion status has not been determined yet.
1298 *
1299 * To compute g, first compute the weighted mean (super center sc)
1300 * of the approximations in the cluster, where the weight are the
1301 * radii, then compute the radius (super radius sr) of the disk
1302 * centered in the super center containing all the disks of the cluster.
1303 * Apply few steps of Newton's iteration to the (m-1)-st derivative
1304 * of the polynomial starting from the super center and obtain the
1305 * point g where to shift the variable.
1306 * If g is outside the super disk of center sc and radius sr
1307 * output a warning message.
1308 */
1309 MPS_PRIVATE void
mps_drestart(mps_context * s)1310 mps_drestart (mps_context * s)
1311 {
1312 int k, j, l;
1313 rdpe_t sr, rtmp, rtmp1;
1314 cdpe_t sc, corr, ctmp;
1315 mps_boolean tst;
1316 mps_monomial_poly *p = MPS_MONOMIAL_POLY (s->active_poly);
1317
1318 /* Cluster related variables */
1319 mps_cluster_item * c_item;
1320 mps_cluster * cluster;
1321 mps_root * root;
1322 mps_approximation * g = NULL;
1323
1324 s->operation = MPS_OPERATION_SHIFT;
1325
1326 /* For user's polynomials skip the restart stage (not yet implemented) */
1327 if (!MPS_IS_MONOMIAL_POLY (s->active_poly))
1328 return;
1329
1330 for (c_item = s->clusterization->first; c_item != NULL; c_item = c_item->next)
1331 { /* loop1: */
1332 cluster = c_item->cluster;
1333
1334 if (cluster->n == 1)
1335 continue;
1336
1337 tst = true;
1338 for (root = cluster->first; root != NULL; root = root->next)
1339 { /* looptst: */
1340 l = root->k;
1341 if (!s->root[l]->again)
1342 goto loop1;
1343 if (s->output_config->goal == MPS_OUTPUT_GOAL_COUNT)
1344 {
1345 if (s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1346 s->root[l]->inclusion == MPS_ROOT_INCLUSION_UNKNOWN)
1347 {
1348 tst = false;
1349 break;
1350 }
1351 }
1352 else if ((s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1353 s->root[l]->inclusion == MPS_ROOT_INCLUSION_UNKNOWN) ||
1354 (s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1355 s->root[l]->inclusion == MPS_ROOT_INCLUSION_IN))
1356 {
1357 tst = false;
1358 break;
1359 }
1360 } /* for */
1361 if (tst)
1362 goto loop1;
1363
1364 /* Compute super center sc and super radius sr */
1365 mps_dsrad (s, cluster, sc, sr);
1366
1367 /* Check the relative width of the cluster
1368 * If it is greater than 1 do not shift
1369 * and set statu(:1)='c' that means
1370 * keep iterating Aberth's step. */
1371 cdpe_mod (rtmp, sc);
1372 if (rdpe_gt (sr, rtmp))
1373 {
1374 for (root = cluster->first; root != NULL; root = root->next)
1375 {
1376 s->root[root->k]->status = MPS_ROOT_STATUS_CLUSTERED;
1377 /* err(clust[j])=true */
1378 }
1379 MPS_DEBUG (s, "cluster rel. large: skip to the next component");
1380 goto loop1;
1381 }
1382
1383 /* Now check the Newton isolation of the cluster */
1384 mps_cluster_item * c_item2;
1385 mps_cluster * cluster2;
1386
1387 for (c_item2 = s->clusterization->first; c_item2 != NULL; c_item2 = c_item2->next)
1388 if (c_item2 != c_item)
1389 {
1390 cluster2 = c_item2->cluster;
1391 for (root = cluster2->first; root != NULL; root = root->next)
1392 {
1393 cdpe_sub (ctmp, sc, s->root[root->k]->dvalue);
1394 cdpe_mod (rtmp, ctmp);
1395 rdpe_add (rtmp1, sr, s->root[root->k]->drad);
1396 rdpe_mul_eq_d (rtmp1, 2.0 * s->n);
1397 if (rdpe_lt (rtmp, rtmp1))
1398 {
1399 for (root = cluster->first; root != NULL; root = root->next)
1400 s->root[root->k]->status = MPS_ROOT_STATUS_CLUSTERED;
1401 MPS_DEBUG (s, "Cluster not Newton isolated: skip to the next component.");
1402 goto loop1;
1403 }
1404 }
1405 }
1406
1407 /* Compute the coefficients of the derivative of p(x) having order
1408 * equal to the multiplicity of the cluster -1. */
1409 for (j = 0; j <= s->n; j++)
1410 cdpe_set (s->dpc2[j], p->dpc[j]);
1411 for (j = 1; j < cluster->n; j++)
1412 {
1413 for (k = 0; k <= s->n - j; k++)
1414 cdpe_mul_d (s->dpc2[k], s->dpc2[k + 1], (double)(k + 1));
1415 }
1416 for (j = 0; j < s->n - cluster->n + 2; j++)
1417 cdpe_mod (s->dap1[j], s->dpc2[j]);
1418
1419 /* Apply at most max_newt_it steps of Newton's iterations
1420 * to the above derivative starting from the super center
1421 * of the cluster. */
1422 g = mps_approximation_new (s);
1423 cdpe_set (g->dvalue, sc);
1424 for (j = 0; j < s->max_newt_it; j++)
1425 { /* loop_newt: */
1426 mps_dnewton (s, MPS_POLYNOMIAL (p), g, corr);
1427 cdpe_sub_eq (g->dvalue, corr);
1428 if (!g->again)
1429 {
1430 break;
1431 }
1432 }
1433 if (j == s->max_newt_it)
1434 {
1435 MPS_DEBUG (s, "Exceeded maximum Newton iterations in frestart");
1436 mps_approximation_free (s, g);
1437 return;
1438 }
1439 cdpe_sub (ctmp, sc, g->dvalue);
1440 cdpe_mod (rtmp, ctmp);
1441 if (rdpe_gt (rtmp, sr))
1442 {
1443 MPS_DEBUG (s, "The gravity center falls outside the cluster");
1444 mps_approximation_free (s, g);
1445 return;
1446 }
1447 /* Shift the variable and compute new approximations */
1448 MPS_DEBUG_CALL (s, "mps_dshift");
1449 mps_dshift (s, cluster->n, c_item, sr, g->dvalue, s->eps_out);
1450 cdpe_mod (rtmp, g->dvalue);
1451 rdpe_mul_eq_d (rtmp, DBL_EPSILON * 2);
1452 for (root = cluster->first; root != NULL; root = root->next)
1453 {
1454 l = root->k;
1455
1456 /* Choose as new incl. radius 2*multiplicity*(radius of the circle) */
1457 cdpe_mod (s->root[l]->drad, s->root[l]->dvalue);
1458 rdpe_mul_eq_d (s->root[l]->drad,
1459 (double)(2 * cluster->n));
1460 cdpe_add_eq (s->root[l]->dvalue, g->dvalue);
1461 if (rdpe_lt (s->root[l]->drad, rtmp))
1462 rdpe_set (s->root[l]->drad, rtmp);
1463 }
1464
1465 loop1:
1466 if (g)
1467 {
1468 mps_approximation_free (s, g);
1469 g = NULL;
1470 }
1471 ;
1472 }
1473 }
1474
1475 /**
1476 * @brief Walk through the clusterization to find clusters that were detached
1477 * from the given one. Check if they were detached correctly according to the
1478 * new inclusion radii computed by mps_mrestart().
1479 *
1480 * @param ctx The context in which MPSolve is running.
1481 * @param clusterization The clusterization where clusters should be looked up
1482 * @param cluster_item The item containing the cluster from which the other clusters
1483 * may have been detached.
1484 */
1485 static void
mps_cluster_check_detachment(mps_context * ctx,mps_clusterization * clusterization,mps_cluster_item * cluster_item)1486 mps_cluster_check_detachment (mps_context * ctx, mps_clusterization * clusterization,
1487 mps_cluster_item * cluster_item)
1488 {
1489 MPS_DEBUG_CALL (ctx, "mps_cluster_check_detachment");
1490
1491 mps_cluster_item * item = NULL;
1492 mps_cluster * cluster = cluster_item->cluster;
1493 rdpe_t radius, distance_abs;
1494 mpc_t center, distance;
1495
1496 mpc_init2 (center, ctx->mpwp);
1497 mpc_init2 (distance, ctx->mpwp);
1498
1499 mps_msrad (ctx, cluster, center, radius);
1500 rdpe_mul_eq_d (radius, 2.0 * ctx->n);
1501
1502 for (item = clusterization->first; item != NULL; item = item->next)
1503 {
1504 mps_cluster * detached_cluster = item->cluster;
1505 if (item->detached == cluster_item)
1506 {
1507 int k = detached_cluster->first->k;
1508 mpc_sub (distance, center, ctx->root[k]->mvalue);
1509 mpc_rmod (distance_abs, distance);
1510
1511 if (rdpe_lt (distance_abs, radius))
1512 {
1513 mps_cluster_item * next = item->next;
1514 if (ctx->debug_level & MPS_DEBUG_CLUSTER)
1515 MPS_DEBUG (ctx,
1516 "Cluster containing root %d has not been correctly detached, reattaching.", k);
1517
1518 mps_cluster_insert_root (ctx, cluster, k);
1519 mps_clusterization_remove_cluster (ctx, clusterization, item);
1520
1521 item = next;
1522 }
1523 else
1524 {
1525 if (ctx->debug_level & MPS_DEBUG_CLUSTER)
1526 MPS_DEBUG (ctx,
1527 "Cluster containing root %d was successfuly detached.", k);
1528
1529 /* We need to stop marking this cluster as detached, that means
1530 * "experimental" in this context. */
1531 item->detached = NULL;
1532 }
1533 }
1534 }
1535
1536 mpc_clear (center);
1537 mpc_clear (distance);
1538 }
1539
1540 static void
mps_cluster_reattach_all_detached_clusters(mps_context * ctx,mps_clusterization * clusterization,mps_cluster_item * cluster_item)1541 mps_cluster_reattach_all_detached_clusters (mps_context * ctx, mps_clusterization * clusterization,
1542 mps_cluster_item * cluster_item)
1543 {
1544 MPS_DEBUG_CALL (ctx, "mps_cluster_reattach_all_detached_clusters");
1545
1546 mps_cluster_item * item = NULL;
1547 mps_cluster * cluster = cluster_item->cluster;
1548
1549 /* Find cluster that have been detached from this and
1550 * attach them. */
1551 for (item = clusterization->first; item != NULL; )
1552 {
1553 /* Pre-cache the next cluster item se we don't ruin the
1554 * loop if we need to free the current one. */
1555 mps_cluster *d_cluster = item->cluster;
1556 mps_cluster_item * next = item->next;
1557
1558 if (item->detached == cluster_item)
1559 {
1560 if (ctx->debug_level & MPS_DEBUG_CLUSTER)
1561 MPS_DEBUG (ctx,
1562 "Reattaching root %ld to its original cluster", d_cluster->first->k);
1563
1564 mps_cluster_insert_root (ctx, cluster, d_cluster->first->k);
1565 mps_clusterization_remove_cluster (ctx, clusterization, item);
1566 }
1567
1568 item = next;
1569 }
1570 }
1571
1572 /**
1573 * @brief This function scans the existing clusters and selects the
1574 * ones where shift in the gravity center must be done.
1575 * Then computes the gravity center g, performs the shift of
1576 * the variable and compute new starting approximations in the
1577 * cluster (MP version).
1578 *
1579 * @param s A pointer to the current mps_context.
1580 *
1581 * Shift in g is perfomed if the approximation is included in the search set
1582 * or its inclusion status has not been determined yet.
1583 *
1584 * To compute g, first compute the weighted mean (super center sc)
1585 * of the approximations in the cluster, where the weight are the
1586 * radii, then compute the radius (super radius sr) of the disk
1587 * centered in the super center containing all the disks of the cluster.
1588 * Apply few steps of Newton's iteration to the (m-1)-st derivative
1589 * of the polynomial starting from the super center and obtain the
1590 * point g where to shift the variable.
1591 * If g is outside the super disk of center sc and radius sr
1592 * output a warning message.
1593 */
1594 MPS_PRIVATE void
mps_mrestart(mps_context * s)1595 mps_mrestart (mps_context * s)
1596 {
1597 mps_boolean tst;
1598 int j, k, l;
1599 rdpe_t sr, rtmp, rtmp1, rtmp2, ag;
1600 cdpe_t tmp;
1601 mpf_t rea, srmp;
1602 mpc_t sc, corr, temp;
1603 mps_monomial_poly* p = MPS_MONOMIAL_POLY (s->active_poly);
1604
1605 /* Variables for cluster iteration */
1606 mps_cluster_item * c_item;
1607 mps_cluster * cluster;
1608 mps_root * root;
1609 mps_approximation * g = NULL;
1610
1611 long int starting_wp = s->mpwp;
1612
1613 MPS_DEBUG_THIS_CALL (s);
1614
1615 s->operation = MPS_OPERATION_SHIFT;
1616
1617 /* For user's polynomials skip the restart stage (not yet implemented) */
1618 if (!MPS_IS_MONOMIAL_POLY (s->active_poly))
1619 {
1620 MPS_DEBUG_WITH_INFO (s, "Skipping restart on user polynomial");
1621 return;
1622 }
1623
1624 mpf_init2 (rea, s->mpwp);
1625 mpf_init2 (srmp, s->mpwp);
1626 mpc_init2 (sc, s->mpwp);
1627 mpc_init2 (corr, s->mpwp);
1628 mpc_init2 (temp, s->mpwp);
1629
1630 /* Try to detach quasi-convergent elements from the clusters */
1631 mps_clusterization_detach_clusters (s, s->clusterization);
1632
1633 k = 0;
1634 for (c_item = s->clusterization->first; c_item != NULL; c_item = c_item->next)
1635 k = MAX (k, c_item->cluster->n);
1636
1637 for (c_item = s->clusterization->first; c_item != NULL; c_item = c_item->next)
1638 {
1639 /* loop1: */
1640 cluster = c_item->cluster;
1641
1642 if (cluster->n == 1)
1643 continue;
1644
1645 tst = true;
1646 for (root = cluster->first; root != NULL; root = root->next)
1647 { /* looptst: */
1648 l = root->k;
1649
1650 if (s->output_config->goal == MPS_OUTPUT_GOAL_COUNT)
1651 {
1652 if (s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1653 s->root[l]->inclusion == MPS_ROOT_INCLUSION_UNKNOWN)
1654 {
1655 tst = false;
1656 break;
1657 }
1658 }
1659 else if ((s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1660 s->root[l]->inclusion == MPS_ROOT_INCLUSION_UNKNOWN) ||
1661 (s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1662 s->root[l]->inclusion == MPS_ROOT_INCLUSION_IN))
1663 {
1664 tst = false;
1665 break;
1666 }
1667 } /* for */
1668
1669 if (tst)
1670 goto clean_detached_cluster;
1671
1672 /* Compute super center sc and super radius sr */
1673 mps_msrad (s, cluster, sc, sr);
1674
1675 /* MPS_DEBUG (s, "Clust = %d", i); */
1676 MPS_DEBUG_MPC (s, 10, sc, "Super center");
1677 MPS_DEBUG_RDPE (s, sr, "Super radius");
1678
1679 /* Check the relative width of the cluster
1680 * If it is greater than 1 do not shift
1681 * and set status[:1)='c' that means
1682 * keep iterating Aberth's step.
1683 * Check also the Newton-isolation of the cluster */
1684
1685 mpc_get_cdpe (tmp, sc);
1686 cdpe_mod (rtmp, tmp);
1687
1688 if (s->DOLOG)
1689 {
1690 rdpe_div (rtmp2, sr, rtmp);
1691 }
1692 if (s->debug_level & MPS_DEBUG_CLUSTER)
1693 MPS_DEBUG_RDPE (s, rtmp2, "Relative width");
1694
1695 if (rdpe_gt (sr, rtmp))
1696 {
1697 for (root = cluster->first; root != NULL; root = root->next)
1698 s->root[root->k]->status = MPS_ROOT_STATUS_CLUSTERED;
1699 MPS_DEBUG (s, "Cluster relat. large: skip to the next component");
1700
1701 goto clean_detached_cluster;
1702 }
1703
1704 /* Now check the Newton isolation of the cluster */
1705 mps_cluster_item * c_item2;
1706 mps_cluster * cluster2;
1707
1708 rdpe_set (rtmp2, rdpe_zero);
1709 for (c_item2 = s->clusterization->first; c_item2 != NULL; c_item2 = c_item2->next)
1710 {
1711 if (c_item2 != c_item)
1712 {
1713 cluster2 = c_item2->cluster;
1714 for (root = cluster2->first; root != NULL; root = root->next)
1715 {
1716 mpc_sub (temp, sc, s->root[root->k]->mvalue);
1717 mpc_get_cdpe (tmp, temp);
1718 cdpe_mod (rtmp, tmp);
1719 rdpe_sub_eq (rtmp, s->root[root->k]->drad);
1720 rdpe_sub_eq (rtmp, sr);
1721 rdpe_inv_eq (rtmp);
1722 rdpe_add_eq (rtmp2, rtmp);
1723 }
1724 }
1725 }
1726 rdpe_mul_eq (rtmp2, sr);
1727 rdpe_set_d (rtmp1, 0.3);
1728
1729 if (rdpe_gt (rtmp2, rtmp1))
1730 {
1731 for (root = cluster->first; root != NULL; root = root->next)
1732 s->root[root->k]->status = MPS_ROOT_STATUS_CLUSTERED;
1733 MPS_DEBUG (s, "Cluster not Newton isolated: skip to the next component");
1734 goto clean_detached_cluster;
1735 }
1736
1737 mps_debug_cluster_structure (s);
1738
1739 /* Create the approximation used to find the derivative zero */
1740 g = mps_approximation_new (s);
1741 mpc_set (g->mvalue, sc);
1742
1743 /* Compute the coefficients of the derivative of p(x) having order
1744 * equal to the multiplicity of the cluster -1. */
1745 mps_monomial_poly * der = mps_monomial_poly_derive (s, p, cluster->n - 1, s->mpwp);
1746 rdpe_t epsilon;
1747 rdpe_t error;
1748
1749 rdpe_set (epsilon, rdpe_zero);
1750
1751 /* Evaluate the necessary precision to perform the Newton iterations. */
1752 do
1753 {
1754 long int prec = mps_monomial_poly_get_precision (s, der);
1755 mps_mhorner_with_error2 (s, der, g->mvalue, corr, error, prec);
1756
1757 mpc_rmod (rtmp, corr);
1758 rdpe_set_2dl (epsilon, 1.0, -s->mpwp);
1759 rdpe_mul_eq (epsilon, rtmp);
1760
1761 if (rdpe_gt (error, epsilon))
1762 {
1763 if (s->debug_level & MPS_DEBUG_CLUSTER)
1764 MPS_DEBUG (s, "Raising precision of the derivative to %ld bits", prec * 2);
1765 mps_polynomial_raise_data (s, MPS_POLYNOMIAL (der), prec * 2);
1766 mpc_set_prec (g->mvalue, prec * 2);
1767 }
1768 } while (rdpe_gt (error, epsilon) && mps_monomial_poly_get_precision (s, der) <= cluster->n * s->mpwp);
1769
1770 /* Apply at most max_newt_it steps of Newton's iterations */
1771 /* to the above derivative starting from the super center */
1772 /* of the cluster. */
1773
1774 if (s->debug_level & MPS_DEBUG_CLUSTER)
1775 MPS_DEBUG_MPC (s, 30, g->mvalue, "g before newton");
1776
1777 for (j = 0; j < s->max_newt_it; j++)
1778 { /* loop_newt: */
1779 mps_mnewton (s, MPS_POLYNOMIAL (der), g, corr,
1780 mpc_get_prec (g->mvalue));
1781 mpc_sub_eq (g->mvalue, corr);
1782
1783 if (s->debug_level & MPS_DEBUG_CLUSTER)
1784 {
1785 MPS_DEBUG_RDPE (s, g->drad, "Radius of the cluster");
1786 MPS_DEBUG_MPC (s, 100, g->mvalue, "Iteration %d on the derivative", j);
1787 }
1788
1789 mpc_rmod (ag, g->mvalue);
1790
1791 if (!g->again || rdpe_Esp (g->drad) < rdpe_Esp (ag) - s->mpwp)
1792 break;
1793 }
1794
1795 mps_polynomial_free (s, MPS_POLYNOMIAL (der));
1796
1797 if (s->debug_level & MPS_DEBUG_CLUSTER)
1798 MPS_DEBUG (s, "Performed %d Newton iterations", j + 1);
1799 if (j == s->max_newt_it)
1800 {
1801 if (s->debug_level & MPS_DEBUG_CLUSTER)
1802 MPS_DEBUG (s, "Exceeded maximum number of Newton iterations.");
1803
1804 goto clean_detached_cluster;
1805 }
1806
1807 mpc_sub (temp, sc, g->mvalue);
1808 mpc_get_cdpe (tmp, temp);
1809 cdpe_mod (rtmp, tmp);
1810 if (rdpe_gt (rtmp, sr))
1811 {
1812 if (s->debug_level & MPS_DEBUG_CLUSTER)
1813 MPS_DEBUG (s, "The gravity center falls outside the cluster");
1814 goto clean_detached_cluster;
1815 }
1816
1817 /* shift the variable and compute new approximations */
1818 MPS_DEBUG_CALL (s, "mps_mshift");
1819 for (root = cluster->first; root != NULL; root = root->next)
1820 {
1821 l = root->k;
1822 mpc_get_cdpe (s->root[l]->dvalue, s->root[l]->mvalue);
1823 }
1824
1825 /*#D perform shift only if the new computed sr is smaller than old*0.25 */
1826 /* if (s->algorithm == MPS_ALGORITHM_SECULAR_GA) */
1827 /* rdpe_set (rtmp, sr); */
1828 /* else */
1829 rdpe_mul_d (rtmp, sr, 0.25);
1830
1831 /*#D AGO99 Factors: 0.1 (MPS2.0), 0.5 (GIUGN98) */
1832 mps_mshift (s, c_item->cluster->n, c_item, sr, g->mvalue);
1833 if (rdpe_le (sr, rtmp))
1834 { /* Perform shift only if the new clust is smaller */
1835 mpc_get_cdpe (tmp, g->mvalue);
1836 cdpe_mod (rtmp, tmp);
1837 if (s->lastphase == mp_phase)
1838 rdpe_set_2dl (rtmp1, 1.0, -starting_wp);
1839 else
1840 rdpe_set_d (rtmp1, DBL_EPSILON);
1841 rdpe_mul_eq (rtmp, rtmp1);
1842 rdpe_mul_eq_d (rtmp, 2);
1843
1844 MPS_DEBUG_RDPE (s, rtmp, "rtmp");
1845 MPS_DEBUG_RDPE (s, rtmp1, "rtmp1");
1846
1847 for (root = cluster->first; root != NULL; root = root->next)
1848 {
1849 l = root->k;
1850 mpc_set_cdpe (s->root[l]->mvalue, s->root[l]->dvalue);
1851 mpc_add_eq (s->root[l]->mvalue, g->mvalue);
1852 cdpe_mod (rtmp1, s->root[l]->dvalue);
1853 rdpe_mul_d (s->root[l]->drad, rtmp1,
1854 2.0 * cluster->n);
1855 if (rdpe_lt (s->root[l]->drad, rtmp))
1856 rdpe_set (s->root[l]->drad, rtmp);
1857
1858 if (s->debug_level & MPS_DEBUG_CLUSTER)
1859 {
1860 MPS_DEBUG_MPC (s, s->mpwp, s->root[l]->mvalue, "Repositioned root %4d", l);
1861 MPS_DEBUG_RDPE (s, s->root[l]->drad, "Radius for root %4d", l);
1862 }
1863 }
1864
1865 /* Check if the clusters that have been detached from this have been
1866 * detached for a good reason. */
1867 mps_cluster_check_detachment (s, s->clusterization, c_item);
1868 }
1869 else
1870 {
1871 if (s->debug_level & MPS_DEBUG_CLUSTER)
1872 MPS_DEBUG (s, "DO NOT PERFORM RESTART, "
1873 "new radius of the cluster is larger");
1874
1875 goto clean_detached_cluster;
1876 }
1877
1878 clean_detached_cluster:
1879 mps_cluster_reattach_all_detached_clusters (s, s->clusterization, c_item);
1880
1881 if (g != NULL)
1882 {
1883 mps_approximation_free (s, g);
1884 g = NULL;
1885 }
1886 }
1887
1888 mpc_clear (temp);
1889 mpc_clear (corr);
1890 mpc_clear (sc);
1891 mpf_clear (srmp);
1892 mpf_clear (rea);
1893 }
1894
1895 /**
1896 * @brief Check if the clusters are Newton isolated, in a way that we can
1897 * apply the shift in the gravity center without problems deriving from
1898 * other roots.
1899 *
1900 * @param s A pointer to the current mps_context.
1901 */
1902 MPS_PRIVATE void
mps_mnewtis(mps_context * s)1903 mps_mnewtis (mps_context * s)
1904 {
1905 mps_boolean tst;
1906 int k, l;
1907 rdpe_t sr, rtmp, rtmp1;
1908 cdpe_t tmp;
1909 mpf_t rea, srmp;
1910 mpc_t sc, temp;
1911 rdpe_t rtmp2;
1912
1913 /* Cluster variables */
1914 mps_cluster_item * c_item;
1915 mps_cluster * cluster;
1916 mps_root * root;
1917
1918 /* For user's polynomials skip the restart stage (not yet implemented) */
1919 if (!MPS_IS_MONOMIAL_POLY (s->active_poly))
1920 return;
1921 mpf_init2 (rea, s->mpwp);
1922 mpf_init2 (srmp, s->mpwp);
1923 mpc_init2 (sc, s->mpwp);
1924 mpc_init2 (temp, s->mpwp);
1925
1926 k = 0;
1927 for (c_item = s->clusterization->first; c_item != NULL; c_item = c_item->next)
1928 k = MAX (k, c_item->cluster->n);
1929
1930 for (c_item = s->clusterization->first; c_item != NULL; c_item = c_item->next)
1931 { /* loop1: */
1932 cluster = c_item->cluster;
1933
1934 if (cluster->n == 1)
1935 continue;
1936
1937 tst = true;
1938
1939 for (root = cluster->first; root != NULL; root = root->next)
1940 { /* looptst: */
1941 l = root->k;
1942 if (!s->root[l]->again)
1943 goto loop1;
1944 if (s->output_config->goal == MPS_OUTPUT_GOAL_COUNT)
1945 {
1946 if (s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1947 s->root[l]->inclusion == MPS_ROOT_INCLUSION_UNKNOWN)
1948 {
1949 tst = false;
1950 break;
1951 }
1952 }
1953 else if ((s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1954 s->root[l]->inclusion == MPS_ROOT_INCLUSION_UNKNOWN) ||
1955 (s->root[l]->status == MPS_ROOT_STATUS_CLUSTERED &&
1956 s->root[l]->inclusion == MPS_ROOT_INCLUSION_IN))
1957 {
1958 tst = false;
1959 break;
1960 }
1961 } /* for */
1962 if (tst)
1963 goto loop1;
1964
1965 /* Compute super center sc and super radius sr */
1966 mpf_set_ui (srmp, 0);
1967 for (root = cluster->first; root != NULL; root = root->next)
1968 {
1969 l = root->k;
1970 mpf_set_rdpe (rea, s->root[l]->drad);
1971 mpf_add (srmp, srmp, rea);
1972 }
1973 mpc_set_ui (sc, 0, 0);
1974 for (root = cluster->first; root != NULL; root = root->next)
1975 {
1976 l = root->k;
1977 mpf_set_rdpe (rea, s->root[l]->drad);
1978 mpc_mul_f (temp, s->root[l]->mvalue, rea);
1979 mpc_add_eq (sc, temp);
1980 }
1981 mpc_div_eq_f (sc, srmp);
1982 rdpe_set (sr, rdpe_zero);
1983 for (root = cluster->first; root != NULL; root = root->next)
1984 {
1985 l = root->k;
1986 mpc_sub (temp, sc, s->root[l]->mvalue);
1987 mpc_get_cdpe (tmp, temp);
1988 cdpe_mod (rtmp, tmp);
1989 rdpe_add_eq (rtmp, s->root[l]->drad);
1990 if (rdpe_lt (sr, rtmp))
1991 rdpe_set (sr, rtmp);
1992 }
1993
1994 /* Check the relative width of the cluster
1995 * If it is greater than 1 do not shift
1996 * and set status[:1)='c' that means
1997 * keep iterating Aberth's step.
1998 * Check also the Newton-isolation of the cluster */
1999 mpc_get_cdpe (tmp, sc);
2000 cdpe_mod (rtmp, tmp);
2001 rdpe_div (rtmp2, sr, rtmp);
2002 if (rdpe_gt (sr, rtmp))
2003 {
2004 for (root = cluster->first; root != NULL; root = root->next)
2005 s->root[root->k]->status = MPS_ROOT_STATUS_CLUSTERED;
2006 MPS_DEBUG (s, "Custer relatively large: "
2007 "skip to the next compoent");
2008 goto loop1;
2009 }
2010
2011 /* Now check the Newton isolation of the cluster */
2012 mps_cluster_item * c_item2;
2013 mps_cluster * cluster2;
2014 rdpe_set (rtmp2, rdpe_zero);
2015
2016 for (c_item2 = s->clusterization->first; c_item2 != NULL; c_item2 = c_item2->next)
2017 {
2018 if (c_item2 != c_item)
2019 {
2020 cluster2 = c_item2->cluster;
2021 for (root = cluster2->first; root != NULL; root = root->next)
2022 {
2023 mpc_sub (temp, sc, s->root[root->k]->mvalue);
2024 mpc_get_cdpe (tmp, temp);
2025 cdpe_mod (rtmp, tmp);
2026 rdpe_sub_eq (rtmp, s->root[root->k]->drad);
2027 rdpe_sub_eq (rtmp, sr);
2028 rdpe_inv_eq (rtmp);
2029 rdpe_add_eq (rtmp2, rtmp);
2030 }
2031 }
2032 }
2033 rdpe_mul_eq (rtmp2, sr);
2034 rdpe_set_d (rtmp1, 0.3);
2035
2036 if (rdpe_gt (rtmp2, rtmp1))
2037 {
2038 for (root = cluster->first; root != NULL; root = root->next)
2039 s->root[root->k]->status = MPS_ROOT_STATUS_CLUSTERED;
2040 MPS_DEBUG (s, "Cluster not Newton isolated: "
2041 "skip to the next component");
2042 goto loop1;
2043 }
2044 s->newtis = 1;
2045
2046 loop1:
2047 ;
2048 }
2049
2050 mpc_clear (temp);
2051 mpc_clear (sc);
2052 mpf_clear (srmp);
2053 mpf_clear (rea);
2054 }
2055