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
14 /********************************************************
15 This file contains a sample of a user-defined polynomial.
16 In this sample the (Mandelbrot) polynomial is defined by
17 the relation
18 p_{i+1}(x)=x*p_i(x)^2+1 where p_0(x)=1, (1)
19 and its derivative by
20 p'_{i+1}(x)=p_i(x)^2+x*2*p_i(x)*p'_i(x).
21 The degree is 1, 3, 7, 15, 31,... (2^i-1).
22 The three programs below compute the values of p_i(x) and
23 p'_i(x), by means of the above formulae in the float version,
24 dpe version and mp version.
25 The programs compute also the auxiliary variable apeps needed
26 for testing the stop condition apeps>|p|, and for computing
27 the error bound rad (radius of the inclusion disk).
28 The formulae for the computation of apeps and rad are
29 obtained by means of a rounding error analysis of (1).
30 **********************************************************/
31
32 #include <float.h>
33 #include <math.h>
34 #include <mps/mps.h>
35
36 /**
37 * @brief User-defined program for the computation of \f$p\f$, \f$p'\f$.
38 *
39 * @param s The current mps_context
40 * @param poly The mps_polynomial being solved.
41 * @param root The approximation whose Newton correction shall be computed.
42 * @param corr The output value where the newton correction will be stored.
43 *
44 * This sample computes the 'Mandelbrot polynomial by
45 * means of the relation: p=1+x*p**2, starting with p=1
46 */
47 void
mps_fnewton_usr(mps_context * s,mps_polynomial * poly,mps_approximation * root,cplx_t corr)48 mps_fnewton_usr (mps_context * s, mps_polynomial * poly, mps_approximation * root, cplx_t corr)
49 {
50 cplx_t p, pp, pt, tmp;
51 double ap, ax, eps;
52 int i, m;
53
54 cplx_t x;
55
56 cplx_set (x, root->fvalue);
57 double * rad = &root->frad;
58 mps_boolean * again = &root->again;
59
60 m = (int)(log (s->n + 1.0) / LOG2);
61 if ((1 << m) <= s->n)
62 m++;
63 eps = (DBL_EPSILON * 4.0) * s->n;
64 ax = cplx_mod (x);
65
66 cplx_set (p, cplx_one);
67 cplx_set (pp, cplx_zero);
68 ap = 1.0;
69 for (i = 1; i <= m; i++)
70 {
71 cplx_sqr (tmp, p);
72 cplx_mul (pt, x, tmp);
73 cplx_add_eq (pt, cplx_one);
74 cplx_mul_eq (pp, x);
75 cplx_mul_eq (pp, p);
76 cplx_mul_eq_d (pp, 2.0);
77 cplx_add_eq (pp, tmp);
78 cplx_set (p, pt);
79 ap = ap * ax + cplx_mod (p);
80 }
81 ap = ap * ax;
82 cplx_div (corr, p, pp);
83
84 *again = cplx_mod (p) > eps * ap * 3;
85
86 *rad = s->n * (cplx_mod (p) + 3 * ap * eps) / cplx_mod (pp);
87 }
88
89 /**
90 * @brief User-defined program for the computation of \f$p\f$, \f$p'\f$.
91 *
92 * @param s The current mps_context
93 * @param poly The mps_polynomial being solved.
94 * @param root The approximation whose Newton correction shall be computed.
95 * @param corr The output value where the newton correction will be stored.
96 *
97 * This sample computes the 'Mandelbrot polynomial by
98 * means of the relation: p=1+x*p**2, starting with p=1
99 */
100 MPS_PRIVATE void
mps_dnewton_usr(mps_context * s,mps_polynomial * poly,mps_approximation * root,cdpe_t corr)101 mps_dnewton_usr (mps_context * s, mps_polynomial * poly, mps_approximation * root, cdpe_t corr)
102 {
103 cdpe_t p, pp, pt, tmp;
104 rdpe_t ap, ax, eps, temp, apeps, atmp;
105 int i, m;
106 cdpe_t x;
107
108 cdpe_set (x, root->dvalue);
109
110 m = (int)(log (s->n + 1.0) / LOG2);
111 if ((1 << m) <= s->n)
112 m++;
113 rdpe_set_d (eps, DBL_EPSILON);
114 rdpe_mul_eq_d (eps, (double)4 * s->n);
115 cdpe_mod (ax, x);
116
117 cdpe_set (p, cdpe_one);
118 cdpe_set (pp, cdpe_zero);
119 rdpe_set (ap, rdpe_one);
120 for (i = 1; i <= m; i++)
121 {
122 cdpe_sqr (tmp, p);
123 cdpe_mul (pt, x, tmp);
124 cdpe_add_eq (pt, cdpe_one);
125 cdpe_mul_eq (pp, x);
126 cdpe_mul_eq (pp, p);
127 cdpe_mul_eq_d (pp, 2.0);
128 cdpe_add_eq (pp, tmp);
129 cdpe_set (p, pt);
130 rdpe_mul_eq (ap, ax);
131 cdpe_mod (atmp, p);
132 rdpe_add_eq (ap, atmp);
133 }
134 rdpe_mul_eq (ap, ax);
135 cdpe_div (corr, p, pp);
136
137 cdpe_mod (temp, p);
138 rdpe_mul (apeps, ap, eps);
139 rdpe_mul_eq_d (apeps, 3.0);
140 root->again = rdpe_gt (temp, apeps);
141
142 rdpe_add (root->drad, temp, apeps);
143 rdpe_mul_eq_d (root->drad, (double)s->n);
144 cdpe_mod (temp, pp);
145 rdpe_div_eq (root->drad, temp);
146 if (rdpe_eq (root->drad, rdpe_zero))
147 rdpe_mul (root->drad, ax, eps);
148 }
149
150 /**
151 * @brief User-defined program for the computation of \f$p\f$, \f$p'\f$.
152 *
153 * @param s The current mps_context
154 * @param poly The mps_polynomial being solved.
155 * @param root The approximation whose Newton correction shall be computed.
156 * @param corr The output value where the newton correction will be stored.
157 *
158 * This sample computes the 'Mandelbrot polynomial by
159 * means of the relation: p=1+x*p**2, starting with p=1
160 */
161 MPS_PRIVATE void
mps_mnewton_usr(mps_context * s,mps_polynomial * poly,mps_approximation * root,mpc_t corr,long int wp)162 mps_mnewton_usr (mps_context * s, mps_polynomial * poly, mps_approximation * root, mpc_t corr, long int wp)
163 {
164 int i, m;
165 rdpe_t ap, ax, eps, temp, apeps, atmp;
166 cdpe_t ctmp;
167 mpc_t p, pp, pt, tmp;
168
169 mpc_init2 (p, s->mpwp);
170 mpc_init2 (pp, s->mpwp);
171 mpc_init2 (pt, s->mpwp);
172 mpc_init2 (tmp, s->mpwp);
173
174 m = (int)(log (s->n + 1.0) / LOG2);
175 if ((1 << m) <= s->n)
176 m++;
177 rdpe_set (eps, s->mp_epsilon);
178 rdpe_mul_eq_d (eps, (double)4 * s->n);
179 mpc_get_cdpe (ctmp, root->mvalue);
180 cdpe_mod (ax, ctmp);
181
182 mpc_set_ui (p, 1, 0);
183 mpc_set_ui (pp, 0, 0);
184 rdpe_set (ap, rdpe_one);
185 for (i = 1; i <= m; i++)
186 {
187 mpc_sqr (tmp, p);
188 mpc_mul (pt, root->mvalue, tmp);
189 mpc_add_eq_ui (pt, 1, 0);
190 mpc_mul_eq (pp, root->mvalue);
191 mpc_mul_eq (pp, p);
192 mpc_mul_eq_ui (pp, 2);
193 mpc_add_eq (pp, tmp);
194 mpc_set (p, pt);
195 rdpe_mul_eq (ap, ax);
196 mpc_get_cdpe (ctmp, p);
197 cdpe_mod (atmp, ctmp);
198 rdpe_add_eq (ap, atmp);
199 }
200 rdpe_mul_eq (ap, ax);
201 mpc_div (corr, p, pp);
202
203 mpc_get_cdpe (ctmp, p);
204 cdpe_mod (temp, ctmp);
205 rdpe_mul (apeps, ap, eps);
206 rdpe_mul_eq_d (apeps, 3.0);
207 root->again = rdpe_gt (temp, apeps);
208
209 rdpe_add (root->drad, temp, apeps);
210 rdpe_mul_eq_d (root->drad, (double)s->n);
211 mpc_get_cdpe (ctmp, pp);
212 cdpe_mod (temp, ctmp);
213 rdpe_div_eq (root->drad, temp);
214 if (rdpe_eq (root->drad, rdpe_zero))
215 rdpe_mul (root->drad, ax, eps);
216
217 mpc_clear (tmp);
218 mpc_clear (pt);
219 mpc_clear (pp);
220 mpc_clear (p);
221 }
222
223 mps_boolean
mps_feval_usr(mps_context * ctx,mps_polynomial * p,cplx_t x,cplx_t value,double * error)224 mps_feval_usr (mps_context * ctx, mps_polynomial * p, cplx_t x, cplx_t value, double * error)
225 {
226 int i;
227 int m = (int)(log (p->degree + 1.0) / LOG2);
228 double ax = cplx_mod (x);
229 cplx_t tmp;
230
231 if ((1 << m) <= p->degree)
232 m++;
233
234 cplx_set (value, cplx_one);
235
236 if (error)
237 *error = 0.0;
238
239 for (i = 1; i <= m; i++)
240 {
241 cplx_sqr (tmp, value);
242 cplx_mul (value, x, tmp);
243 cplx_add_eq (value, cplx_one);
244
245 if (error)
246 *error = *error * ax + cplx_mod (value);
247 }
248
249 if (error)
250 *error *= DBL_EPSILON;
251
252 return true;
253 }
254
255 mps_boolean
mps_deval_usr(mps_context * ctx,mps_polynomial * p,cdpe_t x,cdpe_t value,rdpe_t error)256 mps_deval_usr (mps_context * ctx, mps_polynomial * p, cdpe_t x, cdpe_t value, rdpe_t error)
257 {
258 int i;
259 int m = (int)(log (p->degree + 1.0) / LOG2);
260 rdpe_t ax, rtmp;
261 cdpe_t tmp;
262
263 if ((1 << m) <= p->degree)
264 m++;
265
266 cdpe_mod (ax, x);
267 cdpe_set (value, cdpe_one);
268 cdpe_mod (error, value);
269
270 for (i = 1; i <= m; i++)
271 {
272 cdpe_sqr (tmp, value);
273 cdpe_mul (value, x, tmp);
274 cdpe_add_eq (value, cdpe_one);
275
276 rdpe_mul_eq (error, ax);
277 cdpe_mod (rtmp, value);
278 rdpe_add_eq (error, rtmp);
279 }
280
281 rdpe_mul_eq_d (error, DBL_EPSILON);
282 return true;
283 }
284
285 mps_boolean
mps_meval_usr(mps_context * ctx,mps_polynomial * p,mpc_t x,mpc_t value,rdpe_t error)286 mps_meval_usr (mps_context * ctx, mps_polynomial * p, mpc_t x, mpc_t value, rdpe_t error)
287 {
288 int i;
289 int m = (int)(log (p->degree + 1.0) / LOG2);
290 rdpe_t ax, rtmp;
291 mpc_t tmp;
292 long int wp = mpc_get_prec (x);
293
294 /* Correct the working precision in case of a limited precision polynomial (quite unlikely
295 * in the Mandelbrot case, but still. */
296 if (p->prec > 0 && p->prec < wp)
297 wp = p->prec;
298
299 if ((1 << m) <= p->degree)
300 m++;
301
302 mpc_init2 (tmp, wp);
303
304 mpc_rmod (ax, x);
305 mpc_set_ui (value, 1U, 0U);
306 mpc_rmod (error, value);
307
308 for (i = 1; i <= m; i++)
309 {
310 mpc_sqr (tmp, value);
311 mpc_mul (value, x, tmp);
312 mpc_add_eq_ui (value, 1U, 0U);
313
314 rdpe_mul_eq (error, ax);
315 mpc_rmod (rtmp, value);
316 rdpe_add_eq (error, rtmp);
317 }
318
319 rdpe_set_2dl (rtmp, 1.0, -wp);
320 rdpe_mul_eq (error, rtmp);
321
322 mpc_clear (tmp);
323
324 return true;
325 }
326