Mercurial > octave-nkf

changeset 11110:0be2d25700a7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
polynomial/polyval.m: Use Horner's method.
author Marco Caliari <marco.caliari@univr.it>
date Mon, 18 Oct 2010 20:13:45 +0200
parents 41d18f6342f9
children 84ad75921e35
files scripts/ChangeLog scripts/polynomial/polyval.m
diffstat 2 files changed, 10 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Mon Oct 18 14:15:06 2010 -0400
+++ b/scripts/ChangeLog	Mon Oct 18 20:13:45 2010 +0200
@@ -1,3 +1,7 @@
+2010-10-18  Marco Caliari <marco.caliari@univr.it>
+
+	* polynomial/polyval.m: Use Horner's method.
+
 2010-10-18  John W. Eaton  <jwe@octave.org>
 
 	* plot/__go_draw_axes__.m: Always use gnuplot to display images.
--- a/scripts/polynomial/polyval.m	Mon Oct 18 14:15:06 2010 -0400
+++ b/scripts/polynomial/polyval.m	Mon Oct 18 20:13:45 2010 +0200
@@ -69,10 +69,11 @@
   endif
 
   n = length (p) - 1;
-  k = numel (x);
   x = (x - mu(1)) / mu(2);
-  A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0));
-  y = A * p(:);
+  y = p(1);
+  for i = 2:n+1
+    y = y .* x(:) + p(i);
+  endfor
   y = reshape (y, size (x));
 
   if (nargout == 2)
@@ -82,6 +83,8 @@
     ##   dy = t * sqrt (1 + sumsq (A/s.R, 2)) * s.normr / sqrt (s.df)
     ## If my inference is correct, then t must equal 1 for polyval.
     ## This is because finv (0.5, n, n) = 1.0 for any n.
+    k = numel (x);
+    A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0));
     dy = sqrt (1 + sumsq (A/s.R, 2)) * s.normr / sqrt (s.df);
     dy = reshape (dy, size (x));
   endif