"Randomize" output of a Linear Feedback Shift Register for the same taps? I'm using a (Galois) LFSR to sample a large array, ensuring that each entry is only visited once. I simply skip past the entries that exceed the array length.
With the same taps then the array entry a is naturally always followed by b.
However, I would like to be able to modify the output with a seed, causing the output for each seed to be different.
A quick and dirty method I used is:
If the current index is a, the "seed" x is 1 or larger and the max index is l:
Calculate a+x and find the next entry in the sequence from that number (a+x), continue until we have a number b in the range x < b < x + l. Calculate b-x, this is the next index.
In other words I shift the sample from the sequence by x.
This sort of works, but isn't very elegant. Are there other possibilities?

Edit: I added a few extra tags, because the comments revealed that in place of LFSRs the OP is also interested in other methods of quickly generating largish sets of permutations of an array with length up to thousands. Hopefully I did not distort the intent, JL.
 A: I add a couple of suggestions to break the ice. These fall into two groups. One is based on the idea of sticking to LFSRs but tweaking them differently. The other is to use modular permutation polynomials. As the latter is new in the contest I also discuss some implementation tricks that may or may not help.
(Mutated) LFSRs:
The taps in a maximal length LFSR (of $n$ bits) guaranteed to cycle through the bit combos from $1$ to $2^n-1$ without repetitions are determined by so called primitive polynomials $p(x)$ of degree $n$ in the ring $\mathbb{F}_2[x]$. For a fixed $n$ there are
$$
N_n=\frac1n\phi(2^n-1)
$$
distinct such polynomials. Here $\phi$ is the Euler totient function. For example, when $n=12$, we have
$$
2^{12}-1=4095=3^2\cdot5\cdot7\cdot13,
$$
so there are
$$
N_{12}=\frac1{12}(3-1)\cdot3\cdot(5-1)\cdot(7-1)\cdot(13-1)=144
$$
different sets of taps you can use here. This is probably not enough for you,
but if you are interested I can describe more details about how to find them.

Given that the LFSR outputs all the non-zero vectors of $n$ bits in sequence you 
can produce more variety by applying a bijective mapping $f:\{0,1\}^n\to\{0,1\}^n$ to all the outputs (still discarding those outputs that fall out of range). It would be relatively simple to use any non-singular linear mapping as $f$: you can specify the images $f(000\ldots001)$,  $f(000\ldots010)$, $f(000\ldots100)$,$\ldots$, $f(100\ldots000)$, and if they are a linearly independent set, you are guaranteed to get a permutation of the non-zero vectors. Calculating the value of $f(k)$ for some bitvector $k$ amounts to taking the bitwise XOR of the images of those bits that are ON in $k$.
As a very simple subset of such functions $f$ you can use a random permutation of the $n$ bits. There are $n!$ ways of permuting all the outputs of the chosen LFSR, which may be a large enough number for you. A potential drawback of sticking to permutations of bits is that a permutation won't change the weight of a binary vector.
A drawback of all LFSR based schemes is that the number $2^n-1$ may be quite a bit larger than the length of your array (in the worst case almost double). Meaning that half the time you will just discard the next entry. The other class
of functions that I discuss is better in that sense.
Modular permutation polynomials:
Here it is simpler to index your array as the range $0,1,\ldots,\ell-1$. I will work in the residue class ring $R=\mathbb{Z}/L\mathbb{Z}$ of integers modulo $L$, where $L$ is an integer that must be $\ge$ than the length of your array. We may have $L=\ell$, but I need to place some requirements on $L$, so I want a bit of freedom here. Anyway, I suggest the use of polynomial mappings
$p:R\to R$. Here $p(x)$ is a polynomial in $x$ with coefficient in $R$ (or just plain integers). There are results of elementary number theory telling us, when such a polynomial gives rise to a permutation of the elements of $R$.
The simplest case is that of linear polynomials (all the calculations are done modulo $L$)
$$
p(x)=ax+b.
$$
This is a permutation, if and only if $\gcd(a,L)=1$. There are no requirements on $b$, it can be anything in the range $0\le b<L$. Altogether there are $L\phi(L)$ such permutations. Implementing such a function is very fast. You simple select a random $b$, set $p(0)=b$, and from that point on use the formula
$$
p(x+1)=p(x)+a,$$
or, if you prefer pseudocode
$$p(x+1)=(p(x)+a)\bmod L.
$$
So using a modular linear polynomial you will start from a random point, and then use a fixed jump length (coprime to the length of the array to guarantee that you will visit all the entries without repetitions). As you expressed a desire to avoid some correlations, this may not be a good choice. 
A slightly more random polynomial would be quadratic. So let us take a look at polynomials of the form
$$
p(x)=ax^2+bx+c.
$$
It is not difficult to show that this is a permutation of $R$, if (this is an "if and only if" for almost all practical purposes) the following two requirements are met: 


*

*$\gcd(b,L)=1$, and

*the coefficient $a$ must be divisible by all the prime numbers $p$ that are factors of $L$.


The second requirement places a serious constraint on us, as we want to use a non-zero $a$ (so no multiple of $L$ will do). For this to be possible, we need $L$ to be divisible by a square of some prime. And even then we still have relatively few choices.
The best case may be to select $L$ to be a square itself, when we can let $a$ to be any multiple of $\sqrt{L}$ (if $L$ is divisible by a fourth power, then we get even more room to play). If $\ell\approx10000$, then the number of "dummy" entries could be up to $200$. Still better than with LFSRs in a worst case.
Another possibility would be to select $L$ that is a multiple of $2^6=64$. If 
$L=2^6M$, then we can use $a=kM$ for any $k$ such that $1\le k<64$. The number of choices for the coefficient $a$ is $\sqrt{L}$ in the first case and $63$ in the second.
Implementing a quadratic polynomial is also easy. We have $p(0)=c$, $p(1)=a+b+c$ (again modulo $L$). In addition to $p(x)$ let us use the difference 
$$\Delta(x):=p(x+1)-p(x).$$ So we know that $\Delta(0)=p(1)-p(0)=a+b$. In general we can calculate that
$$
\Delta(x+1)-\Delta(x)=p(x+2)-2p(x+1)+p(x)=\cdots=2a.
$$
So if you initialize $p(0)=c$, $\Delta(0)=a+b$, and update these two numbers
according to the rules
$$
p(x+1)=(p(x)+\Delta(x))\bmod L
$$
and
$$
\Delta(x+1)=(\Delta(x)+2a)\bmod L,
$$
you only need to allocate memory for two integers to reproduce the entire permutation.
The update rule of $\Delta$ tells us that with a quadratic permutation polynomial the lengths of the jumps vary according to the choice of $a$. The number of different jump lengths will be $L/\gcd(2a,L)$ so in the two example cases $\sqrt{L}/2$ (resp. $32$). A bit better than with linear polynomials, but I'm not sure if that is good enough to kill the correlations you want to avoid.
A: If I understood correctly your question, you do not want to have a random generator, but a family of random generators, so that you can extract from each of them a single permutation of your array. Therefore, you need to parametrize you random generator and produce different streams of random numbers. This is easily feasible with some random generators, but I'm not sure about the LFSR. 
I have been using a family of LCG generators which performs quite well (and the quality of the random numbers is not that bad). I found it in the library SPRNG, which has some nice generators. You can find many informations in the documentation (->Version 4.0 -> User Guid -> Generators).
You should probably use one of the generators mentioned there, since they have been already studied.
Hope it helps!
