## Name:

interpolate computes an interpolation polynomial.

## Library names:

sollya_obj_t sollya_lib_interpolate(sollya_obj_t, sollya_obj_t, ...) sollya_obj_t sollya_lib_v_interpolate(sollya_obj_t, sollya_obj_t, va_list)

## Usage:

interpolate(X, Y) : (list, list) -> function interpolate(X, Y, D) : (list, list, constant) -> function interpolate(X, Y, D, I) : (list, list, constant, range) -> function interpolate(X, f) : (list, function) -> function interpolate(X, f, D) : (list, function, constant) -> function interpolate(X, f, D, I) : (list, function, constant, range) -> function

## Parameters:

• X is a list of constants for the interpolation abscissas
• Y is a list of constants for the interpolation ordinates
• f is a function to be interpolated
• D is a constant defining the maximum permissible error
• I is an interval defining the domain to be considered for the error

## Description:

• interpolate computes an interpolation polynomial p such that for all given abscissa points xi in the list X the interpolation condition is fullfilled with p(xi) = yi, where yi is the corresponding entry in the list of ordinates Y. All entries of the lists X and Y need to be constant expressions; the entries in the list of abscissas X need to be pairwise distinct. If no parameter D is given, the interpolation is exact, i.e. an interpolation polynomial (expression) is formed such that the interpolation condition is satisfied without error.
• interpolate accepts a function f instead of a list of ordinates. When called in this manner, interpolate computes a polynomial p such that for all given abscissa points xi in the list X the interpolation condition is fullfilled with p(xi) = f(xi). All entries in the list X still need to be constant expressions and they must be pairwise distinct. If no parameter D is given, the interpolation stays exact, even if this means forming constant expressions for the ordinates f(xi).
• As exact interpolation is often not needed and in order to speed up computations, interpolate accepts an optional parameter D, a positive constant, that defines the maximum permissible error. In this case, instead of computing the exact interpolation polynomial p as it is defined above, interpolate computes and returns an approximate interpolation polynomial p~ that does not differ from p more than the amount D with supremum norm over the interval I that is given by the intervall of all abscissa points in X. Formally, let I be the interval hull of all abscissa points in X; typically, the lower bound of I is the minimum value amongst all X, the upper bound is the maximum value amongst all X. Let p be the unique polynomial such that p(xi) = yi for all xi in the list X and yi either the corresponding entry in the list Y or yi = f(xi). Then interpolate returns p~ such that ||p~ - p ||^I <= D. The support for D works for any list of distinct abscissa points X and for any list of ordinate points Y or function f. The algorithm is however optimized for the case when X is a list of floating-point (or rational) constants and a function f. When f is a polynomial of degree n and X contains at least n+1 abscissa points, interpolate may return f as-is, even though D has been provided; typically without performing any rounding on the coefficients of the polynomial f.
• When a permissible-error-argument D is given to interpolate, the command additionnally supports an optional interval parameter I. If this parameter is specified, the interval I to consider the error between p~ and p on is not taken as the interval hull of the abscissa points X but assumed to be the interval I. This parameter is supported to both allow for interpolation over intervals where the first (or last) interpolation point does not lie on the lower (resp. upper) bound of the interval and to speed up computation: determining the minimum and maximum values in X may be expensive.
• interpolate is completely insensible to the tool's working precision prec. It adapts its own working precision entirely automatically. However the user should be aware that precision adaption comes at a price: if the interpolation problem requires high precision to be solved, the execution time of interpolate will be high. This can typically be observed when the list of abscissas X contains two values that are mathematically different but for which high precision is needed to distinguish them.
• In case of error, e.g. when the abscissa points in X are not all distinct, when the function does not evaluate at one of the abscissa points or when the permissible error D or domain interval I are not acceptable, interpolate evaluates to error.

## Example 1:

> p = interpolate([| 1, 2, 3 |], [| 3, 5, 9 |]);
> write("p = ", p, "\n");
p = 3 + _x_ * (-1 + _x_)

## Example 2:

