追記

上のコードはortho_basisで基底ベクトルの組を、配列をヒープで作って返しているんだけどこれはタプルで十分ですね。
受け側ではmatch〜withでバインドする。
この変更を行ったところ約9.4秒→約8.3秒程度と高速化しました。
以下変更部分ソース

let ortho_basis(normal:Vec) =
    let tmp = 
        if (normal.x < 0.6) && (normal.x > -0.6) then
            Vec(1.0,0.0,0.0)
        elif (normal.y < 0.6) && (normal.y > -0.6) then
            Vec(0.0,1.0,0.0)
        elif (normal.z < 0.6) && (normal.z > -0.6) then
            Vec(0.0,0.0,1.0)
        else
            Vec(1.0,0.0,0.0)        
                 
    let v0 = vnormalize(vcross(tmp, normal))
    let v1 = vnormalize(vcross(normal, v0))
    (v0,v1,normal)
    
let ambient_occlusion(p:Vec, n:Vec) =
    let ntheta = NAO_SAMPLES
    let nphi   = NAO_SAMPLES
    let eps = 0.0001
    let rp = Vec(p.x+eps*n.x,p.y+eps*n.y,p.z+eps*n.z) 
    let b0,b1,b2 = match ortho_basis(n) with 
                        |(x,y,z) -> x,y,z    
    
    let mutable count = 0
    for j=0 to ntheta-1 do
        for i=0 to nphi-1 do
            let theta = Math.Sqrt(rnd.NextDouble())
            let phi = 2.0 * Math.PI * rnd.NextDouble()
            let x = Math.Cos(phi) * theta
            let y = Math.Sin(phi) * theta
            let z = Math.Sqrt(1.0 - theta * theta)  
            let dx = x*b0.x + y * b1.x + z * b2.x
            let dy = x*b0.y + y * b1.y + z * b2.y
            let dz = x*b0.z + y * b1.z + z * b2.z
            let rd = Vec(dx,dy,dz)
            let r = Ray(rp,rd)
            let h = intersection(r,10000.0)
            let pc = match h with 
                        | Hit(t,p,n) -> 1
                        | NoHit -> 0
            count <- count + pc
            
    let occlusionRatio = double(ntheta * nphi - count) / (double)(ntheta * nphi)
    Vec(occlusionRatio,occlusionRatio,occlusionRatio)