# UPB

 Other toolboxes required UPB Generates an unextendible product basis none IsUPBMinUPBSizeUPBSepDistinguishable Unextendible product bases

UPB is a function that generates an unextendible product basis (UPB). The user may either request a specific UPB from the literature such as 'Tiles' or 'Pyramid'[1], or they may request a minimal UPB of specified dimensions.

## Syntax

• U = UPB(NAME)
• U = UPB(NAME,OPT_PAR)
• [U,V,W,...] = UPB(NAME,OPT_PAR)
• U = UPB(DIM)
• U = UPB(DIM,VERBOSE)
• [U,V,W,...] = UPB(DIM,VERBOSE)

## Argument descriptions

### Input arguments

Important: Do not specify both NAME and DIM: just one or the other!

• NAME: The name of a UPB that is found in the literature. Accepted values are:
• 'GenShifts': A minimal UPB in $(\mathbb{C}^2)^{\otimes p}$ (only valid when p ≥ 3 is odd) constructed in [2]. Note that OPT_PAR must be the number of parties (i.e., the integer p) in this case.
• 'Min4x4': A minimal UPB in $\mathbb{C}^4 \otimes \mathbb{C}^4$ constructed in [3].
• 'Pyramid': A minimal UPB in $\mathbb{C}^3 \otimes \mathbb{C}^3$ constructed in [1].
• 'QuadRes': A minimal UPB in $\mathbb{C}^d \otimes \mathbb{C}^d$ (only valid when 2d-1 is prime and d is odd) constructed in [2]. Note that you must set OPT_PAR equal to d (the local dimension) in this case.
• 'Shifts': A minimal UPB in $\mathbb{C}^2 \otimes \mathbb{C}^2 \otimes \mathbb{C}^2$ introduced in [1].
• 'SixParam': The six-parameter UPB in $\mathbb{C}^3 \otimes \mathbb{C}^3$ introduced in Section IV.A of [2]. Note that OPT_PAR must be a vector containing the six parameters in this case.
• 'Tiles': A minimal UPB in $\mathbb{C}^3 \otimes \mathbb{C}^3$ constructed in [1].
• DIM: A vector containing the local dimensions of the desired UPB. In all cases, the smallest known UPB of the desired dimensionality is returned. If no unextendible product basis is known for the specified dimensions, an error is produced.
• VERBOSE (optional, default 1): A flag (either 1 or 0) indicating whether or not a reference to a journal article that contains the UPB (or a description of how to construct the UPB) returned by this script will be displayed.

### Output arguments

If only one output argument is specified (e.g., U = UPB(DIM)) then U is a matrix whose columns are the product states in the desired UPB.

If multiple output arguments are specified (e.g., [U,V,W,...] = UPB(DIM)) then the unextendible product basis is obtained by tensoring the columns of U, V, W, ... together. That is, U, V, W, ... are the local vectors in the unextendible product basis.

## Examples

### Generating the "Shifts" UPB

The following code returns the "Shifts" UPB [1], which is a UPB of 4 states on $\mathbb{C}^2 \otimes \mathbb{C}^2 \otimes \mathbb{C}^2$:

>> v = UPB('Shifts')

v =

1.0000   -0.0000   -0.0000   -0.0000
0   -0.0000    0.0000   -0.5000
0    0.0000   -0.5000   -0.0000
0    0.0000    0.5000   -0.5000
0   -0.5000   -0.0000    0.0000
0   -0.5000    0.0000    0.5000
0    0.5000   -0.5000    0.0000
0    0.5000    0.5000    0.5000

Alternatively, we can request that the local vectors on each copy of $\mathbb{C}^2$ are returned, rather than the total product vectors on $\mathbb{C}^2 \otimes \mathbb{C}^2 \otimes \mathbb{C}^2$:

>> [u,v,w] = UPB('Shifts')

u =

1.0000    0.0000    0.7071   -0.7071
0    1.0000    0.7071    0.7071

v =

1.0000   -0.7071    0.0000    0.7071
0    0.7071    1.0000    0.7071

w =

1.0000    0.7071   -0.7071    0.0000
0    0.7071    0.7071    1.0000

### Generating bound entangled states

As noted in [1], if $$\big\{|v_i\rangle\big\}$$ is an unextendible product basis, then $I - \sum_i |v_i\rangle\langle v_i|$ is (up to scaling) a bound entangled state. The following code illustrates this fact in $\mathbb{C}^3 \otimes \mathbb{C}^5$ by first constructing a UPB in this space, then constructing the corresponding state, and then verifying that this state is bound entangled.

>> v = UPB([3,5]);
Generated a minimal 7-state UPB from:
N. Alon and L. Lovasz. Unextendible product bases. J. Combinatorial Theory, Ser. A, 95:169-179, 2001.

>> rho = eye(3*5);
>> for j = 1:7
rho = rho - v(:,j)*v(:,j)';
end
rho = rho/7; % we are now done constructing the bound entangled state

>> IsSeparable(rho,[3,5]) % show that the state is indeed entangled
Determined to be entangled by not having a 2-copy symmetric extension. Reference:
A. C. Doherty, P. A. Parrilo, and F. M. Spedalieri. A complete family of separability criteria.
Phys. Rev. A, 69:022308, 2004.

ans =

0
>> IsPPT(rho,2,[3,5]) % verify that this state has positive partial transpose and is thus bound entangled

ans =

1

### Generating multipartite UPBs

Many multipartite minimal UPBs can be constructed via this script. For example, the following code generates a minimal UPB (of 6 states) in $$\mathbb{C}^2 \otimes \mathbb{C}^2 \otimes \mathbb{C}^3$$:

>> [u,v,w] = UPB([2,2,3])
Generated a minimal 6-state UPB from:
K. Feng. Unextendible product bases and 1-factorization of complete graphs.
Discrete Applied Mathematics, 154:942–949, 2006.

u =

1.0000         0    1.0000    0.7071    0.7071    0.7071
0    1.0000         0    0.7071   -0.7071    0.7071

v =

1.0000    0.7071         0         0    0.7071    1.0000
0    0.7071    1.0000    1.0000   -0.7071         0

w =

1.0000    0.7071    0.5345    0.5774         0         0
0         0   -0.8018    0.5774    0.3162    1.0000
0    0.7071   -0.2673   -0.5774   -0.9487         0

However, the minimum size of UPBs is still unknown in many multipartite cases – an error is returned in these cases:

>> [u,v,w,x] = UPB([2,2,3,7])
??? Error using ==> UPB at 132
No minimal UPB is currently known in the specified dimensions.

## Source code

Click on "expand" to the right to view the MATLAB source code for this function.

