1 /*
2 * This file is part of MPSolve 3.0
3 *
4 * Copyright (C) 2001-2013, Dipartimento di Matematica "L. Tonelli", Pisa.
5 * License: http://www.gnu.org/licenses/gpl.html GPL version 3 or higher
6 *
7 * Authors:
8 * Leonardo Robol <robol@mail.dm.unipi.it>
9 */
10
11 #include "quadratic-poly.h"
12 #include <stdlib.h>
13 #include <string.h>
14
15 extern char * starting_file;
16
17 void
mps_quadratic_poly_dstart(mps_context * ctx,mps_polynomial * p,mps_approximation ** approximations)18 mps_quadratic_poly_dstart (mps_context *ctx, mps_polynomial *p, mps_approximation ** approximations)
19 {
20 mps_quadratic_poly *mp = MPS_QUADRATIC_POLY (p);
21 int n = mps_context_get_degree (ctx);
22
23 if (mp->level <= 4)
24 {
25 mps_general_dstart (ctx, p, approximations);
26 }
27 else
28 {
29 int i;
30 cplx_t *roots = NULL;
31 cdpe_t rot;
32
33 cdpe_set_d (rot, 1 + p->degree * DBL_EPSILON, DBL_EPSILON);
34
35 mps_context * new_ctx = mps_context_new ();
36 mps_quadratic_poly *new_mp = mps_quadratic_poly_new (new_ctx, mp->level - 1);
37
38 mps_context_set_input_poly (new_ctx, MPS_POLYNOMIAL (new_mp));
39 mps_context_select_algorithm (new_ctx, MPS_ALGORITHM_SECULAR_GA);
40 mps_context_set_starting_phase (new_ctx, float_phase);
41 mps_context_set_avoid_multiprecision (ctx, true);
42 mps_mpsolve (new_ctx);
43
44 mps_context_get_roots_d (new_ctx, &roots, NULL);
45
46 for (i = 0; i < n - 1; i++)
47 {
48 cdpe_t starting_point;
49 int j = i / 2;
50
51 cdpe_set_x (starting_point, roots[j]);
52
53 if (i % 2 == 0)
54 cdpe_mul_eq (starting_point, rot);
55
56 mps_approximation_set_dvalue (ctx, approximations[i], starting_point);
57 }
58
59 mps_approximation_set_dvalue (ctx, approximations[n - 1], cdpe_one);
60
61 cplx_vfree (roots);
62 mps_quadratic_poly_free (new_ctx, MPS_POLYNOMIAL (new_mp));
63 mps_context_free (new_ctx);
64 }
65 }
66
67 void
mps_quadratic_poly_fstart(mps_context * ctx,mps_polynomial * p,mps_approximation ** approximations)68 mps_quadratic_poly_fstart (mps_context *ctx, mps_polynomial *p, mps_approximation ** approximations)
69 {
70 mps_quadratic_poly *mp = MPS_QUADRATIC_POLY (p);
71 int n = mps_context_get_degree (ctx);
72
73 mps_context * new_ctx = NULL;
74 mps_quadratic_poly * new_mp = NULL;
75
76 if (mp->level <= 7 && ! starting_file)
77 {
78 mps_general_fstart (ctx, p, approximations);
79 }
80 else
81 {
82 int i;
83 cplx_t *roots = cplx_valloc (p->degree / 2);
84 cplx_t rot_up, rot_down;
85
86 double eps = 1e-1 / p->degree;
87 if (eps < 4.0 * DBL_EPSILON)
88 eps = 4.0 * DBL_EPSILON;
89
90 cplx_set_d (rot_up, 1 + eps, 1e4 * DBL_EPSILON);
91 cplx_set_d (rot_down, 1 + eps, DBL_EPSILON);
92
93 if (starting_file != NULL)
94 {
95 char * line = NULL;
96 size_t length = 0;
97 char *p1, *p2;
98 FILE * input = fopen (starting_file, "r");
99
100 for (i = 0; i < pow (2, mp->level - 1) - 1; i++)
101 {
102 if (getline (&line, &length, input) == -1)
103 {
104 fprintf (stderr, "Failure to read a line\n");
105 return;
106 }
107
108 p1 = strstr (line, "(");
109 p2 = strstr (line, ",");
110 *p2 = ' ';
111
112 cplx_set_d (roots[i], strtod (p1 + 1, NULL), strtod (p2 + 1, NULL));
113 }
114
115 starting_file = NULL;
116 }
117 else
118 {
119 new_ctx = mps_context_new ();
120 new_mp = mps_quadratic_poly_new (new_ctx, mp->level - 1);
121
122 mps_context_set_input_poly (new_ctx, MPS_POLYNOMIAL (new_mp));
123 mps_context_select_algorithm (new_ctx, MPS_ALGORITHM_SECULAR_GA);
124 mps_context_set_starting_phase (new_ctx, float_phase);
125 mps_context_set_avoid_multiprecision (ctx, true);
126 mps_mpsolve (new_ctx);
127 mps_context_get_roots_d (new_ctx, &roots, NULL);
128 }
129
130 for (i = 0; i < n - 1; i++)
131 {
132 cplx_t starting_point;
133 int j = i / 2;
134 cplx_set (starting_point, roots[j]);
135
136 if (i % 2 == 0)
137 cplx_mul_eq (starting_point, (cplx_Im (roots[j]) > 0) ? rot_up : rot_down);
138
139 mps_approximation_set_fvalue (ctx, approximations[i], starting_point);
140 mps_approximation_set_frad (ctx, approximations[i], DBL_MAX);
141 }
142
143 if (mp->level % 2 == 1)
144 {
145 cplx_set_d (rot_up, -1.0, 0.0);
146 mps_approximation_set_fvalue (ctx, approximations[i], rot_up);
147 }
148 else
149 mps_approximation_set_fvalue (ctx, approximations[i], cplx_one);
150
151 cplx_vfree (roots);
152
153 if (new_ctx)
154 {
155 mps_quadratic_poly_free (new_ctx, MPS_POLYNOMIAL (new_mp));
156 mps_context_free (new_ctx);
157 }
158 }
159 }
160