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 #include <stdarg.h>
14 #include <stdio.h>
15 #include <string.h>
16 #include <mps/mps.h>
17 #include <ctype.h>
18 #include <locale.h>
19
20 #define ISZERO -1
21
22 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25
26 #ifndef __WINDOWS
27 #include <unistd.h>
28 #else
29 #include <io.h>
30 #endif
31
32 /**
33 * @brief Read the approximations from the file opened in the
34 * rtstr member of the mps_context.
35 *
36 * @param s A pointer to the current mps_context.
37 *
38 * NOTE: This function is not used anywhere in the code, atm.
39 */
40 void
mps_readroots(mps_context * s)41 mps_readroots (mps_context * s)
42 {
43 long digits;
44 int i, read_elements;
45
46 if (s->DOLOG)
47 fprintf (s->logstr, "Reading roots...\n");
48
49 read_elements = fscanf (s->rtstr, "%ld", &digits);
50 if (!read_elements)
51 {
52 mps_error (s, "Error while reading roots, aborting.");
53 }
54
55 /* precision setup code goes here */
56
57 for (i = 0; i < s->n; i++)
58 mpc_inp_str_u (s->root[i]->mvalue, s->rtstr, 10);
59 }
60
61 /**
62 * @brief Count the roots that are included in the search set, excluded
63 * form it, or have an undetermined inclusion state.
64 *
65 * @param s A pointer to the current mps_context.
66 */
67 void
mps_countroots(mps_context * s)68 mps_countroots (mps_context * s)
69 {
70 int k;
71
72 if (s->DOLOG)
73 fprintf (s->logstr, "Counting roots...\n");
74
75 s->count[0] = s->count[1] = s->count[2] = 0;
76
77 for (k = 0; k < s->n; k++)
78 switch (s->root[k]->inclusion)
79 {
80 case MPS_ROOT_INCLUSION_IN:
81 s->count[0]++;
82 break;
83
84 case MPS_ROOT_INCLUSION_OUT:
85 s->count[1]++;
86 break;
87
88 default:
89 s->count[2]++;
90 break;
91 }
92
93 if (s->output_config->search_set == MPS_SEARCH_SET_UNITARY_DISC_COMPL)
94 s->count[1] += s->zero_roots;
95 else
96 s->count[0] += s->zero_roots;
97 }
98
99 /**
100 * @brief Print a summary of the count of the roots that can be obtained
101 * through a call to mps_countroots() to the outstr member of the
102 * current mps_context.
103 *
104 * @param s A pointer to the current mps_context.
105 */
106 void
mps_outcount(mps_context * s)107 mps_outcount (mps_context * s)
108 {
109 mps_countroots (s);
110
111 fprintf (s->outstr, "%d roots are inside;\n", s->count[0]);
112 fprintf (s->outstr, "%d roots are outside;\n", s->count[1]);
113 fprintf (s->outstr, "%d roots are uncertain.\n", s->count[2]);
114 if (s->DOLOG)
115 {
116 fprintf (s->logstr, "%d roots are inside;\n", s->count[0]);
117 fprintf (s->logstr, "%d roots are outside;\n", s->count[1]);
118 fprintf (s->logstr, "%d roots are uncertain.\n", s->count[2]);
119 }
120 }
121
122 /**
123 * @brief Print a float to stdout (or whatever the output stream is
124 * atm) respecting the given options, and only with the significant
125 * digits.
126 *
127 * @param s A pointer to the current mps_context.
128 * @param f The float approximation that should be printed.
129 * @param rad The current inclusion radius for that approximation.
130 * @param out_digit The number of output digits required.
131 * @param sign The sign of the approximation.
132 */
133 MPS_PRIVATE void
mps_outfloat(mps_context * s,mpf_t f,rdpe_t rad,long out_digit,mps_boolean sign)134 mps_outfloat (mps_context * s, mpf_t f, rdpe_t rad, long out_digit,
135 mps_boolean sign)
136 {
137 mpf_t t;
138 rdpe_t r, ro;
139 double d;
140 long l, digit, true_digit;
141
142 if (s->output_config->format == MPS_OUTPUT_FORMAT_FULL)
143 {
144 mpf_init2 (t, mpf_get_prec (f));
145 mpf_set (t, f);
146 mpf_out_str (s->outstr, 10, 0, t);
147 mpf_clear (t);
148 return;
149 }
150
151 mpf_init2 (t, s->output_config->prec);
152
153 mpf_get_rdpe (ro, f);
154 if (s->output_config->format == MPS_OUTPUT_FORMAT_GNUPLOT ||
155 s->output_config->format == MPS_OUTPUT_FORMAT_GNUPLOT_FULL)
156 rdpe_out_str_u (s->outstr, ro);
157 else
158 {
159 rdpe_abs_eq (ro);
160 if (rdpe_ne (ro, rdpe_zero))
161 rdpe_div (r, rad, ro);
162 else
163 rdpe_set_d (r, 1.0e-10);
164 digit = (long)(-rdpe_log10 (r) + 1.5);
165 if (digit <= 0)
166 {
167 rdpe_get_dl (&d, &l, ro);
168 fprintf (s->outstr, "0.e%ld", l);
169 }
170 else
171 {
172 true_digit = (long)(LOG10_2 * mpf_get_prec (f)) + 1;
173 true_digit = MIN (digit, true_digit);
174 true_digit = MIN (true_digit, out_digit);
175 if (sign)
176 mpf_set (t, f);
177 else
178 mpf_abs (t, f);
179 mpf_out_str (s->outstr, 10, true_digit, t);
180 }
181 }
182
183 mpf_clear (t);
184 }
185
186 /**
187 * @brief Print an approximation to stdout (or whatever the output
188 * stream currently selected in the mps_context is).
189 *
190 * @param s A pointer to the current mps_context.
191 * @param i The index of the approxiomation that shall be printed.
192 * @param num The number of zero roots.
193 */
194 MPS_PRIVATE void
mps_outroot(mps_context * s,int i,int num)195 mps_outroot (mps_context * s, int i, int num)
196 {
197 long out_digit;
198
199 out_digit = (long)(LOG10_2 * s->output_config->prec) + 10;
200
201 /* print format header */
202 switch (s->output_config->format)
203 {
204 case MPS_OUTPUT_FORMAT_COMPACT:
205 case MPS_OUTPUT_FORMAT_FULL:
206 fprintf (s->outstr, "(");
207 break;
208
209 case MPS_OUTPUT_FORMAT_VERBOSE:
210 fprintf (s->outstr, "Root(%d) = ", num);
211 break;
212
213 default:
214 break;
215 }
216
217 /* print real part */
218 if (i == ISZERO || s->root[i]->attrs == MPS_ROOT_ATTRS_IMAG)
219 fprintf (s->outstr, "0");
220 else
221 mps_outfloat (s, mpc_Re (s->root[i]->mvalue), s->root[i]->drad, out_digit, true);
222
223 /* print format middle part */
224 switch (s->output_config->format)
225 {
226 case MPS_OUTPUT_FORMAT_BARE:
227 fprintf (s->outstr, " ");
228 break;
229
230 case MPS_OUTPUT_FORMAT_GNUPLOT:
231 case MPS_OUTPUT_FORMAT_GNUPLOT_FULL:
232 fprintf (s->outstr, "\t");
233 break;
234
235 case MPS_OUTPUT_FORMAT_COMPACT:
236 case MPS_OUTPUT_FORMAT_FULL:
237 fprintf (s->outstr, ", ");
238 break;
239
240 case MPS_OUTPUT_FORMAT_VERBOSE:
241 if (i == ISZERO || mpf_sgn (mpc_Im (s->root[i]->mvalue)) >= 0)
242 fprintf (s->outstr, " + I * ");
243 else
244 fprintf (s->outstr, " - I * ");
245 break;
246
247 default:
248 break;
249 }
250
251 /* print imaginary part */
252 if (i == ISZERO || s->root[i]->attrs == MPS_ROOT_ATTRS_REAL)
253 fprintf (s->outstr, "0");
254 else
255 mps_outfloat (s, mpc_Im (s->root[i]->mvalue), s->root[i]->drad, out_digit,
256 s->output_config->format != MPS_OUTPUT_FORMAT_VERBOSE);
257
258 /* If the output format is GNUPLOT_FORMAT_FULL, print out also the radius */
259 if (s->output_config->format == MPS_OUTPUT_FORMAT_GNUPLOT_FULL)
260 {
261 fprintf (s->outstr, "\t");
262 rdpe_out_str_u (s->outstr, s->root[i]->drad);
263 fprintf (s->outstr, "\t");
264 rdpe_out_str_u (s->outstr, s->root[i]->drad);
265 }
266
267 /* print format ending */
268 switch (s->output_config->format)
269 {
270 case MPS_OUTPUT_FORMAT_COMPACT:
271 fprintf (s->outstr, ")");
272 break;
273
274 case MPS_OUTPUT_FORMAT_FULL:
275 fprintf (s->outstr, ")\n");
276 if (i != ISZERO)
277 {
278 rdpe_outln_str (s->outstr, s->root[i]->drad);
279 fprintf (s->outstr, "Status: %s, %s, %s\n",
280 MPS_ROOT_STATUS_TO_STRING (s->root[i]->status),
281 MPS_ROOT_ATTRS_TO_STRING (s->root[i]->attrs),
282 MPS_ROOT_INCLUSION_TO_STRING (s->root[i]->inclusion));
283 }
284 else
285 fprintf (s->outstr, " 0\n ---\n");
286 break;
287
288 default:
289 break;
290 }
291 fprintf (s->outstr, "\n");
292
293 /* debug info */
294 if (s->DOLOG)
295 {
296 if (i == ISZERO)
297 fprintf (s->logstr, "zero root %-4d = 0", num);
298 else
299 {
300 fprintf (s->logstr, "Root %-4d = ", i);
301 mpc_out_str_2 (s->logstr, 10, 0, 0, s->root[i]->mvalue);
302 fprintf (s->logstr, "\n");
303 fprintf (s->logstr, " Radius = ");
304 rdpe_outln_str (s->logstr, s->root[i]->drad);
305 fprintf (s->logstr, " Prec = %ld\n",
306 (long)(mpc_get_prec (s->root[i]->mvalue) / LOG2_10));
307 fprintf (s->logstr, " Approximation = %s\n",
308 MPS_ROOT_STATUS_TO_STRING (s->root[i]->status));
309 fprintf (s->logstr, " Attributes = %s\n",
310 MPS_ROOT_ATTRS_TO_STRING (s->root[i]->attrs));
311 fprintf (s->logstr, " Inclusion = %s\n",
312 MPS_ROOT_INCLUSION_TO_STRING (s->root[i]->inclusion));
313 fprintf (s->logstr, "--------------------\n");
314 }
315 }
316 }
317
318 /**
319 * @brief Print the approximations to stdout (or whatever the output
320 * stream currently selected in the mps_context is).
321 *
322 * @param s A pointer to the current mps_context.
323 */
324 void
mps_output(mps_context * s)325 mps_output (mps_context * s)
326 {
327 int i, ind, num = 0;
328
329 if (s->DOLOG)
330 fprintf (s->logstr, "--------------------\n");
331
332 if (s->output_config->format != MPS_OUTPUT_FORMAT_GNUPLOT &&
333 s->output_config->format != MPS_OUTPUT_FORMAT_GNUPLOT_FULL)
334 {
335 if (s->over_max)
336 {
337 mps_warn (s, "Warning: Input precision has been reached during computation, "
338 "so not all the required digits may have been computed.");
339 }
340 }
341
342 /* Start with plotting instructions in the case of
343 * MPS_OUTPUT_GNUPLOT_FULL, so the output can be
344 * piped directly to gnuplot */
345 if (s->output_config->format == MPS_OUTPUT_FORMAT_GNUPLOT_FULL)
346 {
347 fprintf (s->outstr, "# MPSolve output for GNUPLOT\n");
348 fprintf (s->outstr, "# Make user that this output is piped into gnuplot using a command like\n");
349 fprintf (s->outstr, "# mpsolve -Ogf | gnuplot \n");
350 fprintf (s->outstr, "set pointsize 0.3\n");
351 fprintf (s->outstr, "plot '-' title 'Computed roots' with %s\n", s->gnuplot_format);
352 }
353
354 if (s->output_config->goal == MPS_OUTPUT_GOAL_COUNT)
355 mps_outcount (s);
356 else
357 {
358 if (s->output_config->search_set != MPS_SEARCH_SET_UNITARY_DISC_COMPL)
359 for (i = 0; i < s->zero_roots; i++)
360 mps_outroot (s, ISZERO, num++);
361 for (ind = 0; ind < s->n; ind++)
362 {
363 i = s->order[ind];
364 if (s->root[i]->inclusion == MPS_ROOT_INCLUSION_OUT)
365 continue;
366 mps_outroot (s, i, num++);
367 }
368 }
369
370 if (s->output_config->format == MPS_OUTPUT_FORMAT_GNUPLOT_FULL)
371 {
372 fprintf (s->outstr, "e\n");
373 fprintf (s->outstr, "pause mouse close\n");
374 fprintf (s->outstr, "# End of MPSolve GNUPLOT output. If you are seeing this maybe\n");
375 fprintf (s->outstr, "# you forgot to pipe the ***solve command into gnuplot?\n");
376 }
377 }
378
379 /**
380 * @brief Update the MP version of the roots to the latest and greatest approximations.
381 *
382 * @param s A pointer to the current mps_context.
383 */
384 MPS_PRIVATE void
mps_copy_roots(mps_context * s)385 mps_copy_roots (mps_context * s)
386 {
387 int i;
388
389 MPS_DEBUG_THIS_CALL (s);
390
391 switch (s->lastphase)
392 {
393 case no_phase:
394 mps_error (s, "Nothing to copy");
395 break;
396
397 case float_phase:
398 if (s->DOSORT)
399 mps_fsort (s);
400 for (i = 0; i < s->n; i++)
401 {
402 mpc_set_prec (s->root[i]->mvalue, DBL_MANT_DIG);
403 mpc_set_cplx (s->root[i]->mvalue, s->root[i]->fvalue);
404 rdpe_set_d (s->root[i]->drad, s->root[i]->frad);
405 }
406 break;
407
408 case dpe_phase:
409 if (s->DOSORT)
410 mps_dsort (s);
411 for (i = 0; i < s->n; i++)
412 {
413 mpc_set_prec (s->root[i]->mvalue, DBL_MANT_DIG);
414 mpc_set_cdpe (s->root[i]->mvalue, s->root[i]->dvalue);
415 }
416 break;
417
418 case mp_phase:
419 if (s->DOSORT)
420 mps_msort (s);
421 break;
422 }
423 }
424
425 /**
426 * @brief Dump all the current approximation to the logstr selected
427 * in the current mps_context.
428 *
429 * @param s A pointer to the current mps_context.
430 *
431 * This function is tipically used when encountering some errors.
432 */
433 void
mps_dump(mps_context * s)434 mps_dump (mps_context * s)
435 {
436 int i;
437 FILE * dmpstr = s->logstr;
438
439 MPS_DEBUG (s, "Dumping the approximations:");
440
441 /* output current status */
442 /* fprintf (dmpstr, */
443 /* "Phase=%d, In=%d, Out=%d, Uncertain=%d, Zero=%d, Clusters=%ld\n", */
444 /* s->lastphase, s->count[0], s->count[1], s->count[2], s->zero_roots, */
445 /* s->clusterization->n); */
446
447 MPS_DEBUG (s, "Phase = %s, In = %d, Out = %d, Uncertain = %d, Zero = %d, Cluster = %ld",
448 MPS_PHASE_TO_STRING (s->lastphase), s->count[0], s->count[1], s->count[2],
449 s->zero_roots, s->clusterization->n);
450
451 /* output current approximations */
452 /* fprintf (dmpstr, "\nCurrent approximations:\n"); */
453 MPS_DEBUG (s, "Current approximations:");
454 for (i = 0; i < s->n; i++)
455 {
456 switch (s->lastphase)
457 {
458 case no_phase:
459 case float_phase:
460 MPS_DEBUG_CPLX (s, s->root[i]->fvalue, "Approximation %4d", i);
461 break;
462
463 case dpe_phase:
464 MPS_DEBUG_CDPE (s, s->root[i]->dvalue, "Approximation %4d", i);
465 break;
466
467 case mp_phase:
468 MPS_DEBUG_MPC (s, s->mpwp, s->root[i]->mvalue, "Approximation %4d", i);
469 break;
470 }
471 }
472
473 /* output radii */
474 MPS_DEBUG (s, "Current radii:");
475 for (i = 0; i < s->n; i++)
476 {
477 switch (s->lastphase)
478 {
479 case no_phase:
480 case float_phase:
481 MPS_DEBUG (s, "Radius of root %4d = %e", i, s->root[i]->frad);
482 break;
483
484 case dpe_phase:
485 case mp_phase:
486 MPS_DEBUG_RDPE (s, s->root[i]->drad, "Radius of root %4d", i);
487 break;
488 }
489 }
490
491 MPS_DEBUG (s, " ");
492 mps_dump_status (s, dmpstr);
493 }
494
495 /**
496 * @brief Dump cluster structure to <code>outstr</code>.
497 *
498 * @param s the mps_context struct pointer.
499 * @param outstr The output stream where the cluster structure
500 * will be dumped.
501 */
502 MPS_PRIVATE void
mps_dump_cluster_structure(mps_context * s,FILE * outstr)503 mps_dump_cluster_structure (mps_context * s, FILE * outstr)
504 {
505 fprintf (outstr,
506 " MPS_DUMP_CLUSTER_STRUCTURE: Dumping cluster structure\n");
507
508 mps_cluster_item * cluster_item;
509 mps_cluster * cluster;
510 mps_root * root;
511
512 for (cluster_item = s->clusterization->first; cluster_item != NULL;
513 cluster_item = cluster_item->next)
514 {
515 cluster = cluster_item->cluster;
516 fprintf (outstr, " Cluster contains %ld roots:\n", cluster->n);
517
518 /* Dump cluster roots, but not more than 15 for line, to make
519 * the output readable. */
520 int j = 0;
521 for (root = cluster->first; root != NULL; root = root->next)
522 {
523 /* Go to a newlint if 15 roots are printed out */
524 if (j % 15 == 0)
525 {
526 fprintf (outstr, "\n ");
527 }
528
529 fprintf (outstr, " %4ld", root->k);
530 j++;
531 }
532
533 /* Make space untile the next cluster */
534 fprintf (outstr, "\n\n");
535 }
536 }
537
538 /**
539 * @brief Dump status of all the root approximations
540 */
541 MPS_PRIVATE void
mps_dump_status(mps_context * s,FILE * outstr)542 mps_dump_status (mps_context * s, FILE * outstr)
543 {
544 int i;
545
546 MPS_DEBUG (s, " Approximation Attributes Inclusion");
547 for (i = 0; i < s->n; i++)
548 {
549 MPS_DEBUG (s, "Status %4d: %-25s %-15s %-15s", i,
550 MPS_ROOT_STATUS_TO_STRING (s->root[i]->status),
551 MPS_ROOT_ATTRS_TO_STRING (s->root[i]->attrs),
552 MPS_ROOT_INCLUSION_TO_STRING (s->root[i]->inclusion));
553 }
554 }
555
556 /**
557 * @brief Print a warning to the user.
558 *
559 * @param st The current mps_context
560 * @param format The printf-like format for the data to print
561 */
562 void
mps_warn(mps_context * st,char * format,...)563 mps_warn (mps_context * st, char *format, ...)
564 {
565 char * exclamation_mark = "";
566
567 va_list ap;
568
569 va_start (ap, format);
570
571 if (mps_is_a_tty (st->logstr))
572 exclamation_mark = "\033[33;1m!\033[0m";
573
574 if (st->DOWARN)
575 {
576 if (format[strlen (format)] == '\n')
577 {
578 fprintf (stderr, "%s ", exclamation_mark);
579 vfprintf (stderr, format, ap);
580 }
581 else
582 {
583 fprintf (stderr, "%s ", exclamation_mark);
584 vfprintf (stderr, format, ap);
585 fprintf (stderr, "\n");
586 }
587 }
588
589 va_end (ap);
590 }
591
592 /**
593 * @brief Check if the file descriptor associated to stream
594 * is bounded to a tty.
595 *
596 * @param stream the stream to check
597 */
598 MPS_PRIVATE mps_boolean
mps_is_a_tty(FILE * stream)599 mps_is_a_tty (FILE * stream)
600 {
601 #ifndef __WINDOWS
602 return isatty (fileno (stream));
603 #else
604 return _isatty (_fileno (stream));
605 #endif
606 }
607
608 /* Prepare a implicit definition if not provided by the compiler but
609 * available as a non-conformant extension. */
610 #ifndef vsnprintf
611 #ifdef HAVE_VSNPRINTF
612 int snprintf (char *str, size_t size, const char *format, ...);
613 #endif
614 #endif
615
616 /**
617 * @brief Record an error happened during the computation.
618 * This will set the internal status error to on, and the actual
619 * errors will printed with the first call to mps_print_errors().
620 * @param s A pointer to the current mps_context.
621 */
622 void
mps_error(mps_context * s,const char * format,...)623 mps_error (mps_context * s, const char * format, ...)
624 {
625 va_list ap;
626 int buffer_size = 32;
627 int missing_characters = 0;
628
629 va_start (ap, format);
630
631 s->error_state = true;
632 if (s->last_error == NULL)
633 s->last_error = mps_newv (char, buffer_size);
634
635 /* Measure space needed for the string, if our initial guess for the space neede
636 * is not enough */
637 while ((missing_characters = vsnprintf (s->last_error, buffer_size, format, ap)) > buffer_size)
638 {
639 buffer_size += missing_characters + 1;
640 s->last_error = mps_realloc (s->last_error, buffer_size);
641 }
642
643 va_end (ap);
644 }
645
646 /**
647 * @brief Print all the errors that have been recorded up to now.
648 * This function should be called only if mps_context_has_errors()
649 * returns true.
650 *
651 * @param s A pointer to the current mps_context.
652 */
653 void
mps_print_errors(mps_context * s)654 mps_print_errors (mps_context * s)
655 {
656 const char *error = s->last_error;
657 size_t length = strlen (error);
658
659 if (s->logstr == NULL)
660 s->logstr = stderr;
661
662 const char *exclamation_mark = "!";
663 if (mps_is_a_tty (s->logstr))
664 exclamation_mark = "\033[31;1m!\033[0m";
665
666 if (error[length] == '\n')
667 {
668 fprintf (stderr, "%s %s %s", exclamation_mark, "MPSolve encountered an error:", error);
669 }
670 else
671 {
672 fprintf (stderr, "%s %s %s\n", exclamation_mark, "MPSolve encountered an error:", error);
673 }
674
675 /* Dump approximations, but only if they are present */
676 if (s->root && s->lastphase)
677 mps_dump (s);
678 }
679