> p = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |]);
> write("p = ", p, "\n");
p = 21 + _x_ * (-100 / 3 + _x_ * (15.5 + _x_ * -13 / 6))
> p = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |], 2^-17);
> write("p = ", p, "\n");
p = 21 + _x_ * (-33.33333301544189453125 + _x_ * (15.5 + _x_ * (-2.16666662693023681640625)))
> p = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |], 2^-17, [0; 2^42]);
> write("p = ", p, "\n");
p = 21 + _x_ * (-33.33333333333333333347789362299806725786766037345 + _x_ * (15.5 + _x_ * (-2.1666666666666666666666666666666666666666666591931)))

## Example 3:

> p = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |]);
> q = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |], 2^-17);
> delta = horner(p - q);
> write("delta = ", delta, "\n");
delta = _x_ * (-1 / 3145728 + _x_^2 * -1 / 25165824)
> Delta = dirtyinfnorm(delta,[1;4]);
> "Delta = ", Delta, " = 2^(", log2(Delta), ")";
Delta = 3.814697265625e-6 = 2^(-18)
> p = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |]);
> q = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |], 2^-2000);
> delta = horner(p - q);
> write("delta = ", delta, "\n");
delta = _x_ * (1 / 2.7555136686582108581587996828264367616535624850129e603 + _x_^2 * 1 / 2.2044109349265686865270397462611494093228499880103e604)
> Delta = dirtyinfnorm(delta,[1;4]);
> "Delta = ", Delta, " = 2^(", log2(Delta), ")";
Delta = 4.3549049081086083377880977473894361479295518713527e-603 = 2^(-2001)

## Example 4:

> p = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |]);
> q = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |], 2^-17, [0;2^10]);
> delta = horner(p - q);
> write("delta = ", delta, "\n");
delta = _x_ * (1 / 1610612736 + _x_^2 * -1 / 1688849860263936)
> Delta = dirtyinfnorm(delta,[0;2^10]);
> "Delta = ", Delta, " = 2^(", log2(Delta), ")";
Delta = 2.4471294368728034316556666925301338011322095712879e-7 = 2^(-21.962406251802890453634347359869541271899536019231)
> p = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |]);
> q = interpolate([| 1, 2, 3, 4 |], [| 1, -1, 2, -3 |], 2^-2000, [0;2^10]);
> delta = horner(p - q);
> write("delta = ", delta, "\n");
delta = _x_ * (-1 / 1.4108229983530039593773054376071356219666239923266e606 + _x_^2 * 1 / 1.47935513632099947970801742654433984193927471937785e612)
> Delta = dirtyinfnorm(delta,[0;2^10]);
> "Delta = ", Delta, " = 2^(", log2(Delta), ")";
Delta = 2.7936728011019194195581558596099001590721546862147e-604 = 2^(-2004.9624062518028904536343473598695412718995360192)

## Example 5:

> p = interpolate([| 1, 2, 3, 4 |], 17 + _x_ * (42 + _x_ * (23 + _x_ * 1664)));
> write("p = ", p, "\n");
p = 17 + _x_ * (42 + _x_ * (23 + _x_ * 1664))

## Example 6:

