Too bad, in the sense that for smaller L, say on the order of 8 or so, this is not that difficult to do exactly. (I have written code to solve exactly this problem, so I know that to be true.) For a higher number of dimensions, here 300, it gets nasty to do for uniformly distributed samples. Or, if the distribution of the numbers was intended to be normally distributed, it would be again pretty easy.
I might suggest one simple solution, that won't be truly correct in terms of the distribution of the numbers, but it won't be terrible. And it will ensure positively the result will be in the interval [0,1]. You don't even need to use randfixedsum.
L = 300;
y = (rand(1,L) > 0.5)*2 - 1;
yp = y > 0;
yn = ~yp;
a = rand(1,L);
S = a*y';
if S > 0
sumyp = sum(a(yp));
a(yp) = a(yp)*(1 - S/sumyp);
sumyn = sum(a(yn));
a(yn) = a(yn)*(1 + S/sumyn);
The necessary bias will generally not be huge. Did it work?
Well, yes. The requisite sum will be identically zero, except for floating point trash.
The vector a lies always in the interval [0,1]. No exceptions, ever. And it is quite efficient.
The only downside is the necessary (usually tiny) bias applied to roughly half of the elements will make the result not perfectly uniformly distributed. It is not perfect, but possibly acceptable.