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 * Giuseppe Fiorentino <fiorent@dm.unipi.it>
9 * Leonardo Robol <leonardo.robol@unipi.it>
10 */
11
12 #include <mps/mps.h>
13
14 static void
mps_raisetemp(mps_context * s,unsigned long int digits)15 mps_raisetemp (mps_context * s, unsigned long int digits)
16 {
17 int i;
18
19 for (i = 0; i <= s->n; i++)
20 {
21 mpc_set_prec (s->mfpc1[i], digits);
22 mpc_set_prec (s->mfppc1[i], digits);
23 }
24 }
25
26 /**
27 * @brief This routine computes the first \f$m+1\f$ coefficients of the shifted
28 * polynomial \f$p(x+g)\f$, by performing \f$m+1\f$ Horner divisions.
29 * This if the floating point version of this function.
30 *
31 * @param s The current mps_context.
32 * @param m The size of the cluster.
33 * @param mps_cluster_item A pointer to the cluster that shall be shifted.
34 * @param clust_rad A bound for the radius of the cluster.
35 * @param g The gravity center of the cluster.
36 * @parma eps The current value of epsilon that should be used as a treshold.
37 *
38 * Then it computes the new starting approximations for the
39 * cluster selected by applying mps_fstart() and by updating the approximations.
40 */
41 MPS_PRIVATE void
mps_fshift(mps_context * s,int m,mps_cluster_item * cluster_item,double clust_rad,cplx_t g,rdpe_t eps)42 mps_fshift (mps_context * s, int m, mps_cluster_item * cluster_item, double clust_rad,
43 cplx_t g, rdpe_t eps)
44 {
45 int i, j;
46 double ag;
47 cplx_t t;
48 mps_monomial_poly *p = MPS_MONOMIAL_POLY (s->active_poly);
49
50 /* Perform divisions */
51 ag = cplx_mod (g);
52 for (i = 0; i <= s->n; i++)
53 cplx_set (s->fppc1[i], p->fpc[i]);
54 for (i = 0; i <= m; i++)
55 {
56 cplx_set (t, s->fppc1[s->n]);
57 for (j = s->n - 1; j >= i; j--)
58 {
59 cplx_mul_eq (t, g);
60 cplx_add_eq (t, s->fppc1[j]);
61 cplx_set (s->fppc1[j], t);
62 }
63 cplx_set (p->fppc[i], t);
64 }
65
66 /* start */
67 for (i = 0; i <= m; i++)
68 s->fap1[i] = cplx_mod (p->fppc[i]);
69
70 /* If there is a custom starting point function use it, otherwise
71 * use the default one */
72 mps_fstart (s, m, cluster_item, clust_rad, ag, eps, s->fap1);
73 }
74
75 /**
76 * @brief This routine computes the first \f$m+1\f$ coefficients of the shifted
77 * polynomial \f$p(x+g)\f$, by performing \f$m+1\f$ Horner divisions.
78 * This if the DPE version of this function.
79 *
80 * @param s The current mps_context.
81 * @param m The size of the cluster.
82 * @param mps_cluster_item A pointer to the cluster that shall be shifted.
83 * @param clust_rad A bound for the radius of the cluster.
84 * @param g The gravity center of the cluster.
85 * @parma eps The current value of epsilon that should be used as a treshold.
86 *
87 * Then it computes the new starting approximations for the
88 * cluster selected by applying mps_fstart() and by updating the approximations.
89 */
90 MPS_PRIVATE void
mps_dshift(mps_context * s,int m,mps_cluster_item * cluster_item,rdpe_t clust_rad,cdpe_t g,rdpe_t eps)91 mps_dshift (mps_context * s, int m, mps_cluster_item * cluster_item, rdpe_t clust_rad,
92 cdpe_t g, rdpe_t eps)
93 {
94 int i, j;
95 rdpe_t ag;
96 cdpe_t t;
97 mps_monomial_poly * p = MPS_MONOMIAL_POLY (s->active_poly);
98
99 cdpe_mod (ag, g);
100 for (i = 0; i <= s->n; i++)
101 cdpe_set (s->dpc1[i], p->dpc[i]);
102 for (i = 0; i <= m; i++)
103 {
104 cdpe_set (t, s->dpc1[s->n]);
105 for (j = s->n - 1; j >= i; j--)
106 {
107 cdpe_mul_eq (t, g);
108 cdpe_add_eq (t, s->dpc1[j]);
109 cdpe_set (s->dpc1[j], t);
110 }
111 cdpe_set (s->dpc2[i], t);
112 }
113
114 /* start */
115 for (i = 0; i <= m; i++)
116 cdpe_mod (s->dap1[i], s->dpc2[i]);
117
118 mps_dstart (s, m, cluster_item, clust_rad, ag, eps, s->dap1);
119 }
120
121 /**
122 * @brief This routine computes the first \f$m+1\f$ coefficients of the shifted
123 * polynomial \f$p(x+g)\f$, by performing \f$m+1\f$ Horner divisions.
124 * This if the MP version of this function.
125 *
126 * @param s The current mps_context.
127 * @param m The size of the cluster.
128 * @param mps_cluster_item A pointer to the cluster that shall be shifted.
129 * @param clust_rad A bound for the radius of the cluster.
130 * @param g The gravity center of the cluster.
131 * @parma eps The current value of epsilon that should be used as a treshold.
132 *
133 * Then it computes the new starting approximations for the
134 * cluster selected by applying mps_fstart() and by updating the approximations.
135 */
136 MPS_PRIVATE void
mps_mshift(mps_context * s,int m,mps_cluster_item * cluster_item,rdpe_t clust_rad,mpc_t g)137 mps_mshift (mps_context * s, int m, mps_cluster_item * cluster_item, rdpe_t clust_rad, mpc_t g)
138 {
139 int i, j, k;
140 long int mpwp_temp, mpwp_max;
141 rdpe_t ag, ap, abp, as, mp_ep;
142 cdpe_t abd;
143 mpc_t t;
144 mps_monomial_poly *p = MPS_MONOMIAL_POLY (s->active_poly);
145
146 /* mps_cluster * cluster = cluster_item->cluster; */
147
148 mpc_init2 (t, s->mpwp);
149
150 /* Perform divisions
151 * In the mp version of the shift stage the computation
152 * is performed with increasing levels of working precision
153 * until the coefficients of the shifted polynomial have at
154 * least one correct bit. */
155 rdpe_set (mp_ep, s->mp_epsilon);
156 mpc_get_cdpe (abd, g);
157 cdpe_mod (ag, abd);
158 for (i = 0; i <= s->n; i++)
159 mpc_set (s->mfpc1[i], p->mfpc[i]);
160 rdpe_set (as, rdpe_zero);
161 rdpe_set (ap, rdpe_one);
162 mpc_set_ui (t, 0, 0);
163 k = 0;
164
165 /* store the current working precision mpnw into mpnw_tmp */
166 mpwp_temp = s->mpwp;
167 mpwp_max = m * s->mpwp;
168
169 do
170 { /* loop */
171 mpc_set (t, s->mfpc1[MPS_POLYNOMIAL (p)->degree]);
172 mpc_get_cdpe (abd, p->mfpc[s->n]);
173 cdpe_mod (ap, abd);
174 for (j = s->n - 1; j >= 0; j--)
175 {
176 mpc_get_cdpe (abd, p->mfpc[j]);
177 cdpe_mod (abp, abd);
178 rdpe_mul_eq (ap, ag);
179 rdpe_mul_eq_d (abp, (double)j);
180 rdpe_add_eq (ap, abp);
181 mpc_mul_eq (t, g);
182 mpc_add_eq (t, s->mfpc1[j]);
183 mpc_set (s->mfpc1[j], t);
184 }
185
186 mpc_set (s->mfppc1[0], t);
187 mpc_get_cdpe (abd, t);
188 cdpe_mod (as, abd);
189 rdpe_mul_eq (ap, mp_ep);
190 rdpe_mul_eq_d (ap, 4.0 * (s->n + 1));
191 k++;
192
193 if (rdpe_lt (as, ap))
194 {
195 mpwp_temp += s->mpwp;
196
197 /* if ((mpwp_temp > mpwp_max || mpwp_temp > s->output_config->prec * m * 2)) */
198 /* { */
199 /* MPS_DEBUG (s, "Reached the maximum allowed precision in mshift"); */
200 /* break; */
201 /* } */
202
203 rdpe_set_2dl (mp_ep, 1.0, 1 - mpwp_temp);
204 mps_raisetemp (s, mpwp_temp);
205 mpc_set_prec (t, (unsigned long int)mpwp_temp);
206 mpc_set_prec (g, (unsigned long int)mpwp_temp);
207 if (mpwp_max < mpwp_temp)
208 mpwp_max = mpwp_temp;
209
210 for (j = 0; j <= s->n; j++)
211 mpc_set (s->mfpc1[j], p->mfpc[j]);
212 }
213 } while (rdpe_lt (as, ap) && (k <= m)); /* loop */
214
215 mps_raisetemp (s, 1 * mpwp_temp);
216
217 for (i = 1; i <= m; i++)
218 {
219 /* mpwp_temp = MAX (mpwp_temp - s->mpwp, s->mpwp); */
220 /* mps_raisetemp (s, mpwp_temp); */
221 /* mpc_set_prec (t, (unsigned long int) mpwp_temp); */
222 /* mpc_set_prec (g, (unsigned long int) mpwp_temp); */
223 mpc_set (t, s->mfpc1[s->n]);
224
225 for (j = s->n - 1; j >= i; j--)
226 {
227 mpc_mul_eq (t, g);
228 mpc_add_eq (t, s->mfpc1[j]);
229 mpc_set (s->mfpc1[j], t);
230 }
231 mpc_set (s->mfppc1[i], t);
232 }
233 /*
234 raisetemp_raw(mpwp);
235 mpc_set_prec_raw(s, (unsigned long int) mpwp);
236 mpc_set_prec_raw(g, (unsigned long int) mpwp);
237
238 segue alternativa
239 */
240 /* mps_raisetemp (s, 2 * mpwp_max); */
241 /* mpc_set_prec (t, (unsigned long int) 2 * mpwp_max); */
242 /* mpc_set_prec (g, (unsigned long int) 2 * mpwp_max); */
243 mps_raisetemp (s, 2 * mpwp_temp);
244 mpc_set_prec (t, (unsigned long int)s->mpwp);
245 mpc_set_prec (g, (unsigned long int)s->mpwp);
246
247 if (rdpe_lt (as, ap))
248 {
249 for (j = 0; j < m; j++)
250 rdpe_set (s->dap1[j], ap);
251 mpc_get_cdpe (abd, s->mfppc1[m]);
252 cdpe_mod (s->dap1[m], abd);
253 }
254 else
255 for (i = 0; i <= m; i++)
256 {
257 mpc_get_cdpe (abd, s->mfppc1[i]);
258 cdpe_mod (s->dap1[i], abd);
259 }
260
261 /* Debug the coefficients of the shifted polynomial */
262 if (s->debug_level & MPS_DEBUG_CLUSTER)
263 for (i = 0; i <= m; i++)
264 MPS_DEBUG_MPC (s, mpc_get_prec (s->mfppc1[i]), s->mfppc1[i],
265 "P(x + g), coefficient of degree %d", i);
266
267 mps_mstart (s, m, cluster_item, clust_rad, ag, s->dap1, g);
268
269 mpc_clear (t);
270 }
271