# Efficient algorithm to generate a basis for exact diagonalization

The problem is described as follows:

I need to generate a basis(matrix) in lexicographic order.

For two different basis vector $\{n_1,n_2,\cdots,n_M\}$ and $\{t_1,t_2,\cdots,t_M\}$, there must exist a certain index $1\leq k \leq M-1$ such that $n_i=t_i$ for $1\leq i \leq k-1$, while $n_k \neq t_k$. We say $\{n_1,n_2,\cdots,n_M\}$ is superior to $\{t_1,t_2,\cdots,t_M\}$ if $n_k>t_k$.

The whole algorithm can be described as follows, supposing $n_1+n_2+\cdots+n_M=N$,

Starting from $\{N,0,0,\cdots,0\}$, supposing $n_k \neq 0$ while ni=0 for all $k+1 \leq i \leq M-1$, then the next basis vector is $\{t_1,t_2,\cdots,t_M\}$ with

(1) $t_i=n_i$, for $1\leq i \leq k-1$;

(2) $t_k=n_k -1$;

(3) $t_{(k+1)}$=N-Sum[ni,{i,1,k}]

(4) $t_i=0$ for $i\geq k+2$

The generating procedure shall stop at the final vector {0,0,...,0,N}

My code is presented as follows

baseGenerator[M_Integer, N_Integer] :=
NestList[
Block[{k, base},
k = Part[{M}-FirstPosition[Reverse[#[[1 ;; M - 1]]], x_ /; x > 0], 1];
base = #;
base[[k]] = #[[k]] - 1;
base[[k + 1]] = N - Total[#[[1 ;; k]]] + 1;
base[[k + 2 ;; M]] = 0;
base] &,
Join[{N}, ConstantArray[0, M - 1]], (N + M - 1)!/(N!*(M - 1)!) - 1]


Each line corresponds to a specfic step described above. I don't think my code is efficient since it cost much more time than MATLAB writing with loop. Is there any method to improve it a lot ?

In fact, most of the time is costed when calculating position $k$ with command FirstPosisiton.

• AchillesJJ: Thanks for accept, glad it helped - and I hope my Matlab comment was not taken the wrong way - the style of the OP code makes perfect sense in the iterative/matrix oriented world of Matlab, and no slight on your proficiency intended. I simply meant that like Matlab, MMA has idioms that can make a huge performance difference, and vice versa...
– ciao
Jun 17, 2015 at 6:27

baseGenerator2[m_Integer, n_Integer] :=

And for your own sanity, don't use uppercase initials on symbols - you may very well clash with built-ins and/or create debugging nightmares (e.g. N is a built-in, by happenstance you did not have an issue there).