1########################################################################
2##
3## Copyright (C) 1993-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25##
26## Original version by Paul Kienzle distributed as free software in the
27## public domain.
28
29## -*- texinfo -*-
30## @deftypefn {} {} hadamard (@var{n})
31## Construct a Hadamard matrix (@nospell{Hn}) of size @var{n}-by-@var{n}.
32##
33## The size @var{n} must be of the form @math{2^k * p} in which p is one of
34## 1, 12, 20 or 28.  The returned matrix is normalized, meaning
35## @w{@code{Hn(:,1) == 1}} and @w{@code{Hn(1,:) == 1}}.
36##
37## Some of the properties of Hadamard matrices are:
38##
39## @itemize @bullet
40## @item
41## @code{kron (Hm, Hn)} is a Hadamard matrix of size @var{m}-by-@var{n}.
42##
43## @item
44## @code{Hn * Hn' = @var{n} * eye (@var{n})}.
45##
46## @item
47## The rows of @nospell{Hn} are orthogonal.
48##
49## @item
50## @code{det (@var{A}) <= abs (det (Hn))} for all @var{A} with
51## @w{@code{abs (@var{A}(i, j)) <= 1}}.
52##
53## @item
54## Multiplying any row or column by -1 and the matrix will remain a Hadamard
55## matrix.
56## @end itemize
57## @seealso{compan, hankel, toeplitz}
58## @end deftypefn
59
60## Reference [1] contains a list of Hadamard matrices up to n=256.
61## See code for h28 in hadamard.m for an example of how to extend
62## this function for additional p.
63##
64## Reference:
65## [1] A Library of Hadamard Matrices, N. J. A. Sloane
66##     http://www.research.att.com/~njas/hadamard/
67
68function h = hadamard (n)
69
70  if (nargin != 1)
71    print_usage ();
72  endif
73
74  ## Find k if n = 2^k*p.
75  k = 0;
76  while (n > 1 && fix (n/2) == n/2)
77    k += 1;
78    n /= 2;
79  endwhile
80
81  ## Find base hadamard.
82  ## Except for n=2^k, need a multiple of 4.
83  if (n != 1)
84    k -= 2;
85  endif
86
87  ## Trigger error if not a multiple of 4.
88  if (k < 0)
89    n =- 1;
90  endif
91
92  switch (n)
93    case 1
94      h = 1;
95    case 3
96      h = h12 ();
97    case 5
98      h = h20 ();
99    case 7
100      h = h28 ();
101    otherwise
102      error ("hadamard: N must be 2^k*p, for p = 1, 12, 20 or 28");
103  endswitch
104
105  ## Build H(2^k*n) from kron(H(2^k),H(n)).
106  h2 = [1,1;1,-1];
107  while (true)
108    if (fix (k/2) != k/2)
109      h = kron (h2, h);
110    endif
111    k = fix (k/2);
112    if (k == 0)
113      break;
114    endif
115    h2 = kron (h2, h2);
116  endwhile
117
118endfunction
119
120function h = h12 ()
121  tu = [-1,+1,-1,+1,+1,+1,-1,-1,-1,+1,-1];
122  tl = [-1,-1,+1,-1,-1,-1,+1,+1,+1,-1,+1];
123  ## Note: assert (tu(2:end), tl(end:-1:2)).
124  h = ones (12);
125  h(2:end,2:end) = toeplitz (tu, tl);
126endfunction
127
128function h = h20 ()
129  tu = [+1,-1,-1,+1,+1,+1,+1,-1,+1,-1,+1,-1,-1,-1,-1,+1,+1,-1,-1];
130  tl = [+1,-1,-1,+1,+1,-1,-1,-1,-1,+1,-1,+1,-1,+1,+1,+1,+1,-1,-1];
131  ## Note: assert (tu(2:end), tl(end:-1:2)).
132  h = ones (20);
133  h(2:end,2:end) = fliplr (toeplitz (tu, tl));
134endfunction
135
136function h = h28 ()
137  ## Williamson matrix construction from
138  ## http://www.research.att.com/~njas/hadamard/had.28.will.txt
139  ## Normalized so that each row and column starts with +1
140  h = [1 1  1  1  1  1  1  1  1 1  1  1  1 1 1 1 1 1  1 1 1 1 1  1 1  1 1  1
141       1 1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 1 1 1 1 -1 1 1 1 1 -1 1 -1 1 -1
142       1 -1 1 -1 -1 -1 -1 1 -1 1 1 -1 -1 1 -1 -1 -1 -1 1 1 -1 1 -1 1 1 1 1 1
143       1 -1 -1 1 -1 -1 -1 1 1 1 1 1 -1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 1 -1 1
144       1 -1 -1 -1 1 -1 -1 1 1 -1 1 1 1 1 1 -1 -1 -1 1 1 1 -1 1 -1 1 -1 -1 -1
145       1 -1 -1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 1 1 1 1 -1
146       1 -1 -1 -1 -1 -1 1 -1 1 -1 -1 -1 1 -1 1 1 1 -1 1 1 1 1 -1 1 -1 1 -1 1
147       1 -1 1 1 1 1 -1 -1 1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1 1 -1 -1 1 1 1 1 -1
148       1 -1 -1 1 1 1 1 1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 1 -1 -1 1 -1 1
149       1 1 1 1 -1 -1 -1 -1 -1 -1 1 1 1 -1 1 -1 -1 -1 -1 -1 1 1 1 -1 -1 1 1 1
150       1 -1 1 1 1 -1 -1 -1 -1 1 -1 1 1 -1 -1 -1 1 1 -1 1 1 -1 -1 1 1 -1 -1 1
151       1 -1 -1 1 1 1 -1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 1 1 -1 -1 1 -1 -1 1 1 -1
152       1 -1 -1 -1 1 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1
153       1 1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 1 -1 -1 1
154       1 1 -1 1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 -1 -1 1 -1
155       1 1 -1 1 -1 -1 1 -1 1 1 1 -1 -1 -1 -1 -1 -1 1 1 -1 1 -1 1 1 1 1 -1 -1
156       1 1 -1 1 -1 1 1 1 1 1 -1 -1 1 -1 1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 1 1
157       1 1 -1 1 -1 1 -1 1 -1 1 -1 1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 -1 -1
158       1 -1 1 -1 1 -1 1 1 1 1 1 -1 1 -1 -1 1 -1 1 -1 -1 1 1 -1 -1 -1 -1 1 -1
159       1 1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 1 -1 -1
160       1 1 -1 -1 1 -1 1 -1 -1 -1 -1 1 1 1 -1 1 -1 1 1 -1 -1 -1 -1 -1 1 1 1 1
161       1 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1
162       1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 -1 1 -1 -1 1 -1 -1 1 1
163       1 -1 1 -1 -1 1 1 -1 1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 -1 -1 -1 1 1
164       1 1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 -1 1 1 1 -1 1 -1 -1 -1 1 -1 -1
165       1 -1 1 1 -1 1 1 -1 -1 -1 1 -1 1 1 1 -1 1 1 1 -1 -1 1 -1 -1 1 -1 -1 -1
166       1 1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 1 -1 1 -1 1 -1 1 -1 -1 1 1 -1 -1 -1 1
167       1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 1 1 -1 -1 1 1 -1 -1 1 -1];
168endfunction
169
170
171%!assert (hadamard (1), 1)
172%!assert (hadamard (2), [1,1;1,-1])
173%!test
174%! for n = [1,2,4,8,12,24,48,20,28,2^9]
175%!   h = hadamard (n);
176%!   assert (norm (h*h' - n*eye (n)), 0);
177%! endfor
178
179%!error hadamard ()
180%!error hadamard (1,2)
181%!error <N must be 2\^k\*p> hadamard (5)
182