# HG changeset patch
# User Marco Caliari <marco.caliari@univr.it>
# Date 1287425625 -7200
# Node ID 0be2d25700a7f6c04d4619bed5290abac772932b
# Parent  41d18f6342f9f88147842137ae4dd22e0e75d8fa
polynomial/polyval.m: Use Horner's method.

diff -r 41d18f6342f9 -r 0be2d25700a7 scripts/ChangeLog
--- 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.
diff -r 41d18f6342f9 -r 0be2d25700a7 scripts/polynomial/polyval.m
--- 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