Not sure what you want but if you choose to use Mathematica, here's the code I'd use to numerically integrate the function over the upper hemisphere of the unit sphere.

Code:

NIntegrate[f[x, y, Sqrt[1 - x^2 - y^2]]*
Sqrt[1/(1 - x^2 - y^2)], {x, -1, 1},
{y, -Sqrt[1 - x^2], Sqrt[1 - x^2]}]

I realize not everyone has Mathematica or even wants to use it. Just sayin' if you're interested in a brute force approach to the answer using Mathematica, which I'd think would be rather quick (at 2GHz) for a reasonable function, that's what I'd start with. Here's a check of over the upper hemisphere of the unit sphere, the answer of which is , just to get some indication if the coding is correct:

Code:

f[x_, y_, z_] := x^2;
NIntegrate[f[x, y, Sqrt[1 - x^2 - y^2]]*
Sqrt[1/(1 - x^2 - y^2)], {x, -1, 1},
{y, -Sqrt[1 - x^2], Sqrt[1 - x^2]}]
N[2*(Pi/3)]
2.0944
2.0944

That took lest than 1/10 of a second. Can you post your function?