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 * Leonardo Robol <leonardo.robol@unipi.it>
9 */
10
11 #include <mps/mps.h>
12
13 mps_boolean
mps_chebyshev_poly_meval(mps_context * ctx,mps_polynomial * poly,mpc_t x,mpc_t value,rdpe_t error)14 mps_chebyshev_poly_meval (mps_context * ctx, mps_polynomial * poly, mpc_t x, mpc_t value, rdpe_t error)
15 {
16 long int wp = mpc_get_prec (x);
17
18 /* Lower the working precision in case of limited precision coefficients
19 * in the input polynomial. */
20 if (poly->prec > 0 && poly->prec < wp)
21 wp = poly->prec;
22
23 mps_chebyshev_poly * cpoly = MPS_CHEBYSHEV_POLY (poly);
24 int i;
25
26 mpc_t t0, t1, ctmp, ctmp2;
27 rdpe_t ax, rtmp, rtmp2;
28
29 mpc_rmod (ax, x);
30 rdpe_set (error, rdpe_zero);
31
32 /* Make sure that we have sufficient precision to perform the computation */
33 mps_polynomial_raise_data (ctx, poly, wp);
34
35 mpc_init2 (t0, wp);
36 mpc_init2 (t1, wp);
37 mpc_init2 (ctmp, wp);
38 mpc_init2 (ctmp2, wp);
39
40 mpc_set (value, cpoly->mfpc[0]);
41 mpc_set_ui (t0, 1U, 0U);
42 if (poly->degree == 0)
43 {
44 return true;
45 }
46
47 mpc_set (t1, x);
48 mpc_mul (ctmp, cpoly->mfpc[1], x);
49 mpc_add_eq (value, ctmp);
50
51 mpc_rmod (rtmp, ctmp);
52 rdpe_add_eq (error, rtmp);
53
54 for (i = 2; i <= poly->degree; i++)
55 {
56 mpc_mul (ctmp, x, t1);
57 mpc_mul_eq_ui (ctmp, 2U);
58 mpc_rmod (rtmp, ctmp);
59 mpc_sub_eq (ctmp, t0);
60
61 mpc_rmod (rtmp2, t0);
62 rdpe_add_eq (rtmp, rtmp2);
63
64 mpc_mul (ctmp2, ctmp, cpoly->mfpc[i]);
65 mpc_add_eq (value, ctmp2);
66
67 rdpe_mul_eq (rtmp, ax);
68 rdpe_add_eq (error, rtmp);
69
70 mpc_set (t0, t1);
71 mpc_set (t1, ctmp);
72 }
73
74 mpc_clear (t0);
75 mpc_clear (t1);
76 mpc_clear (ctmp);
77 mpc_clear (ctmp2);
78
79 rdpe_set_2dl (rtmp, 2.0, -wp);
80 rdpe_mul_eq (error, rtmp);
81
82 return true;
83 }
84