Faster recursions in sphere decoding

© Erik Agrell, May–Dec. 2010. This software may be copied and modified, if you give proper credit to the author in derived works and publications.

Purpose: The eight sphere decoding algorithms described in [1] are here implemented and tested.

References:
[1] A. Ghasemmehdi and E. Agrell, "Faster recursions in sphere decoding," IEEE Transactions on Information Theory, to appear, 2011.
[2] E. Agrell and T. Eriksson, "Optimization of lattices for quantization," IEEE Transactions on Information Theory, vol. 44, no. 5, pp. 1814-1828, Sept. 1998.

Software: The implementations are intended for use with Mathematica 6.0.

Defining algorithms

In the following, the eight algorithms are first listed interlaced with each other, as in [1, Fig. 2], and they are  defined separately by extracting appropriate lines from this interlaced representation.

Global variables: n, Umin, Umax.

Utilities

sign[x_]:=If[x≤0,-1,1];
round[x_]:=Round[x];
roundc[x_]:=Min[Max[Round[x],Umin],Umax];

Eight algorithms in one figure

This is [1, Fig. 2], translated line by line into Mathematica. It is not an algorithm, but a graphical representation of all eight algorithms, to facilitate comparisons between their designs. It does not make sense to evaluate this representation. As explained in [1], each of the eight algorithms 1,...,8 is defined by extracting the lines labeled with the algorithm's digit.

(*12345678*) C = Infinity;
(*12  5678*) i = n+1;
(*  34    *) i = n;
(*    5678*) d = Table[n,{j,1,n}];
(*12345678*) lambda[[n+1]] = 0;
(*  34  78*) Do[E[[n,j]] = Sum[r[[k]]*H[[k,j]],{k,j,n}], {j,1,n}];
(*    56  *) Do[F[[n,j]] = r[[j], {j,1,n}];
(*  3     *) u[[n]] = round[E[[n,n]]];
(*   4    *) u[[n]] = roundc[E[[n,n]]];
(*  34    *) y = (E[[n,n]]-u[[n]])/H[[n,n]];
(*  34    *) Delta[[n]] = sign[y];
(*  34    *) lambda[[n]] = y^2;
(*12345678*) Label[LOOP];
(*12345678*) While[lambda[[i]] < C,
(*12345678*)   If[i != 1,
(*12345678*)     i = i-1;
(*  34    *)     Do[E[[i,j]] = E[[i+1,j]]-y*H[[i+1,j]], {j,1,i}];
(*    56  *)     Do[F[[j-1,i]] = F[[j,i]]-u[[j]]*G[[j,i]], {j,d[[i]],i+1,-1}];
(*      78*)     Do[E[[j-1,i]] = E[[j,i]]-y[[j]]*H[[j,i]], {j,d[[i]],i+1,-1}];
(*12      *)     p[[i]] = (r[[i]]-Sum[u[[j]]*G[[j,i]],{j,i+1,n}])/G[[i,i]];
(*    56  *)     p[[i]] = F[[i,i]]/G[[i,i]];
(*1   5   *)     u[[i]] = round[p[[i]]];
(* 2   6  *)     u[[i]] = roundc[p[[i]]];
(*  3   7 *)     u[[i]] = round[E[[i,i]]];
(*   4   8*)     u[[i]] = roundc[E[[i,i]]];
(*12  56  *)     y = (p[[i]]-u[[i]])*G[[i,i]];
(*  34    *)     y = (E[[i,i]]-u[[i]])/H[[i,i]];
(*      78*)     y[[i]] = (E[[i,i]]-u[[i]])/H[[i,i]];
(*123456  *)     Delta[[i]] = sign[y];
(*      78*)     Delta[[i]] = sign[y[[i]]];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*      78*)     lambda[[i]] = lambda[[i+1]]+y[[i]]^2;
(*12345678*)   , (* else *)
(*12345678*)     uhat = u;
(*12345678*)     C = lambda[[1]];
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) m = i;
(*12345678*) While[lambda[[i]] >= C,
(*12345678*)   If[i == n,
(*12345678*)     Return[uhat];
(*12345678*)   , (* else *)
(*12345678*)     i = i+1;
(* 2 4 6  *)     y = Infinity;
(*       8*)     y[[i]] = Infinity;
(*12345678*)     u[[i]] = u[[i]]+Delta[[i]];
(*12345678*)     Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)     If[Umin <= u[[i]] <= Umax,
(*12  56  *)       y = (p[[i]]-u[[i]])*G[[i,i]];
(*  34    *)       y = (E[[i,i]]-u[[i]])/H[[i,i]];
(*      78*)       y[[i]] = (E[[i,i]]-u[[i]])/H[[i,i]];
(* 2 4 6 8*)     , (* else *)
(* 2 4 6 8*)       u[[i]] = u[[i]]+Delta[[i]];
(* 2 4 6 8*)       Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)       If[Umin <= u[[i]] <= Umax,
(* 2   6  *)         y = (p[[i]]-u[[i]])*G[[i,i]]];
(*   4    *)         y = (E[[i,i]]-u[[i]])/H[[i,i]]];
(*       8*)         y[[i]] = (E[[i,i]]-u[[i]])/H[[i,i]]];
(* 2 4 6 8*)     ]; (* end If *)
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*      78*)     lambda[[i]] = lambda[[i+1]]+y[[i]]^2;
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) Do[d[[j]] = i, {j,m,i-1}];
(*    5678*) Do[
(*    5678*)   If[d[[j]] < i,
(*    5678*)     d[[j]] = i;
(*    5678*)   , (* else *)
(*    5678*)     Goto[LOOP]];
(*    5678*) , {j,m-1,1,-1}]; (* end Do *)
(*12345678*) Goto[LOOP];

