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 /**
16 * @brief Evaluate a secular equation <code>sec</code> in the point <code>x</code>
17 *
18 * @param s The <code>mps_context</code> of the computation.
19 * @param p The secular equation to evaluate.
20 * @param x The point in which the secular equation must be evaluated.
21 * @param value The value of the secular equation in the pointer <code>x</code>.
22 */
23 mps_boolean
mps_secular_feval(mps_context * s,mps_polynomial * p,cplx_t x,cplx_t value)24 mps_secular_feval (mps_context * s, mps_polynomial * p, cplx_t x, cplx_t value)
25 {
26 mps_secular_equation * sec = MPS_SECULAR_EQUATION (p);
27 cplx_t ctmp;
28 int i;
29
30 cplx_set (value, cplx_zero);
31
32 for (i = 0; i < s->n; i++)
33 {
34 cplx_sub (ctmp, x, sec->bfpc[i]);
35 if (cplx_eq_zero (ctmp))
36 return false;
37 cplx_div (ctmp, sec->afpc[i], ctmp);
38 cplx_add_eq (value, ctmp);
39 }
40
41 cplx_sub_eq (value, cplx_one);
42 return true;
43 }
44
45 mps_boolean
mps_secular_feval_derivative(mps_context * s,mps_polynomial * p,cplx_t x,cplx_t value)46 mps_secular_feval_derivative (mps_context * s, mps_polynomial * p, cplx_t x, cplx_t value)
47 {
48 mps_secular_equation * sec = MPS_SECULAR_EQUATION (p);
49 cplx_t ctmp;
50 int i;
51
52 cplx_set (value, cplx_zero);
53
54 for (i = 0; i < s->n; i++)
55 {
56 /* Compute 1 / (x - b_i) */
57 cplx_sub (ctmp, x, sec->bfpc[i]);
58
59 if (cplx_eq_zero (ctmp))
60 return false;
61
62 cplx_inv_eq (ctmp);
63 cplx_mul_eq (ctmp, ctmp);
64
65 /* Compute a_i / (x - b_i) */
66 cplx_mul_eq (ctmp, sec->afpc[i]);
67
68 /* Sum to the secular eqation */
69 cplx_sub_eq (value, ctmp);
70 }
71
72 return true;
73 }
74
75 /**
76 * @brief Evaluate a secular equation <code>sec</code> in the point <code>x</code>,
77 * estimating the error on the evaluation.
78 *
79 * @param s The <code>mps_context</code> of the computation.
80 * @param p The secular equation to evaluate.
81 * @param x The point in which the secular equation must be evaluated.
82 * @param value The value of the secular equation in the pointer <code>x</code>.
83 * @param error The absolute error on the evaluation.
84 */
85 mps_boolean
mps_secular_feval_with_error(mps_context * s,mps_polynomial * p,cplx_t x,cplx_t value,double * error)86 mps_secular_feval_with_error (mps_context * s, mps_polynomial * p, cplx_t x, cplx_t value,
87 double * error)
88 {
89 mps_secular_equation * sec = MPS_SECULAR_EQUATION (p);
90 cplx_t ctmp;
91 int i;
92
93 *error = 0.0f;
94 cplx_set (value, cplx_zero);
95
96 for (i = 0; i < s->n; i++)
97 {
98 cplx_sub (ctmp, x, sec->bfpc[i]);
99
100 if (cplx_eq_zero (ctmp))
101 return false;
102
103 cplx_div (ctmp, sec->afpc[i], ctmp);
104 cplx_add_eq (value, ctmp);
105 *error += cplx_mod (ctmp) * (i + 2);
106 }
107
108 cplx_sub_eq (value, cplx_one);
109 *error += 1.0f;
110
111 *error *= 4.0f * DBL_EPSILON;
112
113 return true;
114 }
115
116 /**
117 * @brief Evaluate a secular equation <code>sec</code> in the point <code>x</code>
118 *
119 * @param s The <code>mps_context</code> of the computation.
120 * @param p The secular equation to evaluate.
121 * @param x The point in which the secular equation must be evaluated.
122 * @param value The value of the secular equation in the point <code>x</code>.
123 */
124 mps_boolean
mps_secular_deval(mps_context * s,mps_polynomial * p,cdpe_t x,cdpe_t value)125 mps_secular_deval (mps_context * s, mps_polynomial * p, cdpe_t x, cdpe_t value)
126 {
127 mps_secular_equation *sec = MPS_SECULAR_EQUATION (p);
128 cdpe_t ctmp;
129 int i;
130
131 cdpe_set (value, cdpe_zero);
132
133 for (i = 0; i < s->n; i++)
134 {
135 cdpe_sub (ctmp, x, sec->bdpc[i]);
136
137 if (cdpe_eq_zero (ctmp))
138 return false;
139
140 cdpe_div (ctmp, sec->adpc[i], ctmp);
141 cdpe_add_eq (value, ctmp);
142 }
143
144 cdpe_sub_eq (value, cdpe_one);
145
146 return true;
147 }
148
149 mps_boolean
mps_secular_deval_derivative(mps_context * s,mps_polynomial * p,cdpe_t x,cdpe_t value)150 mps_secular_deval_derivative (mps_context * s, mps_polynomial * p, cdpe_t x, cdpe_t value)
151 {
152 mps_secular_equation *sec = MPS_SECULAR_EQUATION (p);
153 cdpe_t ctmp;
154 int i;
155
156 cdpe_set (value, cdpe_zero);
157
158 for (i = 0; i < s->n; i++)
159 {
160 /* Compute 1 / (x - b_i) */
161 cdpe_sub (ctmp, x, sec->bdpc[i]);
162
163 if (cdpe_eq_zero (ctmp))
164 return false;
165
166 cdpe_inv_eq (ctmp);
167 cdpe_mul_eq (ctmp, ctmp);
168
169 /* Compute a_i / (x - b_i) */
170 cdpe_mul_eq (ctmp, sec->adpc[i]);
171
172 /* Sum to the secular eqation */
173 cdpe_sub_eq (value, ctmp);
174 }
175
176 return true;
177 }
178
179
180 /**
181 * @brief Evaluate a secular equation <code>sec</code> in the point <code>x</code>
182 *
183 * @param s The <code>mps_context</code> of the computation.
184 * @param p The secular equation to evaluate.
185 * @param x The point in which the secular equation must be evaluated.
186 * @param value The value of the secular equation in the point <code>x</code>.
187 * @param error A bound to the module of the relative error occurred in the computation.
188 */
189 mps_boolean
mps_secular_deval_with_error(mps_context * s,mps_polynomial * p,cdpe_t x,cdpe_t value,rdpe_t error)190 mps_secular_deval_with_error (mps_context * s, mps_polynomial * p,
191 cdpe_t x, cdpe_t value, rdpe_t error)
192 {
193 mps_secular_equation * sec = MPS_SECULAR_EQUATION (p);
194 cdpe_t ctmp;
195 rdpe_t rtmp;
196 int i;
197
198 cdpe_set (value, cdpe_zero);
199 rdpe_set (error, rdpe_zero);
200
201 for (i = 0; i < s->n; i++)
202 {
203 cdpe_sub (ctmp, x, sec->bdpc[i]);
204 if (cdpe_eq_zero (ctmp))
205 return false;
206 cdpe_div (ctmp, sec->adpc[i], ctmp);
207 cdpe_mod (rtmp, ctmp);
208 cdpe_add_eq (value, ctmp);
209 rdpe_mul_eq_d (rtmp, i + 2);
210 rdpe_add_eq (error, rtmp);
211 }
212
213 cdpe_sub_eq (value, cdpe_one);
214 rdpe_add_eq (error, rdpe_one);
215
216 rdpe_mul_eq_d (error, 4.0f * DBL_EPSILON);
217
218 return true;
219 }
220
221
222 /**
223 * @brief Evaluate a secular equation <code>sec</code> in the point <code>x</code>.
224 *
225 * @param s The <code>mps_context</code> of the computation.
226 * @param p The secular equation to evaluate.
227 * @param x The point in which the sceular equation must be evaluated.
228 * @param value The value of the secular equation in the point <code>x</code>.
229 */
230 mps_boolean
mps_secular_meval(mps_context * s,mps_polynomial * p,mpc_t x,mpc_t value)231 mps_secular_meval (mps_context * s, mps_polynomial * p, mpc_t x, mpc_t value)
232 {
233 mps_secular_equation * sec = MPS_SECULAR_EQUATION (p);
234 mps_boolean success = true;
235 mpc_t ctmp;
236 unsigned int wp = mpc_get_prec (x);
237 int i;
238
239 /* Lower the working precision in case of limited precision coefficients
240 * in the input polynomial. */
241 if (p->prec > 0 && p->prec < wp)
242 wp = p->prec;
243
244 mpc_init2 (ctmp, wp);
245 mpc_set_ui (value, 0U, 0U);
246
247 for (i = 0; i < s->n; ++i)
248 {
249 mpc_sub (ctmp, x, sec->bmpc[i]);
250 if (mpc_eq_zero (ctmp))
251 {
252 success = false;
253 goto cleanup;
254 }
255
256 mpc_div (ctmp, sec->ampc[i], ctmp);
257 mpc_add_eq (value, ctmp);
258 }
259
260 mpc_sub_eq_ui (value, 1U, 0U);
261
262 cleanup:
263 mpc_clear (ctmp);
264 return success;
265 }
266
267 /**
268 * @brief Evaluate a secular equation <code>sec</code> in the point <code>x</code>.
269 *
270 * @param s The <code>mps_context</code> of the computation.
271 * @param p The secular equation to evaluate.
272 * @param x The point in which the sceular equation must be evaluated.
273 * @param value The value of the secular equation in the point <code>x</code>.
274 * @param error A bound to the absolute value of the error introduced in the computation.
275 */
276 mps_boolean
mps_secular_meval_with_error(mps_context * s,mps_polynomial * p,mpc_t x,mpc_t value,rdpe_t error)277 mps_secular_meval_with_error (mps_context * s, mps_polynomial * p, mpc_t x, mpc_t value, rdpe_t error)
278 {
279 mps_secular_equation * sec = MPS_SECULAR_EQUATION (p);
280 mpc_t ctmp;
281 rdpe_t rtmp, ax;
282 cdpe_t cdtmp;
283 unsigned int wp = mpc_get_prec (x);
284 int i;
285
286 mps_boolean successful_evaluation = true;
287
288 /* Lower the working precision in case of limited precision coefficients
289 * in the input polynomial. */
290 if (p->prec > 0 && p->prec < wp)
291 wp = p->prec;
292
293 if (mpc_get_prec (sec->ampc[0]) < wp)
294 mps_polynomial_raise_data (s, p, wp);
295
296 mpc_init2 (ctmp, wp);
297 mpc_set_ui (value, 0U, 0U);
298 mpc_set_prec (value, wp);
299
300 /* Get |x| */
301 mpc_rmod (ax, x);
302
303 rdpe_set (error, rdpe_zero);
304
305 for (i = 0; i < s->n; i++)
306 {
307 mpc_sub (ctmp, x, sec->bmpc[i]);
308
309 if (mpc_eq_zero (ctmp))
310 {
311 successful_evaluation = false;
312 goto cleanup;
313 }
314
315 mpc_div (ctmp, sec->ampc[i], ctmp);
316 mpc_add_eq (value, ctmp);
317
318 mpc_get_cdpe (cdtmp, ctmp);
319 cdpe_mod (rtmp, cdtmp);
320 rdpe_mul_eq_d (rtmp, i + 2);
321 rdpe_add_eq (error, rtmp);
322 }
323
324 mpc_sub_eq_ui (value, 1U, 0U);
325 rdpe_add_eq (error, rdpe_one);
326
327 if (p->prec < wp)
328 rdpe_set_2dl (rtmp, 4.0, 1 - p->prec);
329 else
330 rdpe_set_2dl (rtmp, 4.0, 1 - (long int)wp);
331
332 rdpe_mul_eq (error, rtmp);
333
334 cleanup:
335
336 mpc_clear (ctmp);
337
338 return successful_evaluation;
339 }
340
341