1## Copyright 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz 2## Copyright 1990-2000 Institut für Angewandte Mathematik, 3## Universität Karlsruhe, Germany 4## Copyright 2000-2014 Wissenschaftliches Rechnen/Softwaretechnologie, 5## Universität Wuppertal, Germany 6## Copyright 2015-2016 Oliver Heimlich 7## 8## This program is derived from RPolyEval in CXSC, C++ library for eXtended 9## Scientific Computing (V 2.5.4), which is distributed under the terms of 10## LGPLv2+. Migration to Octave code has been performed by Oliver Heimlich. 11## 12## This program is free software; you can redistribute it and/or modify 13## it under the terms of the GNU General Public License as published by 14## the Free Software Foundation; either version 3 of the License, or 15## (at your option) any later version. 16## 17## This program is distributed in the hope that it will be useful, 18## but WITHOUT ANY WARRANTY; without even the implied warranty of 19## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20## GNU General Public License for more details. 21## 22## You should have received a copy of the GNU General Public License 23## along with this program; if not, see <http://www.gnu.org/licenses/>. 24 25## -*- texinfo -*- 26## @documentencoding UTF-8 27## @defmethod {@@infsup} polyval (@var{P}, @var{X}) 28## 29## Evaluate polynomial @var{P} with argument @var{X}. 30## 31## Horner's scheme is used to evaluate a first approximation. The result is 32## improved with iterative refinement. 33## 34## Accuracy: The result is a tight enclosure for polynomials of degree 1 or 35## less. For polynomials of higher degree the result is a valid enclosure. 36## For @var{X} being no singleton interval, the algorithm suffers from the 37## dependency problem. 38## 39## @example 40## @group 41## output_precision (16, 'local') 42## polyval (infsup ([3 4 2 1]), 42) # 3x^3 + 4x^2 + 2x^1 + 1 | x = 42 43## @result{} [229405] 44## polyval (infsup ([3 4 2 1]), "42?") # ... | x = 41.5 .. 42.5 45## @result{} [221393.125, 237607.875] 46## @end group 47## @end example 48## @seealso{@@infsup/fzero} 49## @end defmethod 50 51## Author: Oliver Heimlich 52## Keywords: interval 53## Created: 2015-05-29 54 55function result = polyval (p, x) 56 57 if (nargin ~= 2) 58 print_usage (); 59 return 60 endif 61 62 if (not (isa (x, "infsup"))) 63 x = infsup (x); 64 endif 65 if (not (isa (p, "infsup"))) 66 p = infsup (p); 67 endif 68 69 if (not (isscalar (x))) 70 error ('point of evaluation X must be a scalar') 71 endif 72 73 if (not (isvector (p))) 74 error ('polynomial P must be a vector of coefficients') 75 endif 76 77 if (isempty (x)) 78 result = x; 79 endif; 80 81 n = numel (p); 82 switch (n) 83 case 0 # empty sum 84 result = infsup (0); 85 return 86 case 1 # p .* x.^0 87 result = p; 88 return 89 case 2 # p(1) .* x.^1 + p(2) .* x.^0 90 result = fma (x, subsref (p, substruct ("()", {1})), ... 91 subsref (p, substruct ("()", {2}))); 92 return 93 endswitch 94 95 if (x == 0) 96 result = subsref (p, substruct ("()", {n})); 97 return 98 elseif (x == 1) 99 result = sum (p); 100 return 101 elseif (x == -1) 102 result = dot (p, (-1) .^ ((n : -1 : 1) - 1)); 103 return 104 endif 105 106 kMax = 20; 107 108 idxNext.type = '()'; 109 idxLast.type = '()'; 110 111 ## Compute first approximation using Horner's scheme 112 yy = infsup (zeros (n, 1)); 113 idxNext.subs = {1}; 114 result = subsref (p, idxNext); 115 yy = subsasgn (yy, idxNext, result); 116 for i = 2 : n 117 idxNext.subs = {i}; 118 result = fma (result, x, subsref (p, idxNext)); 119 yy = subsasgn (yy, idxNext, result); 120 endfor 121 122 y = zeros (n, kMax); 123 for k = 1 : kMax 124 lastresult = result; 125 126 if (isempty (result)) 127 break 128 endif 129 130 ## Iterative refinement 131 ## Store middle of residual as the next correction of y 132 y(:, k) = mid (yy); 133 134 ## Computation of the residual [r] and 135 ## evaluation of the interval system A*[y] = [r] 136 yy = infsup (zeros (n, 1)); 137 for i = 2 : n 138 idxNext.subs = {i}; 139 idxLast.subs = {i-1}; 140 141 coef = subsref (p, idxNext); 142 yy = subsasgn (yy, idxNext, dot ([subsref(yy, idxLast), ... 143 coef, y(i, 1 : k), ... 144 y(i - 1, 1 : k)], ... 145 [x, ... 146 1, ... 147 -ones(1, k), x.*ones(1, k)])); 148 endfor 149 150 ## Determination of a new enclosure of p (x) 151 idxLast.subs = {n}; 152 result = intersect (result, sum ([subsref(yy, idxLast), y(n, 1 : k)])); 153 if (eq (result, lastresult)) 154 ## No improvement 155 break 156 endif 157 if (mpfr_function_d ('plus', +inf, inf (result), pow2 (-1074)) >= ... 158 sup (result)) 159 ## 1 ULP accuracy reached 160 break 161 endif 162 endfor 163 164endfunction 165 166%!assert (polyval (infsup (42), 0) == 42); 167%!assert (polyval (infsup ([42 42]), 0) == 42); 168%!assert (polyval (infsup ([42 42]), 1) == 84); 169%!assert (polyval (infsup ([42 42]), -1) == 0); 170%!assert (polyval (infsup ([-42 42 42]), .5) == -42*0.5^2 + 42*0.5 + 42); 171%!assert (polyval (infsup (vec (pascal (3))), 0.1) == "[0X6.502E9A7231A08P+0, 0X6.502E9A7231A0CP+0]"); 172