Random directions on hemisphere oriented by an arbitrary vector Hy, i'm writing a raytracer, and for that I need to generate n random vectors that are inside an hemisphere oriented by the surface normal. Ideally, I would also like being able to restrict the rays so that the angle with the normal is <= some alpha. I don't require mathematics strictness, it can be an approximation.
What I have been doing is creating an arbitrary tangent frame, and a random vector where the x and y coordinate go from -1 to 1 and the z coordinate from 0 to 1, and then multiplying it with the tangent frame, but for some reason this is making the rays to go along a direction, instead of covering the entire hemisphere.
Here's some pseudo code, based on what i'm using, but removed the HLSL clumsy things:
dir = normal

ray.x = GetRandomFloat(-1, 1)
ray.y = GetRandomFloat(-1, 1)
ray.z = GetRandomFloat(0, 1)

N = abs(normal)

if( N.z > N.x and N.z > N.y )
    rt = (1, 0, 0)
else
    rt = (0, 0, 1)

rt = normalize(rt - normal.xzy * dot(rt, In))
rtb = cross(rt, normal.xzy)

rt.xyz = rt.xzy
rtb.xyz = rtb.xzy

ray = normalize(ray.x * rbt + (ray.y * rt + (ray.z * dir)))

If someone can provide some insight into this, I would be grateful. It doesn't need to be code, even a math formula will do. Also, distribution doesn't need to be uniform, just look "good" and 1/0 or 0/0 in some cases are allowed, as they map to INF or NaN, and I can handle that without problems.
EDIT:
While a PDF that fits this is interesting, it will require a random number generator based on a PDF, and that's really expensive, and not actually what I'm looking for.
EDIT 2:
If it can be done, it would be better not to use inverse trigonometric functions, even if it sacrifices uniformity
EDIT 3: Tried Tryss method, but the results are the same, here's a screen to show what I mean:

It is clear that the rays aren't generated along all directions.
I'm using the second method to generate the rays, and then using the same process I had before to orientate it, that is constructing a basis change matrix from the normal and an arbitrary tangent vector. That might be what is wrong, can anyone provide some insight if the code I have posted earlier actually orients the directions from the z up hemisphere to the hemisphere oriented by the vector "normal"
EDIT 4: 
I'll try to make the problem more clear.
I want to generate n random rays, in an hemisphere oriented by an arbitrary direction. I do not require an uniform distribution.
The first part has been answered, taking x = rand(-1,1) y = rand(-1,1) z = rand(0,1) and then normalizing, while not uniform, produces acceptable results and is fast enough.
Now the second part is transforming that vector so that it is oriented by an arbitrary direction, in a way that (0,0,1) -> (nx,ny,nz) , where (nx,ny,nz) is the arbitrary vector.
The coordinate system for the final vector and the arbitrary directions is a left handed system where x is left and right, z is front and back and y is up and down.
I know that this matrix does the transformation I want, because I do the same coordinate change for normal mapping (http://en.wikipedia.org/wiki/Normal_mapping )
(bitangent.x, bitangent.y, bitangent.z)
(tangent.x  , tangent.y  , tangent.z  )
(normal.x   , normal.y   , normal.z   )

Where bitangent is normal X tangent
The problem now is that I only have the "normal" vector, so I need to generate an arbitrary "tangent" vector to construct that matrix. How can I do that, in a way that the rays get properly transformed ( they don't biased to a direction, like in the screen I showed) ?
The code I was using earlier is at the top of the page, but it produces the biasing I want to avoid.
 A: An uniform distribution on a sphere can be calculated like this :


*

*$\theta_0 = 2\pi \text{Rand}()$

*$\theta_1 = \arccos( 1- 2 \text{Rand}())$

*$x = R \sin(\theta_0)\sin(\theta_1)$

*$y = R \sin(\theta_0)\cos(\theta_1)$

*$z = R \sin(\theta_1)$


Source : http://mathproofs.blogspot.fr/2005/04/uniform-random-distribution-on-sphere.html
Another approach that may be more suited to computer is the following algorithm :
do
   x = GetRandomFloat(-1, 1)
   y = GetRandomFloat(-1, 1)
   z = GetRandomFloat(-1, 1)
   d = sqrt(x^2+y^2+z^2)
while(d>1)

x=x/d
y=y/d
z=z/d

And you have your distribution on the sphere. Notice that the number of step may be unbounded, but the probability to not have a correct point after N iterations is less than $\frac{1}{2^N}$, so it's decent : on average only two iterations to get a point on the sphere
A: Have a good idea: 1. generate lot of random points over spherical surface. 2.check cosine value between the normal and each point. If cosine $< 0$, invertir point. And voila!!
A: As I already said in my comment the best choice is the von-Mises/Fisher distribution something similar to the normal distribution on the n-dimensional sphere, see the cited wikipedi article. Since the OP has problems with the creation of such random vectors, here some references:
To simulate random vectors from this distribution see:
http://www.jstor.org/discover/10.2307/2347441?sid=21105943589933&uid=4&uid=2&uid=3737864
or here:
https://stackoverflow.com/questions/4415668/drawing-random-values-from-a-fisher-distribution
or here:
http://luc.devroye.org/rnbookindex.html
@Santiago Pacheco you could have found them all yourself on the net!
A: In order to tackle your second challenge, you may want to use Rodrigues' rotation formula, though it would be computationally intensive (basically aligning - or rotating - your set of randomly generated coordinates/vectors along another vector).
This second link is another answer on this same site for the same question.