1. %%  UPB  Generates an unextendible product basis
2. %   This function may be called in several different ways:
3. %
4. %   U = UPB(NAME) is a matrix containing as its columns the vectors in the
5. %   unextendible product basis specified by the string NAME. NAME must be
6. %   one of: 'GenShifts', 'Min4x4', 'Pyramid', 'QuadRes', 'Shifts',
7. %   'SixParam', or 'Tiles'. See the online documentation for descriptions
8. %   of these different UPBs.
9. %
10. %   [U,V,W,...] = UPB(NAME) is the same as above, except in this case, the
11. %   unextendible product basis is obtained by tensoring the columns of U,
12. %   V, W, ... together. That is, U, V, W, ... are the local vectors in the
13. %   unextendible product basis.
14. %
15. %   U = UPB(DIM) and [U,V,W,...] = UPB(DIM) are as above, except DIM is a
16. %   vector containing the local dimensions of the requested UPB rather than
17. %   the name of the UPB.
18. %
19. %   U = UPB(DIM,VERBOSE) and [U,V,W,...] = UPB(DIM,VERBOSE) are as above,
20. %   where VERBOSE is a flag (either 1 or 0, default 1) that indicates that
21. %   a reference for the returned UPB will or will not be displayed.
22. %
23. %   URL: http://www.qetlab.com/UPB
24. 
25. %   requires: FourierMatrix.m, IsTotallyNonsingular.m, MinUPBSize.m,
26. %             normalize_cols.m, one_factorization.m, opt_args.m,
27. %             opt_disp.m, perm_inv.m, PermuteSystems.m, sporth.m, Swap.m
28. %
29. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
30. %   package: QETLAB
31. %   last updated: November 12, 2014
32. %   URL: http://www.qetlab.com/UPB
33. 
34. function [u,varargout] = UPB(name,varargin)
35. 
36.     show_name = false;
37.     given_dims = false;
38.     revp = -1; % by default, don't permute systems around after we're done constructing the UPB
39.     if(isnumeric(name)) % user provided dimensions, not a name, so find an appropriate UPB
40.         % set optional argument defaults: verbose=1
41.         [verbose] = opt_args({ 1 },varargin{:});
42. 
43.         given_dims = true;
44.         np = length(name);
45. 
46.         % pre-load various references
47.         refs = {'K. Feng. Unextendible product bases and 1-factorization of complete graphs. Discrete Applied Mathematics, 154:942-949, 2006.', ...
48.                 'D.P. DiVincenzo, T. Mor, P.W. Shor, J.A. Smolin, and B.M. Terhal. Unextendible product bases, uncompletable product bases and bound entanglement. Commun. Math. Phys. 238, 379-410, 2003.', ...
49.                 'C.H. Bennett, D.P. DiVincenzo, T. Mor, P.W. Shor, J.A. Smolin, and B.M. Terhal. Unextendible product bases and bound entanglement. Phys. Rev. Lett. 82, 5385-5388, 1999.', ...
50.                 'T.B. Pedersen. Characteristics of unextendible product bases. Thesis, Aarhus Universitet, Datalogisk Institut, 2002.', ...
51.                 'N. Johnston. The minimum size of qubit unextendible product bases. In Proceedings of the 8th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2013). E-print: arXiv:1302.1604 [quant-ph], 2013.', ...
52.                 'The realm of common sense (if there is only a single party, the only UPBs are full bases of the space).', ...
53.                 'N. Alon and L. Lovasz. Unextendible product bases. J. Combinatorial Theory, Ser. A, 95:169-179, 2001.\nSee also: http://www.njohnston.ca/2013/03/how-to-construct-minimal-upbs/', ...
54.                 'J. Chen and N. Johnston. The minimum size of unextendible product bases in the bipartite case (and some multipartite cases). Comm. Math. Phys., 333(1):351-365, 2015.'};
55. 
56.         [name,revp] = sort(name);
57.         if(np == 1)
58.             upbp{1} = eye(name);
59.             name = '';
60.             ref_ind = 6;
61.         elseif(np == 2 && min(name) <= 2)
62.             % In this case, there are no UPBs smaller than a basis of the
63.             % whole space.
64.             upbp{1} = repmat(eye(name(1)),1,name(2));
65.             upbp{2} = repmat(eye(name(2)),1,name(1));
66.             upbp{revp(1)} = Swap(repmat(eye(name(revp(1))),1,name(revp(2))).',[1,2],[name(revp(2)),name(revp(1))],1).';
67. 
68.             ref_ind = 2;
69.             name = '';
70.         elseif(np == 2 && all(name == [3,3]))
71.             name = 'Tiles';
72.             ref_ind = 3;
73.             show_name = true;
74.         elseif(np == 2 && all(name == [4,4]))
75.             name = 'Feng4x4';
76.             ref_ind = 1;
77.         elseif(np == 2 && all(name == name(1)) && mod(name(1),2) == 1 && isprime(2*name(1)-1))
78.             varargin = {name(1)};
79.             name = 'QuadRes';
80.             ref_ind = 2;
81.             show_name = true;
82.         elseif(np == 3 && all(name == [2,2,2]))
83.             name = 'Shifts';
84.             ref_ind = 3;
85.             show_name = true;
86.         elseif(np == 3 && all(name == [2,2,3]))
87.             name = 'Feng2x2x3';
88.             ref_ind = 1;
89.         elseif(np == 3 && all(name == [2,2,5]))
90.             name = 'Feng2x2x5';
91.             ref_ind = 1;
92.         elseif(np == 4 && all(name == [2,2,2,2]))
93.             name = 'Feng2x2x2x2';
94.             ref_ind = 1;
95.         elseif(np == 4 && all(name == [2,2,2,4]))
96.             name = 'Feng2x2x2x4';
97.             ref_ind = 1;
98.         elseif(np == 5 && all(name == [2,2,2,2,5]))
99.             name = 'Feng2x2x2x2x5';
100.             ref_ind = 1;
101.         elseif(np == 8 && all(name == [2,2,2,2,2,2,2,2]))
102.             name = 'John2^8';
103.             ref_ind = 5;
104.         elseif(mod(np,4) == 0 && all(name == 2*ones(1,np)))
105.             varargin = {np};
106.             name = 'John2^4k';
107.             ref_ind = 5;
108.         elseif(mod(np,2) == 1 && all(name == 2*ones(1,np)))
109.             varargin = {np};
110.             name = 'GenShifts';
111.             ref_ind = 2;
112.             show_name = true;
113.         elseif((mod(sum(name)-np,2) == 1 || sum(mod(name,2)==0) == 0) && sum(mod(name,2)==0) <= 1)
114.             varargin = {name};
115.             name = 'AlonLovasz';
116.             ref_ind = 7;
117.         elseif(name(end)-1 == sum(name(1:end-1)-1) && sum(name-1) >= 3 && mod(sum(name)-np,2) == 0)
118.             varargin = {name};
119.             name = 'CJBip';
120.             ref_ind = 8;
121.         elseif(np == 2 && all(name == [4,6]))
122.             name = 'CJBip46';
123.             ref_ind = 8;
124.         elseif(np == 3 && name(1) == 2 && name(2) == 2 && mod(name(3),4) == 1)
125.             varargin = {name(3)};
126.             name = 'CJ4k1';
127.             ref_ind = 8;
128.         else
129.             try
130.                 min_size = MinUPBSize(name,0);
131.             catch err
132.                 if(strcmpi(err.identifier,'MinUPBSize:MinSizeUnknown'))
133.                     error('UPB:MinSizeUnknown','No minimal UPB is currently known in the specified dimensions.');
134.                 else
135.                     rethrow(err);
136.                 end
137.             end
138.             error('UPB:HardToConstruct',['Minimal UPBs are known to have size ',num2str(min_size),' in this case, but their construction is complicated and not implemented by this script.']);
139.         end
140.     end
141. 
142.     if(strcmpi(name,'Shifts')) % GenShifts reduces to Shifts when there are 3 parties
143.         name = 'GenShifts';
144.         varargin = {3};
145.     end
146. 
147.     % the "Pyramid" UPB
148.     if(strcmpi(name,'Pyramid'))
149.         h = sqrt(1 + sqrt(5))/2;
150.         for j = 4:-1:0 % pre-allocate
151.             upbp{1}(:,j+1) = [cos(2*pi*j/5);sin(2*pi*j/5);h];
152.         end
153.         upbp{1} = 2*upbp{1}/sqrt(5+sqrt(5));
154.         upbp{2} = upbp{1}(:,[1,3,5,2,4]);
155. 
156.     % the "Tiles" UPB
157.     elseif(strcmpi(name,'Tiles'))
158.         upbp{1}(:,5) = ones(3,1)/sqrt(3); % pre-allocate
159.         upbp{1}(:,1) = [1;0;0];
160.         upbp{1}(:,2) = [1;-1;0]/sqrt(2);
161.         upbp{1}(:,3) = [0;0;1];
162.         upbp{1}(:,4) = [0;1;-1]/sqrt(2);
163.         upbp{2}(:,5) = ones(3,1)/sqrt(3);
164.         upbp{2}(:,1) = [1;-1;0]/sqrt(2);
165.         upbp{2}(:,2) = [0;0;1];
166.         upbp{2}(:,3) = [0;1;-1]/sqrt(2);
167.         upbp{2}(:,4) = [1;0;0];
168. 
169.     % the "Min4x4" UPB
170.     elseif(strcmpi(name,'Min4x4'))
171.         upbp{1}(:,8) = [0;0;1;0]; % pre-allocate
172.         upbp{1}(:,1) = [1;-3;1;1]/sqrt(12);
173.         upbp{1}(:,2) = [1;0;0;0];
174.         upbp{1}(:,3) = [0;1;2;1]/sqrt(6);
175.         upbp{1}(:,4) = [1;0;0;-1]/sqrt(2);
176.         upbp{1}(:,5) = [0;1;0;0];
177.         upbp{1}(:,6) = [3;1;-1;1]/sqrt(12);
178.         upbp{1}(:,7) = [0;1;1;0]/sqrt(2);
179.         upbp{2}(:,8) = [-1;1+sqrt(2);0;1]/sqrt(5+2*sqrt(2)); % pre-allocate
180.         upbp{2}(:,1) = [0;1;-3-sqrt(2);-1-sqrt(2)]/sqrt(15+8*sqrt(2));
181.         upbp{2}(:,2) = [1;0;0;0];
182.         upbp{2}(:,3) = [1;0;sqrt(2)-1;1]/sqrt(5-2*sqrt(2));
183.         upbp{2}(:,4) = [0;1;0;0];
184.         upbp{2}(:,5) = [-1;1+sqrt(2);0;1]/sqrt(5+2*sqrt(2));
185.         upbp{2}(:,6) = [0;0;1;0];
186.         upbp{2}(:,7) = [1;1;1;-sqrt(2)]/sqrt(5);
187. 
188.     % the "QuadRes" UPB
189.     elseif(strcmpi(name,'QuadRes'))
190.         if(isempty(varargin))
191.             error('UPB:InvalidArguments','When using NAME=QuadRes, you must specify a second input argument that gives the dimension of the desires UPB.')
192.         elseif(~isprime(2*varargin{1}-1))
193.             error('UPB:InvalidArguments','When using NAME=QuadRes, the second input argument DIM must be such that 2*DIM-1 is prime.')
194.         elseif(mod(varargin{1},2) == 0)
195.             error('UPB:InvalidArguments','When using NAME=QuadRes, the second input argument DIM must be odd.')
196.         end
197. 
198.         p = 2*varargin{1}-1;
199.         q = quad_residue(p);
200.         s = setdiff(1:p-1,q);
201.         sm = sum(exp(2i*pi*q/p));
202.         N = max(-sm,1+sm);
203. 
204.         % This UPB is obtained by choosing entries from the Fourier matrix
205.         % carefully.
206.         F = FourierMatrix(p);
207.         F(1,:) = sqrt(N)*F(1,:);
208.         upbp{1} = F([1,q+1],:);
209.         upbp{2} = F([1,mod(s(1)*q,p)+1],:);
210. 
211.         % normalize the output
212.         upbp{1} = upbp{1}./repmat(sqrt(sum(abs(upbp{1}).^2,1)),varargin{1},1);
213.         upbp{2} = upbp{2}./repmat(sqrt(sum(abs(upbp{2}).^2,1)),varargin{1},1);
214. 
215.     % the "SixParam" UPB (i.e., the one from Section IV.A of DMSST03)
216.     elseif(strcmpi(name,'SixParam'))
217.         if(isempty(varargin) || length(varargin{1}) ~= 6)
218.             error('UPB:InvalidArguments','When using NAME=SixParam, you must specify a vector containing the six parameters [gammaA,thetaA,phiA,gammaB,thetaB,phiB].')
219.         elseif(any(abs(sin(varargin{1}([1,2,4,5]))) < 10*eps) || any(abs(cos(varargin{1}([1,2,4,5]))) < 10*eps))
220.             error('UPB:InvalidArguments','When using NAME=SixParam, none of the gammaA, thetaA, gammaB, and thetaB parameters can be a multiple of pi/2.')
221.         end
222. 
223.         arg_cell = num2cell(varargin{1});
224.         [gammaA,thetaA,phiA,gammaB,thetaB,phiB] = arg_cell{:}; % give more memorable names to parameters
225. 
226.         NA = sqrt(cos(gammaA)^2 + sin(gammaA)^2 * cos(thetaA)^2);
227.         NB = sqrt(cos(gammaB)^2 + sin(gammaB)^2 * cos(thetaB)^2);
228. 
229.         upbp{1}(:,5) = [0;sin(gammaA)*cos(thetaA)*exp(1i*phiA);cos(gammaA)]/NA; % pre-allocate
230.         upbp{1}(:,1) = [1;0;0];
231.         upbp{1}(:,2) = [0;1;0];
232.         upbp{1}(:,3) = [cos(thetaA);0;sin(thetaA)];
233.         upbp{1}(:,4) = [sin(gammaA)*sin(thetaA);cos(gammaA)*exp(1i*phiA);-sin(gammaA)*cos(thetaA)];
234.         upbp{2}(:,5) = [0;sin(gammaB)*cos(thetaB)*exp(1i*phiB);cos(gammaB)]/NB;
235.         upbp{2}(:,1) = [0;1;0];
236.         upbp{2}(:,2) = [sin(gammaB)*sin(thetaB);cos(gammaB)*exp(1i*phiB);-sin(gammaB)*cos(thetaB)];
237.         upbp{2}(:,3) = [1;0;0];
238.         upbp{2}(:,4) = [cos(thetaB);0;sin(thetaB)];
239. 
240.     % the "GenShifts" UPB
241.     elseif(strcmpi(name,'GenShifts'))
242.         if(isempty(varargin))
243.             error('UPB:InvalidArguments','When using NAME=GenShifts, you must specify a second input argument that gives the number of parties in the desired UPB.')
244.         elseif(mod(varargin{1},2) == 0)
245.             error('UPB:InvalidArguments','When using NAME=GenShifts, the second input argument P must be odd.')
246.         end
247. 
248.         k = (varargin{1}+1)/2;
249.         upbp{1} = [cos((0:(1/k):(2-1/k))*pi/2);sin((0:(1/k):(2-1/k))*pi/2)];
250.         upbp{1} = upbp{1}(:,[1,k+1:-1:2,k+2:end]);
251. 
252.         for j = varargin{1}-1:-1:1 % pre-allocate memory
253.             upbp{j+1} = upbp{1}(:,[1,circshift(2:(2*k),[0,j])]);
254.         end
255. 
256.     % the "Feng2x2x2x2" UPB
257.     elseif(strcmpi(name,'Feng2x2x2x2'))
258.         b1 = eye(2);
259.         b2 = [1 1;1 -1]/sqrt(2);
260.         b3 = [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)];
261. 
262.         upbp{1} = [b1(:,1),b1(:,2),b1(:,1),b2(:,1),b2(:,2),b2(:,1)];
263.         upbp{2} = [b1(:,1),b2(:,1),b1(:,2),b1(:,2),b2(:,2),b1(:,1)];
264.         upbp{3} = [b1(:,1),b2(:,1),b3(:,1),b2(:,2),b3(:,2),b1(:,2)];
265.         upbp{4} = [b1(:,1),b2(:,1),b3(:,1),b3(:,2),b1(:,2),b2(:,2)];
266. 
267.     % the "John2^8" UPB
268.     elseif(strcmpi(name,'John2^8'))
269.         b1 = eye(2);
270.         b2 = [1 1;1 -1]/sqrt(2);
271.         b3 = [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)];
272.         b4 = [cos(pi/5) sin(pi/5);-sin(pi/5) cos(pi/5)];
273.         b5 = [cos(pi/7) sin(pi/7);-sin(pi/7) cos(pi/7)];
274. 
275.         upbp{1} = [b1(:,1),b2(:,1),b2(:,2),b1(:,2),b3(:,1),b4(:,1),b4(:,1),b3(:,1),b1(:,2),b3(:,2),b4(:,2)];
276.         upbp{2} = [b1(:,1),b1(:,2),b2(:,1),b3(:,1),b4(:,1),b4(:,1),b4(:,2),b4(:,2),b3(:,1),b2(:,2),b3(:,2)];
277.         upbp{3} = [b1(:,1),b2(:,1),b3(:,1),b3(:,2),b1(:,2),b3(:,2),b1(:,2),b4(:,1),b3(:,1),b2(:,2),b4(:,2)];
278.         upbp{4} = [b1(:,1),b2(:,1),b3(:,1),b3(:,1),b3(:,2),b4(:,1),b3(:,2),b2(:,2),b4(:,1),b4(:,2),b1(:,2)];
279.         upbp{5} = [b1(:,1),b2(:,1),b1(:,2),b3(:,1),b4(:,1),b2(:,2),b5(:,1),b3(:,2),b3(:,1),b5(:,2),b4(:,2)];
280.         upbp{6} = [b1(:,1),b2(:,1),b3(:,1),b2(:,2),b4(:,1),b4(:,2),b5(:,1),b5(:,2),b2(:,2),b1(:,2),b3(:,2)];
281.         upbp{7} = [b1(:,1),b2(:,1),b3(:,1),b4(:,1),b2(:,2),b4(:,2),b2(:,2),b1(:,2),b3(:,2),b5(:,1),b5(:,2)];
282.         upbp{8} = [b1(:,1),b2(:,1),b3(:,1),b4(:,1),b5(:,1),b1(:,2),b5(:,1),b3(:,2),b5(:,2),b4(:,2),b2(:,2)];
283. 
284.     % the "John2^4k" UPB
285.     elseif(strcmpi(name,'John2^4k'))
286.         if(isempty(varargin))
287.             error('UPB:InvalidArguments','When using NAME=John2^4k, you must specify a second input argument that gives the number of parties in the desired UPB.')
288.         elseif(mod(varargin{1},4) ~= 0 || varargin{1} <= 7)
289.             error('UPB:InvalidArguments','When using NAME=John2^4k, the second input argument P must equal 0 (mod 4) and must be at least 8.')
290.         end
291. 
292.         % Construct varargin{1}/2+2 distinct orthonormal bases of C^2.
293.         for j = (varargin{1}/2+2):-1:1 % pre-allocate
294.             jb = 2*pi/(2*j-1);
295.             b(:,:,j) = [cos(jb) sin(jb);-sin(jb) cos(jb)];
296.         end
297. 
298.         % the first 3 parties are special, so we do them separately
299.         upbp{1} = [reshape(repmat(reshape(b(:,1,1:varargin{1}/4+1),2,varargin{1}/4+1),2,1),2,varargin{1}/2+2), reshape(repmat(reshape(b(:,2,1:varargin{1}/4+1),2,varargin{1}/4+1),2,1),2,varargin{1}/2+2)];
300.         upbp{2} = [upbp{1}(:,1:varargin{1}/2+2),circshift(upbp{1}(:,varargin{1}/2+3:end),[0,2])];
301.         upbp{3} = [upbp{1}(:,1:varargin{1}/2+2),circshift(upbp{1}(:,varargin{1}/2+3:end),[0,4])];
302. 
303.         % now do the next 2k-4 parties
304.         if(varargin{1} > 8)
305.             upbp{4} = [reshape(b(:,1,:),2,varargin{1}/2+2), circshift(reshape(b(:,2,:),2,varargin{1}/2+2),[0,6])];
306.             upbp{5} = upbp{4}(:,[1:varargin{1}/2+2,reshape(fliplr(reshape(varargin{1}/2+3:varargin{1}+4,2,varargin{1}/4+1).').',1,varargin{1}/2+2)]);
307. 
308.             for j = 3:(varargin{1}/4-1)
309.                 upbp{2*j} = [upbp{4}(:,1:varargin{1}/2+2),circshift(upbp{4}(:,varargin{1}/2+3:end),[0,2*(j-1)])];
310.                 upbp{2*j+1} = [upbp{5}(:,1:varargin{1}/2+2),circshift(upbp{5}(:,varargin{1}/2+3:end),[0,2*(j-1)])];
311.             end
312.         end
313. 
314.         % Finally, do the last 2k+1 parties, which arise from finding a
315.         % 1-factorization of the complete graph.
316.         fac = one_factorization(varargin{1}/2+2);
317.         resb = reshape(b,2,varargin{1}+4);
318.         for j = 1:(varargin{1}/2+1)
319.             upbp{varargin{1}/2+j-1} = resb(:,[fac(j,:),varargin{1}/2+2+fac(j,:)]);
320.         end
321. 
322.     % the UPB from Theorem 3 of reference [8]
323.     elseif(strcmpi(name,'CJ4k1'))
324.         % Construct (varargin{1}+1)/2+1 distinct orthonormal bases of C^2.
325.         for j = ((varargin{1}+1)/2+1):-1:1 % pre-allocate
326.             jb = 2*pi/(2*j-1);
327.             b(:,:,j) = [cos(jb) sin(jb);-sin(jb) cos(jb)];
328.         end
329. 
330.         % the 2-dimensional parties are the same as in the Feng4m2 UPB
331.         upbp{1} = repmat(reshape(b(:,:,1:((varargin{1}+3)/4)),2,(varargin{1}+1)/2+1),1,2);
332.         u_ind = reshape((varargin{1}+1)/2+2:varargin{1}+3,2,(varargin{1}+3)/4);
333.         upbp{1}(:,(varargin{1}+1)/2+2:varargin{1}+3) = upbp{1}(:,reshape(u_ind([2,1],:),1,(varargin{1}+1)/2+1));
334.         b = reshape(b,2,varargin{1}+3);
335.         upbp{3} = CJ_Lemma6((varargin{1}-1)/4);
336.         upbp{2} = b(:,circshift(1:varargin{1}+3,[0,1]));
337.         upbp{2}(:,[(varargin{1}+3)/2,varargin{1}+3]) = upbp{2}(:,[varargin{1}+3,(varargin{1}+3)/2]);
338. 
339.     % the UPB from Theorem 1 of reference [8]
340.     elseif(strcmpi(name,'CJBip'))
341.         maxD = varargin{1}(end);
342.         U = HollowUnitary(maxD);
343.         upbp{np} = [eye(maxD),U];
344. 
345.         if(maxD >= 20 && ~IsTotallyNonsingular(U,2:maxD-2))
346.             error('UPB:HardToConstruct',['Minimal UPBs are known to have size ',num2str(maxD*2),' in this case, but their construction is complicated and not implemented by this script.']);
347.         end
348. 
349.         prevdims = 0;
350.         for j = 1:np-1
351.             upbp{j} = CJ_Lemma5(maxD,varargin{1}(j)-1,1+prevdims);
352.             prevdims = prevdims + varargin{1}(j)-1;
353.         end
354. 
355.     % another UPB from Theorem 1 of reference [8]
356.     elseif(strcmpi(name,'CJBip46'))
357.         u = 1.64451358502312496885542269243;
358.         upbp{1} = normalize_cols([3 3 3 1 1 -1 -1 2 -2 0;2 1 1 2 3 -2 0 0 -2 -1;2 1 1 1 1 5 -3 -4 2 1;2 2 3 2 2 0 2 1 3 0]);
359.         upbp{2} = eye(6);
360.         upbp{2} = [upbp{2}(:,1:5),normalize_cols([0 0 (u-2)/(1-u) 1 (u-1)/(u-2);u*(u-2)/(2*u-3) 0 0 (3-2*u)/(u-2) 1;1 -1-u 0 0 1/(1+u);1 1 -u 0 0;0 u-1 1 1/(1-u) 0;u 1 1 1 1])];
361. 
362.     % the "Feng2x2x3" UPB
363.     elseif(strcmpi(name,'Feng2x2x3'))
364.         b1 = eye(2);
365.         b2 = [1 1;1 -1]/sqrt(2);
366. 
367.         upbp{3} = normalize_cols([1,1,2,1,0,0;0,0,-3,1,1,1;0,1,-1,-1,-3,0]);
368.         upbp{1} = [b1(:,1),b1(:,2),b1(:,1),b2(:,1),b2(:,2),b2(:,1)];
369.         upbp{2} = [b1(:,1),b2(:,1),b1(:,2),b1(:,2),b2(:,2),b1(:,1)];
370. 
371.     % the "Feng2x2x5" UPB
372.     elseif(strcmpi(name,'Feng2x2x5'))
373.         w = exp(pi*2i/3);
374.         b1 = eye(2);
375.         b2 = [1 1;1 -1]/sqrt(2);
376.         b3 = [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)];
377.         b4 = [cos(pi/5) sin(pi/5);-sin(pi/5) cos(pi/5)];
378. 
379.         upbp{3} = normalize_cols([1,conj(w),w,0,0,0,1,0;0,w,conj(w),1,0,1,0,0;0,conj(w),0,0,0,1,w,1;0,0,w,0,1,1,conj(w),0;0,1,1,0,0,1,1,0]);
380.         upbp{1} = [b2(:,2),b2(:,1),b1(:,2),b1(:,1),b1(:,1),b1(:,2),b2(:,1),b2(:,2)];
381.         upbp{2} = [b4(:,2),b1(:,2),b4(:,1),b1(:,1),b3(:,1),b2(:,1),b3(:,2),b2(:,2)];
382. 
383.     % the "Feng4m2" UPB of Theorem 3.2 from the Feng paper
384.     % I don't know how to compute a 1-factorization of the complement graph
385.     % in this UPB's construction yet. The commented section below shows the
386.     % part of the construction that I *do* know how to implement.
387. %    elseif(strcmpi(name,'Feng4m2'))
388. %        if(isempty(varargin))
389. %            error('UPB:InvalidArguments','When using NAME=Feng4m2, you must specify a second input argument that gives the number of parties in the desired UPB.')
390. %        elseif(mod(varargin{1},4) ~= 2)
391. %            error('UPB:InvalidArguments','When using NAME=Feng4m2, the second input argument P must equal 2 (mod 4).')
392. %        end
393. %
394. %        % Construct varargin{1}/2+1 distinct orthonormal bases of C^2.
395. %        for j = (varargin{1}/2+1):-1:1 % pre-allocate
396. %            jb = 2*pi/(2*j-1);
397. %            b(:,:,j) = [cos(jb) sin(jb);-sin(jb) cos(jb)];
398. %        end
399. %
400. %        upbp{1} = repmat(reshape(b(:,:,1:((varargin{1}+2)/4)),2,varargin{1}/2+1),1,2); % these are the a_j's in the proof of the theorem
401. %        u_ind = reshape(varargin{1}/2+2:varargin{1}+2,2,(varargin{1}+2)/4);% we need to move the columns of upbp{1} around a little bit first
402. %        upbp{1}(:,varargin{1}/2+2:varargin{1}+2) = upbp{1}(:,reshape(u_ind([2,1],:),1,varargin{1}/2+1));
403. %        b = reshape(b,2,varargin{1}+2); % these are the vectors that will be used in all subsystems except the first
404. 
405.     % the "Feng4x4" UPB (Theorem 3.3(4) of Feng paper)
406.     elseif(strcmpi(name,'Feng4x4'))
407.         upbp{1}(:,8) = [1;-1;1;0]/sqrt(3); % pre-allocate
408.         upbp{1}(:,1:4) = eye(4);
409.         upbp{1}(:,5) = [0;1;1;1]/sqrt(3);
410.         upbp{1}(:,6) = [1;0;-1;1]/sqrt(3);
411.         upbp{1}(:,7) = [1;1;0;-1]/sqrt(3);
412.         upbp{2}(:,8) = [0;1;1;1]/sqrt(3);
413.         upbp{2}(:,[1,7,6,4]) = eye(4);
414.         upbp{2}(:,5) = [1;-1;1;0]/sqrt(3);
415.         upbp{2}(:,2) = [1;0;-1;1]/sqrt(3);
416.         upbp{2}(:,3) = [1;1;0;-1]/sqrt(3);
417. 
418.     % the "Feng2x2x2x4" UPB (Theorem 3.3(3) of Feng paper)
419.     elseif(strcmpi(name,'Feng2x2x2x4'))
420.         % need four different bases of C^2
421.         b1 = eye(2);
422.         b2 = [1 1;1 -1]/sqrt(2);
423.         b3 = [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)];
424.         b4 = [cos(pi/5) sin(pi/5);-sin(pi/5) cos(pi/5)];
425. 
426.         upbp{4}(:,8) = [1;-1;1;0]/sqrt(3); % pre-allocate
427.         upbp{4}(:,1:4) = eye(4);
428.         upbp{4}(:,5) = [0;1;1;1]/sqrt(3);
429.         upbp{4}(:,6) = [1;0;-1;1]/sqrt(3);
430.         upbp{4}(:,7) = [1;1;0;-1]/sqrt(3);
431. 
432.         upbp{1} = [b1 b2 b3 b4];
433.         upbp{3} = [b1 b2 b3 b4];
434.         upbp{2} = [b1 b2 b3 b4];
435.         upbp{2} = upbp{2}(:,[1,3,7,5,4,2,6,8]);
436.         upbp{3} = upbp{3}(:,[1,5,7,3,4,8,6,2]);
437.         upbp{1} = upbp{1}(:,[1,3,7,5,8,6,2,4]);
438. 
439.     % the "Feng2x2x2x2x5" UPB (Theorem 3.3(5) of Feng paper)
440.     elseif(strcmpi(name,'Feng2x2x2x2x5'))
441.         w = exp(pi*2i/3);
442. 
443.         upbp{5}(:,10) = [1;conj(w);w;1;0]/2; % pre-allocate
444.         upbp{5}(:,1:5) = eye(5);
445.         upbp{5}(:,6) = [0;1;1;1;1]/2;
446.         upbp{5}(:,7) = [1;0;1;w;conj(w)]/2;
447.         upbp{5}(:,8) = [1;1;0;conj(w);w]/2;
448.         upbp{5}(:,9) = [1;w;conj(w);0;1]/2;
449. 
450.         % need five different bases of C^2
451.         upbp{1} = [eye(2), [1 1;1 -1]/sqrt(2), [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)], [cos(pi/5) sin(pi/5);-sin(pi/5) cos(pi/5)], [cos(pi/7) sin(pi/7);-sin(pi/7) cos(pi/7)]];
452.         upbp{4} = upbp{1};
453.         upbp{3} = upbp{1};
454.         upbp{2} = upbp{1};
455. 
456.         upbp{2} = upbp{2}(:,[1,3,5,7,9,10,2,4,6,8]);
457.         upbp{3} = upbp{3}(:,[1,3,5,7,9,8,10,2,4,6]);
458.         upbp{4} = upbp{4}(:,[1,3,5,7,9,6,8,10,2,4]);
459.         upbp{1} = upbp{1}(:,[1,3,5,7,9,4,6,8,10,2]);
460.     elseif(strcmpi(name,'AlonLovasz'))
461.         dim = varargin{1};
462.         min_size = sum(dim) - np + 1;
463.         os = 0;
464. 
465.         if(mod(dim(1),2) == 0)
466.             upbp{1} = AL_even(dim(1),min_size);
467.         else
468.             upbp{1} = AL_odd(dim(1),min_size,os);
469.             os = os + (dim(1)-1)/2;
470.         end
471.         if(length(upbp{1}) == 1 && upbp{1} == -1)
472.             error('UPB:HardToConstruct',['Minimal UPBs are known to have size ',num2str(min_size),' in this case, but their construction is complicated and not implemented by this script.']);
473.         end
474. 
475.         for j = np:-1:2
476.             if(mod(dim(j),2) == 0)
477.                 upbp{j} = AL_even(dim(j),min_size);
478.             else
479.                 upbp{j} = AL_odd(dim(j),min_size,os);
480.                 os = os + (dim(j)-1)/2;
481.             end
482.             if(length(upbp{j}) == 1 && upbp{j} == -1)
483.                 error('UPB:HardToConstruct',['Minimal UPBs are known to have size ',num2str(min_size),' in this case, but their construction is complicated and not implemented by this script.']);
484.             end
485.         end
486.     end
487. 
488.     if(length(revp) == 1)
489.         revp = 1:length(upbp);
490.     end
491.     revp = perm_inv(revp);
492.     u = upbp{revp(1)};
493.     varargout = upbp(revp(2:end));
494. 
495. 
496.     % If the user just requested one output argument, tensor local vectors
497.     % together.
498.     if(nargout <= 1 && length(varargout) >= 1)
499.         num_vec = size(u,2);
500.         % Do a clever column-by-column kronecker product that avoids having to
501.         % loop over every vector within each party. Note that this method
502.         % relies on not using sparse matrices.
503.         u = reshape(u,1,size(u,1),num_vec);
504.         for k = 1:length(varargout)-1
505.             v = reshape(varargout{k},size(varargout{k},1),1,num_vec);
506.             u = reshape(bsxfun(@times,v,u),1,size(u,2)*size(v,1),num_vec);
507.         end % we split the last iteration outside of the loop to reduce the number of required reshapes by 1
508.         v = reshape(varargout{end},size(varargout{end},1),1,num_vec);
509.         u = reshape(bsxfun(@times,u,v), size(u,2)*size(varargout{end},1),num_vec);
510.     end
511. 
512.     if(given_dims)
513.         if(show_name)
514.             opt_disp(['Generated the ',num2str(length(upbp{1})),'-state ''',name,''' UPB from:\n',refs{ref_ind},'\n'],verbose);
515.         else
516.             opt_disp(['Generated a minimal ',num2str(length(upbp{1})),'-state UPB from:\n',refs{ref_ind},'\n'],verbose);
517.         end
518.     end
519. end
520. 
521. % Computes all quadratic residues mod p (needed for the QuadRes UPB)
522. function q = quad_residue(p)
523.     q = unique(mod((1:floor(p/2)).^2,p));
524. end
525. 
526. % Computes the vectors in an Alon-Lovasz UPB in odd dimensions.
527. function W = AL_odd(r,c,os)
528. 
529.     its = 0;
530.     lin_indep = false;
531. 
532.     while ~lin_indep
533.         W = randn(r,c);
534. 
535.         for j = 2:c
536.             ind = intersect(1:j-1,mod([(j-(r-1)/2-os):(j-1-os),(j+1+os):(j+(r-1)/2+os)]-1,c)+1);
537.             tmp = null(W(:,ind)');
538.             W(:,j) = tmp*randn(size(tmp,2),1);
539.         end
540. 
541.         W = normalize_cols(W);
542.         lin_indep = IsTotallyNonsingular(W,r);
543.         its = its + 1;
544. 
545.         if(its >= 2 && ~lin_indep)
546.             W = -1;
547.             return
548.         end
549.     end
550. end
551. 
552. % Computes the vectors in an Alon-Lovasz UPB in one even dimension.
553. function W = AL_even(r,c)
554. 
555.     its = 0;
556.     lin_indep = false;
557. 
558.     while ~lin_indep
559.         r2 = c - r + 1;
560.         W = randn(r,c);
561. 
562.         for j = 2:c
563.             ind = setdiff(1:c,mod((j-(r2-1)/2:j+(r2-1)/2)-1,c)+1);
564.             tmp = null(W(:,ind)');
565.             W(:,j) = tmp*randn(size(tmp,2),1);
566.         end
567. 
568.         W = normalize_cols(W);
569.         lin_indep = IsTotallyNonsingular(W,r);
570.         its = its + 1;
571. 
572.         if(its >= 2 && ~lin_indep)
573.             W = -1;
574.             return
575.         end
576.     end
577. end
578. 
579. %%  CJ_LEMMA5    Construct a matrix W described by Lemma 5 of [CJ]
580. %   This function has three required argument:
581. %     Q: the positive integer q described by the lemma
582. %     R: the positive integer r described by the lemma
583. %     S: the positive integer s described by the lemma
584. %
585. %   W = CJ_Lemma5(Q,R,S) is a (R+1)-by-2Q matrix with columns of unit
586. %   length that satisfy the orthogonality condition (a) and the
587. %   nonsingularity condition (b) of Lemma 5.
588. %
589. %   This function has one optional argument:
590. %     NICE (default 0): a non-negative integer that specifies how "nice"
591. %                       the entries of W should be
592. %
593. %   W = CJ_Lemma5(Q,R,S,NICE) is the same as before if NICE = 0. If NICE is
594. %   positive integer, then W will be a symbolic matrix (rather than a
595. %   numeric matrix) whose entries are integers. In this case, the columns
596. %   of W are no longer scaled to have unit length, but rather are scaled to
597. %   make the entries the smallest integers possible. Lower values of NICE
598. %   lead to smaller entries, but in general take longer to compute (and may
599. %   not be able to be computed at all, if NICE is too small).
600. %
601. %   URL: http://www.njohnston.ca/publications/minimum-upbs/code/
602. %
603. %   References:
604. %   [CJ] J. Chen and N. Johnston. The Minimum Size of Unextendible Product
605. %        Bases in the Bipartite Case (and Some Multipartite Cases).
606. %        Preprint, 2012.
607. 
608. %   requires: IsTotallyNonsingular.m, normalize_cols.m
609. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
610. %   version: 1.00
611. %   last updated: December 19, 2012
612. 
613. function W = CJ_Lemma5(q,r,s,varargin)
614. 
615.     if(q < r + s)
616.         error('CJ_Lemma5:InvalidArgs','The arguments must satisfy q >= r + s.');
617.     end
618. 
619.     if(nargin == 3)
620.         nice = 0;
621.     else
622.         nice = varargin{1};
623.     end
624. 
625.     its = 0;
626.     satisfy_cond_b = false;
627.     while ~satisfy_cond_b
628.         % Construct a matrix W that satisfies condition (a) of Lemma 4.
629.         if(nice == 0)
630.             W = randn(r+1,2*q);
631.             W(:,q+1:(2*q)) = normalize_cols(W(:,q+1:(2*q)));
632.         else
633.             W = sym((1-2*floor(2*rand(r+1,2*q))).*(1+floor(nice*rand(r+1,2*q))));
634.         end
635.         lin_depen = 0;
636. 
637.         % Now generate the remaining conditions according to condition (b).
638.         for j = 0:q-1
639.             tmp_null = null(W(:,q+1+mod(j+(s:s+r-1),q)).');
640.             if(size(tmp_null,2) > 1)
641.                 % Uh-oh, two of the randomly-generated columns are linearly
642.                 % dependent: try again!
643.                 lin_depen = 1;
644.                 break;
645.             end
646.             W(:,j+1) = tmp_null;
647. 
648.             % Get rid of the denominators in the new column, if the user wanted
649.             % all integers.
650.             if(nice > 0)
651.                 [~,cd] = numden(W(1,j+1));
652.                 for k=2:size(W,1)
653.                     [~,den] = numden(W(k,j+1));
654.                     cd = lcm(cd,den);
655.                 end
656. 
657.                 W(:,j+1) = cd*W(:,j+1);
658.             end
659.         end
660. 
661.         % With probability 1, the matrix W also satisfies condition (b) of
662.         % Lemma 4. Nevertheless, we should make sure that this is the case, and
663.         % if it isn't, try again!
664.         if(lin_depen == 0)
665.             W = fliplr(W);
666.             satisfy_cond_b = IsTotallyNonsingular(double(W),r+1);
667.         end
668. 
669.         its = its + 1;
670.         if(its == 10 && nice ~= 0 && ~satisfy_cond_b)
671.             warning('CJ_Lemma5:LinearDependence', 'Tried (and failed) 10 times to find a matrix satisfying all imposed conditions. You may have better luck if you increase (or omit) the NICE argument.');
672.         end
673.     end
674. end
675. 
676. %%  CJ_LEMMA6    Construct a matrix W described by Lemma 6 of [CJ]
677. %   This function has one required argument:
678. %     K: the positive integer k described by the lemma
679. %
680. %   W = CJ_Lemma6(K) is a (4K+1)-by-(4K+4) matrix with columns of unit
681. %   length that satisfy the orthogonality conditions (i), (ii) and (iii)
682. %   and the nonsingularity condition (iv) of Lemma 6.
683. %
684. %   This function has one optional argument:
685. %     NICE (default 0): a non-negative integer that specifies how "nice"
686. %                       the entries of W should be
687. %
688. %   W = CJ_Lemma6(K,NICE) is the same as before if NICE = 0. If NICE is a
689. %   positive integer, then W will be a symbolic matrix (rather than a
690. %   numeric matrix) whose entries are integers. In this case, the columns
691. %   of W are no longer scaled to have unit length, but rather are scaled to
692. %   make the entries the smallest integers possible. Lower values of NICE
693. %   lead to smaller entries, but in general take longer to compute (and may
694. %   not be able to be computed at all, if NICE is too small).
695. %
696. %   URL: http://www.njohnston.ca/publications/minimum-upbs/code/
697. %
698. %   References:
699. %   [CJ] J. Chen and N. Johnston. The Minimum Size of Unextendible Product
700. %        Bases in the Bipartite Case (and Some Multipartite Cases).
701. %        Preprint, 2012.
702. 
703. %   requires: IsTotallyNonsingular.m, normalize_cols.m
704. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
705. %   version: 1.00
706. %   last updated: December 21, 2012
707. 
708. function W = CJ_Lemma6(k,varargin)
709. 
710.     d = 4*k + 1;
711. 
712.     if(nargin == 1)
713.         nice = 0;
714.     else
715.         nice = varargin{1};
716.     end
717. 
718.     its = 0;
719.     satisfy_cond_b = false;
720.     while ~satisfy_cond_b
721.         % Construct a matrix W that satisfies conditions (i), (ii) and (iii) of
722.         % Lemma 6.
723.         if(nice == 0)
724.             W = randn(d,d+3);
725.             W(:,1:2) = normalize_cols(W(:,1:2));
726.         else
727.             W = sym((1-2*floor(2*rand(d,d+3))).*(1+floor(nice*rand(d,d+3))));
728.         end
729. 
730.         for j = 2:d+2
731.             hlf = (d+3)/2;
732.             if(j < hlf)
733.                 ind = setdiff(0:j-2,mod(j+1,hlf))+1;
734.             else
735.                 ind = setdiff(0:j-1,[j-hlf,hlf+mod(j-1,hlf),hlf+mod(j+1,hlf)])+1;
736.             end
737.             tmp = null(W(:,ind).');
738.             W(:,j+1) = tmp(:,1);
739. 
740.             % Get rid of the denominators in the new column, if the user wanted
741.             % all integers.
742.             if(nice > 0)
743.                 [~,cd] = numden(W(1,j+1));
744.                 for k=2:size(W,1)
745.                     [~,den] = numden(W(k,j+1));
746.                     cd = lcm(cd,den);
747.                 end
748. 
749.                 W(:,j+1) = cd*W(:,j+1);
750.             end
751.         end
752. 
753.         % With probability 1, the matrix W also satisfies condition (iv) of
754.         % Lemma 6. Nevertheless, we should make sure that this is the case, and
755.         % if it isn't, try again!
756.         satisfy_cond_b = IsTotallyNonsingular(double(W),d);
757. 
758.         its = its + 1;
759.         if(its == 10 && nice ~= 0 && ~satisfy_cond_b)
760.             warning('CJ_Lemma6:LinearDependence', 'Tried (and failed) 10 times to find a matrix satisfying all imposed conditions. You may have better luck if you increase (or omit) the NICE argument.');
761.         end
762.     end
763. end
764. 
765. %%  HOLLOWUNITARY    Produces a hollow unitary matrix
766. %   This function has one required argument:
767. %     DIM: a positive integer not equal to 1 or 3
768. %
769. %   U = HollowUnitary(DIM) is a unitary matrix with all of its diagonal
770. %   entries equal to 0 and all of its off-diagonal entries non-zero. Such a
771. %   matrix exists if and only if DIM = 2 or DIM >= 4.
772. %
773. %   URL: http://www.njohnston.ca/publications/minimum-upbs/code/
774. %
775. %   References:
776. %   [CJ] J. Chen and N. Johnston. The Minimum Size of Unextendible Product
777. %        Bases in the Bipartite Case (and Some Multipartite Cases).
778. %        Preprint, 2012.
779. 
780. %   requires: FourierMatrix.m
781. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
782. %   version: 1.00
783. %   last updated: December 19, 2012
784. 
785. function U = HollowUnitary(dim)
786. 
787.     if(dim == 3 || dim <= 1)
788.         error('HollowUnitary:InvalidDim','Hollow unitaries only exist when DIM = 2 or DIM >= 4.');
789.     end
790. 
791.     w2 = exp(2i*pi/(dim-2));
792.     F = FourierMatrix(dim-1);
793. 
794.     U2 = zeros(dim-1);
795.     for j = 1:dim-2
796.         U2 = U2 + w2^j*F(:,j+1)*F(:,j+1)';
797.     end
798. 
799.     U = [0,ones(1,dim-1)/sqrt(dim-1);ones(dim-1,1)/sqrt(dim-1),U2];
800. end

## References

1. C.H. Bennett, D.P. DiVincenzo, T. Mor, P.W. Shor, J.A. Smolin, and B.M. Terhal. Unextendible product bases and bound entanglement. Phys. Rev. Lett. 82, 5385–5388, 1999. E-print: arXiv:quant-ph/9808030
2. D.P. DiVincenzo, T. Mor, P.W. Shor, J.A. Smolin, and B.M. Terhal. Unextendible product bases, uncompletable product bases and bound entanglement. Commun. Math. Phys. 238, 379–410, 2003. E-print: arXiv:quant-ph/9908070
3. T.B. Pedersen. Characteristics of unextendible product bases. Master's Thesis, Aarhus Universitet, Datalogisk Institut, 2002.