octave_packages/m/special-matrix/vander.m

   1 ## Copyright (C) 1993-2012 John W. Eaton
   2 ## Copyright (C) 2009 VZLU Prague
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {} vander (@var{c})
  22 ## @deftypefnx {Function File} {} vander (@var{c}, @var{n})
  23 ## Return the Vandermonde matrix whose next to last column is @var{c}.
  24 ## If @var{n} is specified, it determines the number of columns;
  25 ## otherwise, @var{n} is taken to be equal to the length of @var{c}.
  26 ##
  27 ## A Vandermonde matrix has the form:
  28 ## @tex
  29 ## $$
  30 ## \left[\matrix{c_1^{n-1}  & \cdots & c_1^2  & c_1    & 1      \cr
  31 ##               c_2^{n-1}  & \cdots & c_2^2  & c_2    & 1      \cr
  32 ##               \vdots     & \ddots & \vdots & \vdots & \vdots \cr
  33 ##               c_n^{n-1}  & \cdots & c_n^2  & c_n    & 1      }\right]
  34 ## $$
  35 ## @end tex
  36 ## @ifnottex
  37 ##
  38 ## @example
  39 ## @group
  40 ## c(1)^(n-1) @dots{} c(1)^2  c(1)  1
  41 ## c(2)^(n-1) @dots{} c(2)^2  c(2)  1
  42 ##     .     .      .      .    .
  43 ##     .       .    .      .    .
  44 ##     .         .  .      .    .
  45 ## c(n)^(n-1) @dots{} c(n)^2  c(n)  1
  46 ## @end group
  47 ## @end example
  48 ##
  49 ## @end ifnottex
  50 ## @seealso{polyfit}
  51 ## @end deftypefn
  52
  53 ## Author: jwe
  54
  55 function retval = vander (c, n)
  56
  57   if (nargin == 1)
  58     n = length (c);
  59   elseif (nargin != 2)
  60     print_usage ();
  61   endif
  62
  63   if (! isvector (c))
  64     error ("vander: polynomial C must be a vector");
  65   endif
  66
  67   ## avoiding many ^s appears to be faster for n >= 100.
  68   retval = zeros (length (c), n, class (c));
  69   d = 1;
  70   c = c(:);
  71   for i = n:-1:1
  72     retval(:,i) = d;
  73     d .*= c;
  74   endfor
  75
  76 endfunction
  77
  78
  79 %!test
  80 %! c = [0,1,2,3];
  81 %! expect = [0,0,0,1; 1,1,1,1; 8,4,2,1; 27,9,3,1];
  82 %! assert(vander (c), expect);
  83
  84 %!assert (vander (1), 1)
  85 %!assert (vander ([1, 2, 3]), vander ([1; 2; 3]))
  86 %!assert (vander ([1, 2, 3]), [1, 1, 1; 4, 2, 1; 9, 3, 1])
  87 %!assert (vander ([1, 2, 3]*i), [-1, i, 1; -4, 2i, 1; -9, 3i, 1])
  88
  89 %!assert(vander (2, 3), [4, 2, 1])
  90 %!assert(vander ([2, 3], 3), [4, 2, 1; 9, 3, 1])
  91
  92 %!error vander ();
  93 %!error vander (1, 2, 3);
  94 %!error <polynomial C must be a vector> vander ([1, 2; 3, 4]);
  95