> p = interpolate([| 1, 2, 3, 4 |], exp(_x_));
> write("p = ", p, "\n");
p = -1 * (-2 * (-3 * ((exp(4) - exp(3) - (exp(3) - exp(2))) / 2 - (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) / 3 + (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) + exp(2) - exp(1)) + exp(1) + _x_ * (-1 * (-2 * ((exp(4) - exp(3) - (exp(3) - exp(2))) / 2 - (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) / 3 + -3 * ((exp(4) - exp(3) - (exp(3) - exp(2))) / 2 - (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) / 3 + (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) + -2 * (-3 * ((exp(4) - exp(3) - (exp(3) - exp(2))) / 2 - (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) / 3 + (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) + exp(2) - exp(1) + _x_ * (-1 * ((exp(4) - exp(3) - (exp(3) - exp(2))) / 2 - (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) / 3 + -2 * ((exp(4) - exp(3) - (exp(3) - exp(2))) / 2 - (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) / 3 + -3 * ((exp(4) - exp(3) - (exp(3) - exp(2))) / 2 - (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) / 3 + (exp(3) - exp(2) - (exp(2) - exp(1))) / 2 + _x_ * ((exp(4) - exp(3) - (exp(3) - exp(2))) / 2 - (exp(3) - exp(2) - (exp(2) - exp(1))) / 2) / 3))
> p = interpolate([| 1, 2, 3, 4 |], exp(_x_), 2^-17);
> write("p = ", p, "\n");
p = -7.71721172332763671875 + _x_ * (17.914661407470703125 + _x_ * (-9.77757251262664794921875 + _x_ * 2.298404276371002197265625))
> p = interpolate([| 1, 2, 3, 4 |], exp(_x_), 2^-17, [0; 2^42]);
> write("p = ", p, "\n");
p = -7.71721172332763671875 + _x_ * (17.9146616149694119287522076078289501310791820287704 + _x_ * (-9.777572455021435040741686047095908782922898736377 + _x_ * 2.2984042886523568836092778582480142756017855238293))

## Example 7:

> p = interpolate([| 1, 2, 3, 4 |], exp(_x_));
> q = interpolate([| 1, 2, 3, 4 |], exp(_x_), 2^-17);
> delta = horner(p - q);
> Delta = dirtyinfnorm(delta,[1;4]);
> "Delta = ", Delta, " = 2^(", log2(Delta), ")";
Delta = 2.64087128985936026120286087840279073703861406872586e-6 = 2^(-18.5305545798246903139105523606201676108779913563466)
> p = interpolate([| 1, 2, 3, 4 |], exp(_x_));
> q = interpolate([| 1, 2, 3, 4 |], exp(_x_), 2^-2000);
> delta = horner(p - q);
> Delta = dirtyinfnorm(delta,[1;4]);
> "Delta = ", Delta, " = 2^(", log2(Delta), ")";
Delta = 1.4590823977486192906712389246433732155807613636698e-603 = 2^(-2002.5775798592063632271660191460324647977792574035)

## Example 8:

> p = interpolate([| 1, 2, 3, 4 |], exp(_x_));
> q = interpolate([| 1, 2, 3, 4 |], exp(_x_), 2^-17, [0;2^10]);
> delta = horner(p - q);
> Delta = dirtyinfnorm(delta,[0;2^10]);
> "Delta = ", Delta, " = 2^(", log2(Delta), ")";
Delta = 4.75462940203480522770747218406076244573282695611e-7 = 2^(-21.0041637691158024834926102190382130962658000640866)
> p = interpolate([| 1, 2, 3, 4 |], exp(_x_));
> q = interpolate([| 1, 2, 3, 4 |], exp(_x_), 2^-2000, [0;2^10]);
> delta = horner(p - q);
> Delta = dirtyinfnorm(delta,[0;2^10]);
> "Delta = ", Delta, " = 2^(", log2(Delta), ")";
Delta = 9.0435150205780172429401306949628326405643331478968e-604 = 2^(-2003.2676856856633076178139050354211927790570412987)

## Example 9:

> f = exp(_x_);
> I = [-1/2;1/2];
> n = 17;
> X = [||]; for i from 0 to n do X = (~(1/2 * (inf(I) + sup(I)) + 1/2 * (sup(I) - inf(I)) * cos(pi * ((2 * i + 1)/(2 * (n + 1)))))) .: X; X = revert(X);
> p = interpolate(X, f, 2^-110);
> q = remez(f, n, I);
> Delta = dirtyinfnorm(p-f,I);
> "||p-f|| = ", Delta, " = 2^(", log2(Delta), ")";
||p-f|| = 4.682282321340942778487515892686081216293095411976e-27 = 2^(-87.464846636624845038438119985399853470805747522283)
> Delta = dirtyinfnorm(q-f,I);
> "||q-f|| = ", Delta, " = 2^(", log2(Delta), ")";
||q-f|| = 4.5615560104462964431174181483603567519672456911095e-27 = 2^(-87.502532530192383004694428978646385337892102159012)
> Delta = dirtyinfnorm(horner(p-q),I);
> "||p-q|| = ", Delta, " = 2^(", log2(Delta), ")";
||p-q|| = 1.20729158377808689215008324924110505247080203008044e-28 = 2^(-92.742212500498778529126446648943330421384492990147)

## Example 10:

> interpolate([| 1, 1 + 2^-1000 |], 17 + 42 * x);
17 + x * 42
> interpolate([| 1, 1 + 2^-10000 |], 17 + 42 * x);
17 + x * 42
> interpolate([| 1, 1 |], 17 + 42 * x);
error