Mercurial > hg > octave-nkf

comparison scripts/polynomial/polyval.m @ 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 a8ce6bdecce5
children fd0a3ac60b0e
comparison
equal deleted inserted replaced
11109:41d18f6342f9 11110:0be2d25700a7
67 y = p; 67 y = p;
68 return; 68 return;
69 endif 69 endif
70 70
71 n = length (p) - 1; 71 n = length (p) - 1;
72 k = numel (x);
73 x = (x - mu(1)) / mu(2); 72 x = (x - mu(1)) / mu(2);
74 A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0)); 73 y = p(1);
75 y = A * p(:); 74 for i = 2:n+1
75 y = y .* x(:) + p(i);
76 endfor
76 y = reshape (y, size (x)); 77 y = reshape (y, size (x));
77 78
78 if (nargout == 2) 79 if (nargout == 2)
79 ## Note: the F-Distribution is generally considered to be single-sided. 80 ## Note: the F-Distribution is generally considered to be single-sided.
80 ## http://www.itl.nist.gov/div898/handbook/eda/section3/eda3673.htm 81 ## http://www.itl.nist.gov/div898/handbook/eda/section3/eda3673.htm
81 ## t = finv (1-alpha, s.df, s.df); 82 ## t = finv (1-alpha, s.df, s.df);
82 ## dy = t * sqrt (1 + sumsq (A/s.R, 2)) * s.normr / sqrt (s.df) 83 ## dy = t * sqrt (1 + sumsq (A/s.R, 2)) * s.normr / sqrt (s.df)
83 ## If my inference is correct, then t must equal 1 for polyval. 84 ## If my inference is correct, then t must equal 1 for polyval.
84 ## This is because finv (0.5, n, n) = 1.0 for any n. 85 ## This is because finv (0.5, n, n) = 1.0 for any n.
86 k = numel (x);
87 A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0));
85 dy = sqrt (1 + sumsq (A/s.R, 2)) * s.normr / sqrt (s.df); 88 dy = sqrt (1 + sumsq (A/s.R, 2)) * s.normr / sqrt (s.df);
86 dy = reshape (dy, size (x)); 89 dy = reshape (dy, size (x));
87 endif 90 endif
88 91
89 endfunction 92 endfunction