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 <mps/mps.h>
14
15 MPS_PRIVATE void
16 mps_mhorner_sparse (mps_context * s, mps_monomial_poly * p, mpc_t x, mpc_t value);
17
18 /**
19 * @brief Compute the value of the polynomial <code>p</code> in the point <code>x</code>
20 * and save it in <code>value</code>. If you need a bound to the relative error, try
21 * mps_mhorner_with_error().
22 *
23 * @param s The <code>mps_context</code> of the computation.
24 * @param p The <code>monomial_poly</code> to evaluate.
25 * @param x The point where the polynomial will be evaluated.
26 * @param value The multiprecision complex variable where the result will be stored.
27 */
28 MPS_PRIVATE void
mps_mhorner(mps_context * s,mps_monomial_poly * p,mpc_t x,mpc_t value)29 mps_mhorner (mps_context * s, mps_monomial_poly * p, mpc_t x, mpc_t value)
30 {
31 int j;
32
33 if (MPS_DENSITY_IS_SPARSE (s->active_poly->density))
34 {
35 mps_mhorner_sparse (s, p, x, value);
36 }
37 else
38 {
39 mps_with_lock (p->mfpc_mutex[MPS_POLYNOMIAL (p)->degree],
40 mpc_set (value, p->mfpc[MPS_POLYNOMIAL (p)->degree]);
41 );
42
43 for (j = MPS_POLYNOMIAL (p)->degree - 1; j >= 0; j--)
44 {
45 mpc_mul_eq (value, x);
46
47 pthread_mutex_lock (&p->mfpc_mutex[j]);
48 mpc_add_eq (value, p->mfpc[j]);
49 pthread_mutex_unlock (&p->mfpc_mutex[j]);
50 }
51 }
52 }
53
54 /**
55 * @brief Compute the value of the polynomial <code>p</code> in the point <code>x</code>
56 * and save it in <code>value</code>.
57 *
58 * A upper bound to the relative error of the evaluation will be stored in <code>relative_error</code>.
59 * The error is computed using the formula
60 * \f[
61 * n \frac{ap(|x|)}{|p(x)|} u
62 * \f]
63 * where \f$ap(x)\f$ is the polynomial with the coefficients equal to the moduli of the coefficients
64 * of \f$p(x)\f$.
65 *
66 * @param s The <code>mps_context</code> of the computation.
67 * @param p The <code>monomial_poly</code> to evaluate.
68 * @param x The point where the polynomial will be evaluated.
69 * @param value The multiprecision complex variable where the result will be stored.
70 * @param error The <code>RDPE</code> where the absolute error will be saved.
71 * @param wp The working precision to use for the computation. If this value is <code>0</code> then <code>s->mpwp</code>
72 * will be used.
73 */
74 void
mps_mhorner_with_error2(mps_context * s,mps_monomial_poly * p,mpc_t x,mpc_t value,rdpe_t error,long int wp)75 mps_mhorner_with_error2 (mps_context * s, mps_monomial_poly * p, mpc_t x, mpc_t value, rdpe_t error, long int wp)
76 {
77 int i;
78 rdpe_t apol, ax, u;
79 cdpe_t cx;
80
81 pthread_mutex_lock (&p->mfpc_mutex[0]);
82 if (mpc_get_prec (p->mfpc[0]) < wp)
83 {
84 pthread_mutex_unlock (&p->mfpc_mutex[0]);
85 mps_monomial_poly_raise_precision (s, MPS_POLYNOMIAL (p), wp);
86 }
87 else
88 pthread_mutex_unlock (&p->mfpc_mutex[0]);
89
90 if (mpc_get_prec (x) < wp)
91 mpc_set_prec (x, wp);
92
93 /* Set 4 * machine precision in u */
94 rdpe_set_2dl (u, 1.0, 2 - wp);
95
96 /* Compute the polynomial using horner */
97 mps_mhorner (s, p, x, value);
98
99 /* Compute ap(|x|) using horner */
100 mpc_get_cdpe (cx, x);
101 cdpe_mod (ax, cx);
102
103 rdpe_set (apol, p->dap[MPS_POLYNOMIAL (p)->degree]);
104 for (i = MPS_POLYNOMIAL (p)->degree - 1; i >= 0; i--)
105 {
106 rdpe_mul_eq (apol, ax);
107 rdpe_add_eq (apol, p->dap[i]);
108 }
109
110 /* Compute ap(|x|) / |p(x)| */
111 mpc_get_cdpe (cx, value);
112 cdpe_mod (ax, cx);
113
114 rdpe_set (error, apol);
115 rdpe_add_eq (error, ax);
116
117 rdpe_mul_eq (error, u);
118 }
119
120 /**
121 * @brief Compute the value of the polynomial <code>p</code> in the point <code>x</code>
122 * and save it in <code>value</code>.
123 *
124 * A upper bound to the relative error of the evaluation will be stored in <code>relative_error</code>.
125 *
126 * @param s The <code>mps_context</code> of the computation.
127 * @param p The <code>monomial_poly</code> to evaluate.
128 * @param x The point where the polynomial will be evaluated.
129 * @param value The multiprecision complex variable where the result will be stored.
130 * @param relative_error The <code>RDPE</code> where the relative error will be saved.
131 * @param wp The working precision to use for the computation. If this value is <code>0</code> then <code>s->mpwp</code>
132 * will be used.
133 */
134 MPS_PRIVATE void
mps_mhorner_with_error(mps_context * s,mps_monomial_poly * p,mpc_t x,mpc_t value,rdpe_t relative_error,long int wp)135 mps_mhorner_with_error (mps_context * s, mps_monomial_poly * p, mpc_t x, mpc_t value, rdpe_t relative_error, long int wp)
136 {
137 int j, my_wp;
138 mpc_t ss;
139
140 cdpe_t cdtmp;
141 rdpe_t r_ss, r_value, pol_on_ss;
142 rdpe_t my_eps, rtmp;
143
144 if (wp == 0)
145 my_wp = s->mpwp;
146 else
147 my_wp = wp;
148
149 /* Set up precision related variables */
150 rdpe_set_2dl (my_eps, 0.5, -my_wp);
151
152 /* Init multiprecision temporary values */
153 mpc_init2 (ss, my_wp);
154
155 rdpe_set (relative_error, rdpe_zero);
156
157 mpc_set (value, p->mfpc[MPS_POLYNOMIAL (p)->degree]);
158 for (j = MPS_POLYNOMIAL (p)->degree - 1; j >= 0; j--)
159 {
160 /* Normal horner computation */
161 mpc_mul (ss, value, x);
162 mpc_add_eq (ss, p->mfpc[j]);
163
164 /* Error estimate */
165 mpc_get_cdpe (cdtmp, ss);
166 cdpe_mod (r_ss, cdtmp);
167 mpc_get_cdpe (cdtmp, value);
168 cdpe_mod (r_value, cdtmp);
169
170 rdpe_div (pol_on_ss, r_value, r_ss);
171 rdpe_add (rtmp, relative_error, my_eps);
172 rdpe_mul_eq (rtmp, pol_on_ss);
173 rdpe_add_eq (relative_error, rtmp);
174
175 rdpe_div (rtmp, p->dap[j], r_ss);
176 rdpe_mul_eq (rtmp, my_eps);
177 rdpe_add_eq (relative_error, rtmp);
178
179 /* Update horner value */
180 mpc_set (value, ss);
181 }
182
183 mpc_clear (ss);
184 }
185
186 /**
187 * @brief Compute \f$p(z)\f$ by means of the Horner rule.
188 *
189 * @param st The pointer to the mps_context struct.
190 * @param n The degree of the polynomial.
191 * @param x The point in which the polynomial should be evaluated.
192 * @param p Vector of the coefficients of the polynomial.
193 * @param b Sparsity vector. If <code>b[i] == 0</code> then
194 * <code>p[i] == 0</code>.
195 * @param s RDPE pointer in which the result will be stored
196 * @param n_thread The index of the thread that is calling this
197 * routine. This is used to determine a safe memory are
198 * for temporary sparsity vectors.
199 */
200 MPS_PRIVATE void
mps_mhorner_sparse(mps_context * s,mps_monomial_poly * p,mpc_t x,mpc_t value)201 mps_mhorner_sparse (mps_context * s, mps_monomial_poly * p, mpc_t x,
202 mpc_t value)
203 {
204 int m, j, i, i1, i2, q;
205 mpc_t tmp, y;
206 mps_boolean bi;
207
208 /* Degree of the polynomial and sparsity vector */
209 mps_boolean * b = p->spar;
210 int n = MPS_POLYNOMIAL (p)->degree + 1;
211
212 mps_boolean *spar2 = mps_boolean_valloc (MPS_POLYNOMIAL (p)->degree + 2);
213 mpc_t *mfpc2 = mps_newv (mpc_t, MPS_POLYNOMIAL (p)->degree + 1);
214
215 long int wp;
216
217 pthread_mutex_lock (&p->mfpc_mutex[0]);
218 wp = mpc_get_prec (p->mfpc[0]);
219
220 /* Lower the working precision in case of limited precision coefficients
221 * in the input polynomial. */
222 if (MPS_POLYNOMIAL (p)->prec > 0 && MPS_POLYNOMIAL (p)->prec < wp)
223 wp = MPS_POLYNOMIAL (p)->prec;
224
225 mpc_vinit2 (mfpc2, MPS_POLYNOMIAL (p)->degree + 1, wp);
226 pthread_mutex_unlock (&p->mfpc_mutex[0]);
227
228 mpc_init2 (tmp, wp);
229 mpc_init2 (y, wp);
230
231 for (i = 0; i < n; i++)
232 spar2[i] = b[i];
233
234 for (i = 0; i < n; i++)
235 if (b[i])
236 {
237 pthread_mutex_lock (&p->mfpc_mutex[i]);
238 mpc_set (mfpc2[i], p->mfpc[i]);
239 pthread_mutex_unlock (&p->mfpc_mutex[i]);
240 }
241
242 q = mps_intlog2 (n + 1);
243 m = n;
244 mpc_set (y, x);
245 for (j = 0; j < q; j++)
246 {
247 spar2[m] = false;
248 m = (m + 1) >> 1;
249 for (i = 0; i < m; i++)
250 {
251 i2 = (i << 1) + 1;
252 i1 = i2 - 1;
253 bi = spar2[i1] || spar2[i2];
254 if (bi)
255 {
256 if (spar2[i1])
257 if (spar2[i2])
258 {
259 mpc_mul (tmp, y, mfpc2[i2]);
260 mpc_add (mfpc2[i], mfpc2[i1], tmp);
261 }
262 else
263 mpc_set (mfpc2[i], mfpc2[i1]);
264 else
265 mpc_mul (mfpc2[i], y, mfpc2[i2]);
266 }
267 spar2[i] = bi;
268 }
269 spar2[m] = false;
270 mpc_sqr_eq (y);
271 }
272 mpc_set (value, mfpc2[0]);
273
274 mpc_clear (y);
275 mpc_clear (tmp);
276
277 mpc_vclear (mfpc2, MPS_POLYNOMIAL (p)->degree + 1);
278 free (spar2);
279 free (mfpc2);
280 }
281
282 /**
283 * @brief Evaluate the polynomial p in the point x.
284 *
285 * @param s The <code>mps_context</code> of the computation.
286 * @param p The <code>mps_monomial_poly</code> to evaluate.
287 * @param x The point where the polynomial will be evaluated.
288 * @param value The value computed by the function.
289 */
290 void
mps_dhorner(mps_context * s,mps_monomial_poly * p,cdpe_t x,cdpe_t value)291 mps_dhorner (mps_context * s, mps_monomial_poly * p, cdpe_t x, cdpe_t value)
292 {
293 int j;
294
295 cdpe_set (value, p->dpc[MPS_POLYNOMIAL (p)->degree]);
296 for (j = MPS_POLYNOMIAL (p)->degree - 1; j >= 0; j--)
297 {
298 cdpe_mul_eq (value, x);
299 cdpe_add_eq (value, p->dpc[j]);
300 }
301 }
302
303
304
305 /**
306 * @brief Evaluate the polynomial p in the point x, and give also a bound to the
307 * relative error occured in the computation.
308 *
309 * @param s The <code>mps_context</code> of the computation.
310 * @param p The <code>mps_monomial_poly</code> to evaluate.
311 * @param x The point where the polynomial will be evaluated.
312 * @param value The value computed by the function.
313 * @param error A bound to the absolute error of the computation.
314 */
315 void
mps_dhorner_with_error(mps_context * s,mps_monomial_poly * p,cdpe_t x,cdpe_t value,rdpe_t error)316 mps_dhorner_with_error (mps_context * s, mps_monomial_poly * p, cdpe_t x, cdpe_t value, rdpe_t error)
317 {
318 rdpe_t ax;
319 int j;
320
321 mps_dhorner (s, p, x, value);
322
323 cdpe_mod (ax, x);
324 rdpe_set (error, p->dap[MPS_POLYNOMIAL (p)->degree]);
325 for (j = MPS_POLYNOMIAL (p)->degree - 1; j >= 0; j--)
326 {
327 rdpe_mul_eq (error, ax);
328 rdpe_add_eq (error, p->dap[j]);
329 }
330
331 rdpe_mul_eq_d (error, DBL_EPSILON);
332 }
333
334 /**
335 * @brief Evaluate the polynomial p in the point x.
336 *
337 * @param s The <code>mps_context</code> of the computation.
338 * @param p The <code>mps_monomial_poly</code> to evaluate.
339 * @param x The point where the polynomial will be evaluated.
340 * @param value The value computed by the function.
341 */
342 void
mps_fhorner(mps_context * s,mps_monomial_poly * p,cplx_t x,cplx_t value)343 mps_fhorner (mps_context * s, mps_monomial_poly * p, cplx_t x, cplx_t value)
344 {
345 int j;
346
347 cplx_set (value, p->fpc[MPS_POLYNOMIAL (p)->degree]);
348 for (j = MPS_POLYNOMIAL (p)->degree - 1; j >= 0; j--)
349 {
350 cplx_mul_eq (value, x);
351 cplx_add_eq (value, p->fpc[j]);
352 }
353 }
354
355 /**
356 * @brief Evaluate the polynomial p in the point x, and give also a bound to the
357 * relative error occured in the computation.
358 *
359 * @param s The <code>mps_context</code> of the computation.
360 * @param p The <code>mps_monomial_poly</code> to evaluate.
361 * @param x The point where the polynomial will be evaluated.
362 * @param error A pointer to the location when an upper bound to the computation
363 * error will be stored.
364 * @param value The value computed by the function.
365 */
366 void
mps_fhorner_with_error(mps_context * s,mps_monomial_poly * p,cplx_t x,cplx_t value,double * error)367 mps_fhorner_with_error (mps_context * s, mps_monomial_poly * p, cplx_t x, cplx_t value, double * error)
368 {
369 int j;
370 double ax = cplx_mod (x);
371
372 mps_fhorner (s, p, x, value);
373
374 *error = p->fap[MPS_POLYNOMIAL (p)->degree];
375 for (j = MPS_POLYNOMIAL (p)->degree - 1; j >= 0; j--)
376 {
377 *error = *error * ax + p->fap[j];
378 }
379
380 *error *= DBL_EPSILON;
381 }
382