Algorithm 1

sd1[G_,r_]:=Block[{
lambda=Table[,{n+1}],
u=Table[,{n}],
uhat=Table[,{n}],
Delta=Table[,{n}],
p=Table[,{n}],
y,C,i,j
},
(*12345678*) C = Infinity;
(*12  5678*) i = n+1;
(*12345678*) lambda[[n+1]] = 0;
(*12345678*) Label[LOOP];
(*12345678*) While[lambda[[i]] < C,
(*12345678*)   If[i != 1,
(*12345678*)     i = i-1;
(*12      *)     p[[i]] = (r[[i]]-Sum[u[[j]]*G[[j,i]],{j,i+1,n}])/G[[i,i]];
(*1   5   *)     u[[i]] = round[p[[i]]];
(*12  56  *)     y = (p[[i]]-u[[i]])*G[[i,i]];
(*123456  *)     Delta[[i]] = sign[y];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   , (* else *)
(*12345678*)     uhat = u;
(*12345678*)     C = lambda[[1]];
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*12345678*) While[lambda[[i]] >= C,
(*12345678*)   If[i == n,
(*12345678*)     Return[uhat];
(*12345678*)   , (* else *)
(*12345678*)     i = i+1;
(*12345678*)     u[[i]] = u[[i]]+Delta[[i]];
(*12345678*)     Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(*12  56  *)       y = (p[[i]]-u[[i]])*G[[i,i]];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*12345678*) Goto[LOOP]
];

Algorithm 2

sd2[G_,r_]:=Block[{
lambda=Table[,{n+1}],
u=Table[,{n}],
uhat=Table[,{n}],
Delta=Table[,{n}],
p=Table[,{n}],
y,C,i,j
},
(*12345678*) C = Infinity;
(*12  5678*) i = n+1;
(*12345678*) lambda[[n+1]] = 0;
(*12345678*) Label[LOOP];
(*12345678*) While[lambda[[i]] < C,
(*12345678*)   If[i != 1,
(*12345678*)     i = i-1;
(*12      *)     p[[i]] = (r[[i]]-Sum[u[[j]]*G[[j,i]],{j,i+1,n}])/G[[i,i]];
(* 2   6  *)     u[[i]] = roundc[p[[i]]];
(*12  56  *)     y = (p[[i]]-u[[i]])*G[[i,i]];
(*123456  *)     Delta[[i]] = sign[y];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   , (* else *)
(*12345678*)     uhat = u;
(*12345678*)     C = lambda[[1]];
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*12345678*) While[lambda[[i]] >= C,
(*12345678*)   If[i == n,
(*12345678*)     Return[uhat];
(*12345678*)   , (* else *)
(*12345678*)     i = i+1;
(* 2 4 6  *)     y = Infinity;
(*12345678*)     u[[i]] = u[[i]]+Delta[[i]];
(*12345678*)     Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)     If[Umin <= u[[i]] <= Umax,
(*12  56  *)       y = (p[[i]]-u[[i]])*G[[i,i]];
(* 2 4 6 8*)     , (* else *)
(* 2 4 6 8*)       u[[i]] = u[[i]]+Delta[[i]];
(* 2 4 6 8*)       Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)       If[Umin <= u[[i]] <= Umax,
(* 2   6  *)         y = (p[[i]]-u[[i]])*G[[i,i]]];
(* 2 4 6 8*)     ]; (* end If *)
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*12345678*) Goto[LOOP]
];

Algorithm 3

sd3[H_,r_]:=Block[{
E=Table[,{n},{n}],
lambda=Table[,{n+1}],
u=Table[,{n}],
uhat=Table[,{n}],
Delta=Table[,{n}],
y,C,i,j,k
},
(*12345678*) C = Infinity;
(*  34    *) i = n;
(*12345678*) lambda[[n+1]] = 0;
(*  34  78*) Do[E[[n,j]] = Sum[r[[k]]*H[[k,j]],{k,j,n}], {j,1,n}];
(*  3     *) u[[n]] = round[E[[n,n]]];
(*  34    *) y = (E[[n,n]]-u[[n]])/H[[n,n]];
(*  34    *) Delta[[n]] = sign[y];
(*  34    *) lambda[[n]] = y^2;
(*12345678*) Label[LOOP];
(*12345678*) While[lambda[[i]] < C,
(*12345678*)   If[i != 1,
(*12345678*)     i = i-1;
(*  34    *)     Do[E[[i,j]] = E[[i+1,j]]-y*H[[i+1,j]], {j,1,i}];
(*  3   7 *)     u[[i]] = round[E[[i,i]]];
(*  34    *)     y = (E[[i,i]]-u[[i]])/H[[i,i]];
(*123456  *)     Delta[[i]] = sign[y];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   , (* else *)
(*12345678*)     uhat = u;
(*12345678*)     C = lambda[[1]];
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*12345678*) While[lambda[[i]] >= C,
(*12345678*)   If[i == n,
(*12345678*)     Return[uhat];
(*12345678*)   , (* else *)
(*12345678*)     i = i+1;
(*12345678*)     u[[i]] = u[[i]]+Delta[[i]];
(*12345678*)     Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(*  34    *)       y = (E[[i,i]]-u[[i]])/H[[i,i]];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*12345678*) Goto[LOOP]
];

Algorithm 4

sd4[H_,r_]:=Block[{
E=Table[,{n},{n}],
lambda=Table[,{n+1}],
u=Table[,{n}],
uhat=Table[,{n}],
Delta=Table[,{n}],
y,C,i,j,k
},
(*12345678*) C = Infinity;
(*  34    *) i = n;
(*12345678*) lambda[[n+1]] = 0;
(*  34  78*) Do[E[[n,j]] = Sum[r[[k]]*H[[k,j]],{k,j,n}], {j,1,n}];
(*   4    *) u[[n]] = roundc[E[[n,n]]];
(*  34    *) y = (E[[n,n]]-u[[n]])/H[[n,n]];
(*  34    *) Delta[[n]] = sign[y];
(*  34    *) lambda[[n]] = y^2;
(*12345678*) Label[LOOP];
(*12345678*) While[lambda[[i]] < C,
(*12345678*)   If[i != 1,
(*12345678*)     i = i-1;
(*  34    *)     Do[E[[i,j]] = E[[i+1,j]]-y*H[[i+1,j]], {j,1,i}];
(*   4   8*)     u[[i]] = roundc[E[[i,i]]];
(*  34    *)     y = (E[[i,i]]-u[[i]])/H[[i,i]];
(*123456  *)     Delta[[i]] = sign[y];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   , (* else *)
(*12345678*)     uhat = u;
(*12345678*)     C = lambda[[1]];
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*12345678*) While[lambda[[i]] >= C,
(*12345678*)   If[i == n,
(*12345678*)     Return[uhat];
(*12345678*)   , (* else *)
(*12345678*)     i = i+1;
(* 2 4 6  *)     y = Infinity;
(*12345678*)     u[[i]] = u[[i]]+Delta[[i]];
(*12345678*)     Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)     If[Umin <= u[[i]] <= Umax,
(*  34    *)       y = (E[[i,i]]-u[[i]])/H[[i,i]];
(* 2 4 6 8*)     , (* else *)
(* 2 4 6 8*)       u[[i]] = u[[i]]+Delta[[i]];
(* 2 4 6 8*)       Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)       If[Umin <= u[[i]] <= Umax,
(*   4    *)         y = (E[[i,i]]-u[[i]])/H[[i,i]]];
(* 2 4 6 8*)     ]; (* end If *)
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*12345678*) Goto[LOOP]
];

Algorithm 5

sd5[G_,r_]:=Block[{
F=Table[,{n},{n}],
lambda=Table[,{n+1}],
u=Table[,{n}],
uhat=Table[,{n}],
Delta=Table[,{n}],
p=Table[,{n}],
d=Table[,{n}],
y,C,i,j,m
},
(*12345678*) C = Infinity;
(*12  5678*) i = n+1;
(*    5678*) d = Table[n,{j,1,n}];
(*12345678*) lambda[[n+1]] = 0;
(*    56  *) Do[F[[n,j]] = r[[j]], {j,1,n}];
(*12345678*) Label[LOOP];
(*12345678*) While[lambda[[i]] < C,
(*12345678*)   If[i != 1,
(*12345678*)     i = i-1;
(*    56  *)     Do[F[[j-1,i]] = F[[j,i]]-u[[j]]*G[[j,i]], {j,d[[i]],i+1,-1}];
(*    56  *)     p[[i]] = F[[i,i]]/G[[i,i]];
(*1   5   *)     u[[i]] = round[p[[i]]];
(*12  56  *)     y = (p[[i]]-u[[i]])*G[[i,i]];
(*123456  *)     Delta[[i]] = sign[y];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   , (* else *)
(*12345678*)     uhat = u;
(*12345678*)     C = lambda[[1]];
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) m = i;
(*12345678*) While[lambda[[i]] >= C,
(*12345678*)   If[i == n,
(*12345678*)     Return[uhat];
(*12345678*)   , (* else *)
(*12345678*)     i = i+1;
(*12345678*)     u[[i]] = u[[i]]+Delta[[i]];
(*12345678*)     Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(*12  56  *)       y = (p[[i]]-u[[i]])*G[[i,i]];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) Do[d[[j]] = i, {j,m,i-1}];
(*    5678*) Do[
(*    5678*)   If[d[[j]] < i,
(*    5678*)     d[[j]] = i;
(*    5678*)   , (* else *)
(*    5678*)     Goto[LOOP]];
(*    5678*) , {j,m-1,1,-1}]; (* end Do *)
(*12345678*) Goto[LOOP]
];

Algorithm 6

sd6[G_,r_]:=Block[{
F=Table[,{n},{n}],
lambda=Table[,{n+1}],
u=Table[,{n}],
uhat=Table[,{n}],
Delta=Table[,{n}],
p=Table[,{n}],
d=Table[,{n}],
y,C,i,j,m
},
(*12345678*) C = Infinity;
(*12  5678*) i = n+1;
(*    5678*) d = Table[n,{j,1,n}];
(*12345678*) lambda[[n+1]] = 0;
(*    56  *) Do[F[[n,j]] = r[[j]], {j,1,n}];
(*12345678*) Label[LOOP];
(*12345678*) While[lambda[[i]] < C,
(*12345678*)   If[i != 1,
(*12345678*)     i = i-1;
(*    56  *)     Do[F[[j-1,i]] = F[[j,i]]-u[[j]]*G[[j,i]], {j,d[[i]],i+1,-1}];
(*    56  *)     p[[i]] = F[[i,i]]/G[[i,i]];
(* 2   6  *)     u[[i]] = roundc[p[[i]]];
(*12  56  *)     y = (p[[i]]-u[[i]])*G[[i,i]];
(*123456  *)     Delta[[i]] = sign[y];
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   , (* else *)
(*12345678*)     uhat = u;
(*12345678*)     C = lambda[[1]];
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) m = i;
(*12345678*) While[lambda[[i]] >= C,
(*12345678*)   If[i == n,
(*12345678*)     Return[uhat];
(*12345678*)   , (* else *)
(*12345678*)     i = i+1;
(* 2 4 6  *)     y = Infinity;
(*12345678*)     u[[i]] = u[[i]]+Delta[[i]];
(*12345678*)     Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)     If[Umin <= u[[i]] <= Umax,
(*12  56  *)       y = (p[[i]]-u[[i]])*G[[i,i]];
(* 2 4 6 8*)     , (* else *)
(* 2 4 6 8*)       u[[i]] = u[[i]]+Delta[[i]];
(* 2 4 6 8*)       Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)       If[Umin <= u[[i]] <= Umax,
(* 2   6  *)         y = (p[[i]]-u[[i]])*G[[i,i]]];
(* 2 4 6 8*)     ]; (* end If *)
(*123456  *)     lambda[[i]] = lambda[[i+1]]+y^2;
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) Do[d[[j]] = i, {j,m,i-1}];
(*    5678*) Do[
(*    5678*)   If[d[[j]] < i,
(*    5678*)     d[[j]] = i;
(*    5678*)   , (* else *)
(*    5678*)     Goto[LOOP]];
(*    5678*) , {j,m-1,1,-1}]; (* end Do *)
(*12345678*) Goto[LOOP]
];

Algorithm 7

sd7[H_,r_]:=Block[{
E=Table[,{n},{n}],
lambda=Table[,{n+1}],
u=Table[,{n}],
uhat=Table[,{n}],
Delta=Table[,{n}],
d=Table[,{n}],
y=Table[,{n}],
C,i,j,k,m
},
(*12345678*) C = Infinity;
(*12  5678*) i = n+1;
(*    5678*) d = Table[n,{j,1,n}];
(*12345678*) lambda[[n+1]] = 0;
(*  34  78*) Do[E[[n,j]] = Sum[r[[k]]*H[[k,j]],{k,j,n}], {j,1,n}];
(*12345678*) Label[LOOP];
(*12345678*) While[lambda[[i]] < C,
(*12345678*)   If[i != 1,
(*12345678*)     i = i-1;
(*      78*)     Do[E[[j-1,i]] = E[[j,i]]-y[[j]]*H[[j,i]], {j,d[[i]],i+1,-1}];
(*  3   7 *)     u[[i]] = round[E[[i,i]]];
(*      78*)     y[[i]] = (E[[i,i]]-u[[i]])/H[[i,i]];
(*      78*)     Delta[[i]] = sign[y[[i]]];
(*      78*)     lambda[[i]] = lambda[[i+1]]+y[[i]]^2;
(*12345678*)   , (* else *)
(*12345678*)     uhat = u;
(*12345678*)     C = lambda[[1]];
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) m = i;
(*12345678*) While[lambda[[i]] >= C,
(*12345678*)   If[i == n,
(*12345678*)     Return[uhat];
(*12345678*)   , (* else *)
(*12345678*)     i = i+1;
(*12345678*)     u[[i]] = u[[i]]+Delta[[i]];
(*12345678*)     Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(*      78*)       y[[i]] = (E[[i,i]]-u[[i]])/H[[i,i]];
(*      78*)     lambda[[i]] = lambda[[i+1]]+y[[i]]^2;
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) Do[d[[j]] = i, {j,m,i-1}];
(*    5678*) Do[
(*    5678*)   If[d[[j]] < i,
(*    5678*)     d[[j]] = i;
(*    5678*)   , (* else *)
(*    5678*)     Goto[LOOP]];
(*    5678*) , {j,m-1,1,-1}]; (* end Do *)
(*12345678*) Goto[LOOP]
];

Algorithm 8

sd8[H_,r_]:=Block[{
E=Table[,{n},{n}],
lambda=Table[,{n+1}],
u=Table[,{n}],
uhat=Table[,{n}],
Delta=Table[,{n}],
d=Table[,{n}],
y=Table[,{n}],
C,i,j,k,m
},
(*12345678*) C = Infinity;
(*12  5678*) i = n+1;
(*    5678*) d = Table[n,{j,1,n}];
(*12345678*) lambda[[n+1]] = 0;
(*  34  78*) Do[E[[n,j]] = Sum[r[[k]]*H[[k,j]],{k,j,n}], {j,1,n}];
(*12345678*) Label[LOOP];
(*12345678*) While[lambda[[i]] < C,
(*12345678*)   If[i != 1,
(*12345678*)     i = i-1;
(*      78*)     Do[E[[j-1,i]] = E[[j,i]]-y[[j]]*H[[j,i]], {j,d[[i]],i+1,-1}];
(*   4   8*)     u[[i]] = roundc[E[[i,i]]];
(*      78*)     y[[i]] = (E[[i,i]]-u[[i]])/H[[i,i]];
(*      78*)     Delta[[i]] = sign[y[[i]]];
(*      78*)     lambda[[i]] = lambda[[i+1]]+y[[i]]^2;
(*12345678*)   , (* else *)
(*12345678*)     uhat = u;
(*12345678*)     C = lambda[[1]];
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) m = i;
(*12345678*) While[lambda[[i]] >= C,
(*12345678*)   If[i == n,
(*12345678*)     Return[uhat];
(*12345678*)   , (* else *)
(*12345678*)     i = i+1;
(*       8*)     y[[i]] = Infinity;
(*12345678*)     u[[i]] = u[[i]]+Delta[[i]];
(*12345678*)     Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)     If[Umin <= u[[i]] <= Umax,
(*      78*)       y[[i]] = (E[[i,i]]-u[[i]])/H[[i,i]];
(* 2 4 6 8*)     , (* else *)
(* 2 4 6 8*)       u[[i]] = u[[i]]+Delta[[i]];
(* 2 4 6 8*)       Delta[[i]] = -Delta[[i]]-sign[Delta[[i]]];
(* 2 4 6 8*)       If[Umin <= u[[i]] <= Umax,
(*       8*)         y[[i]] = (E[[i,i]]-u[[i]])/H[[i,i]]];
(* 2 4 6 8*)     ]; (* end If *)
(*      78*)     lambda[[i]] = lambda[[i+1]]+y[[i]]^2;
(*12345678*)   ]; (* end If *)
(*12345678*) ]; (* end While *)
(*    5678*) Do[d[[j]] = i, {j,m,i-1}];
(*    5678*) Do[
(*    5678*)   If[d[[j]] < i,
(*    5678*)     d[[j]] = i;
(*    5678*)   , (* else *)
(*    5678*)     Goto[LOOP]];
(*    5678*) , {j,m-1,1,-1}]; (* end Do *)
(*12345678*) Goto[LOOP]
];

Running some examples

The algorithms are here tested in various scenarios. We demonstrate that algorithms 1, 3, 5, and 7 give the same results and that algorithms 2, 4, 6, and 8 also give the same results. Furthermore, we calculate some lattice properties.

Global variables: n, Umin, Umax

Additional definitions

The algorithms need as input a generator matrix in lower triangular form with positive diagonal elements. The built-in QRDecomposition needs to be modified to yield R with positive diagonal elements:

myQRDecomposition[x_?MatrixQ]:=Block[{Q,R,signs},
{Q,R}=QRDecomposition[x];
signs=Sign[Diagonal[R]];
{Q*signs,R*signs}
];

The Gram matrix of a given generator matrix:

gram[x_?MatrixQ]:=x.Transpose[x];

The squared Euclidean distance:

dist2[x_List,y_List]:=(x-y).(x-y);

Random 2-dimensional lattice

n=2;
SeedRandom[1];
G1=Table[RandomReal[NormalDistribution[0,1]],{n},{n}];

Rotate the lattice into lower triangular form, and verify that the rotated generator matrix has the same Gram matrix as the original generator matrix:

{Q,R}=myQRDecomposition[Transpose[G1]];
G=Transpose[R];
gram[G]==gram[G1]

Calculate the H matrix:

H=Chop[Inverse[G]]

Decode two vectors and show that the four lattice algorithms give the same results:

r1={.5,.9};
{sd1[G,r1],sd3[H,r1],sd5[G,r1],sd7[H,r1]}

r1={.5,1.0};
{sd1[G,r1],sd3[H,r1],sd5[G,r1],sd7[H,r1]}

Random 2-dimensional constellation based on 2-PAM

{n,Umin,Umax}={2,0,1};
SeedRandom[1];
G1=Table[RandomReal[NormalDistribution[0,1]],{n},{n}];

Rotate the lattice into lower triangular form, and verify that the rotated generator matrix has the same Gram matrix as the original generator matrix:

{Q,R}=myQRDecomposition[Transpose[G1]];
G=Transpose[R];
gram[G]==gram[G1]

Calculate the H matrix:

H=Chop[Inverse[G]]

Decode two vectors and show that the four algorithms for finite constellations give the same results:

r1={.5,.9};
{sd2[G,r1],sd4[H,r1],sd6[G,r1],sd8[H,r1]}

r1={.5,1.0};
{sd2[G,r1],sd4[H,r1],sd6[G,r1],sd8[H,r1]}

Random 10-dimensional lattices and constellations

Compare 100 random generator matrices, each with 100 received vectors, using both lattices and finite constellations (4-PAM)

{nG,nr}={100,100};
{n,Umin,Umax}={10,0,3};
SeedRandom[1];

Verify that algorithms 1, 3, 5, and 7 give the same results for all generator matrices and all received vectors, and similarly verify that algorithms 2, 4, 6, and 8 give the same results:

Timing[And@@Table[G1=Table[RandomReal[NormalDistribution[0,1]],{n},{n}];
{q,r}=QRDecomposition[Transpose[G1]];
G=Transpose[r*Sign[Diagonal[r]]];
H=Chop[Inverse[G]];
r1list=Table[RandomReal[NormalDistribution[0,1]],{nr},{n}];
Table[sd1[G,r1],{r1,r1list}]==Table[sd3[H,r1],{r1,r1list}]==Table[sd5[G,r1],{r1,r1list}]==Table[sd7[H,r1],{r1,r1list}]&&
Table[sd2[G,r1],{r1,r1list}]==Table[sd4[H,r1],{r1,r1list}]==Table[sd6[G,r1],{r1,r1list}]==Table[sd8[H,r1],{r1,r1list}],{nG}]]

Voronoi diagram of 2-dimensional constellation based on 2-PAM

Define a 2-dimensional constellation as a subset of the hexagonal lattice:

{n,Umin,Umax}={2,0,1};
G=N[{{1,0},{1/2,Sqrt[3]/2}}];
H=Inverse[G];

Plot 5000 random received vectors, indicating by color which of the four constellation points are closest. This illustrates the Voronoi diagram of the constellation:

Graphics[{
Table[r1={RandomReal[{-.5,2}],RandomReal[{-.5,1.4}]};
{Hue[sd8[H,r1].{.5,.25}],Point[r1]},{5000}],
PointSize[.02],Table[Point[{u1,u2}.G],{u1,Umin,Umax},{u2,Umin,Umax}]
}]

E6 lattice

The following is a square generator matrix for the E6 lattice [2, eq. (38)]:

n=6;
E6=N[{{1,0,0,0,0,a}, {0,1,0,0,0,a},{0,0,1,0,0,a},{0,0,0,1,0,a},{0,0,0,0,1,a},{1,1,1,1,1,3a}/2}/.a->Sqrt[3]];

We wish to find the closest lattice point to the following received vector:

r1={6.,5.,4.,3.,2.,1.};

First, rotate the lattice into lower triangular form, and verify that the rotated generator matrix has the same Gram matrix as the original generator matrix:

{Q,R}=myQRDecomposition[Transpose[E6]];
G=Transpose[R];
gram[G]==gram[E6]

Second, rotate the received vector into the same coordinate system as G:

r=r1.Transpose[Q];
Chop[G.Q]==E6
Chop[r.Q]==r1

Third, decode r in the lattice G:

u=sd7[Inverse[G],r1]

Verify that the other algorithms yield the same results:

u==sd1[G,r1]==sd3[Inverse[G],r1]==sd5[G,r1]

Conclusions:
(1) The closest point to r in the lattice spanned by G is u.G.
(2) The closest point to r1 in the lattice spanned by E6 is u.E6.

Norm[r-u.G]
Norm[r1-u.E6]

Monte Carlo integration over a lattice Voronoi cell

We wish to estimate the second moment of the Voronoi cell of D4.

n=4;
D4=N[{{2,0,0,0},{1,1,0,0},{1,0,1,0},{1,0,0,1}}];
H=Inverse[D4];
SeedRandom[2];

There is no need for QR decomposition, since the generator matrix is already lower triangular with positive diagonal elements.

H=Inverse[D4];

Generate 10 000 vectors uniformly inside the fundamental polytope [2], decode them, and calculate the mean square decoding error:

nsamples=10000;
Timing[mse=Sum[With [{x=Table[RandomReal[{0,1}],{n}].D4},
dist2[x,sd7[H,x].D4]],{nsamples}]/nsamples]

This is the second moment of the Voronoi cell. Normalize by the dimension and the volume of the Voronoi cell:

mse/n/Det[gram[D4]]^(1/n)

The true value is 0.076603 [2].