ray-triangle intersection precision - c++

Question

I'm building a software renderer but I think there might be some problem with my ray-triangle intersection accuracy.

I implemented the algorithm referring to pbrt-v3 without some components which I think have little influence on the intersection test.

Below is the code perform intersection test:

bool Triangle::intersectP(const Ray &r, Point3f p0, Point3f p1, Point3f p2) const {
    // translate
    Point3f p0_ = p0 - r._o;
    Point3f p1_ = p1 - r._o;
    Point3f p2_ = p2 - r._o;

    // permute, to make sure z-dimension is not zero
    // permutation is just like to circular move right/left until z is largest dimemsion.
    int kz = r._d.maxDim();
    int kx = kz+1; if (kx == 3) kx = 0;
    int ky = kx+1; if (ky == 3) ky = 0;

    // swap coord
    Vec3f d_ = r._d.permute(kx, ky, kz);
    p0_ = p0_.permute(kx, ky, kz);
    p1_ = p1_.permute(kx, ky, kz);
    p2_ = p2_.permute(kx, ky, kz);
    // shear
    float Sx = -d_.x()/d_.z();
    float Sy = -d_.y()/d_.z();
    float Sz = 1.0f/d_.z();

    // perform transform (p.z doesn't matter now)
    p0_[0] += Sx*p0_[2];
    p0_[1] += Sy*p0_[2];
    p1_[0] += Sx*p1_[2];
    p1_[1] += Sy*p1_[2];
    p2_[0] += Sx*p2_[2];
    p2_[1] += Sy*p2_[2];

    // signed edge function to judge whether origin's projection in projected triangle
    float e0 = p1_.x()*p2_.y() - p2_.x()*p1_.y();   // p1 -> p2
    float e1 = p2_.x()*p0_.y() - p0_.x()*p2_.y();   // p2 -> p0
    float e2 = p0_.x()*p1_.y() - p1_.x()*p0_.y();   // p0 -> p1

    // perform high-precision at edges
    if (e0 == 0.0f || e1 == 0.0f || e2 == 0.0f) {
        double p0xp1y = (double)p0_.x() * (double)p1_.y();
        double p1xp0y = (double)p1_.x() * (double)p0_.y();
        e2 = (float)(p0xp1y - p1xp0y);
        double p1xp2y = (double)p1_.x() * (double)p2_.y();
        double p2xp1y = (double)p2_.x() * (double)p1_.y();
        e0 = (float)(p1xp2y - p2xp1y);
        double p2xp0y = (double)p2_.x() * (double)p0_.y();
        double p0xp2y = (double)p0_.x() * p2_.y();
        e1 = (float)(p2xp0y - p0xp2y);
    }

    // edge detection
    if ( (e0 < 0 || e1 < 0 ||  e2 < 0) && (e0 > 0 || e2 > 0 || e1 > 0) )
        return false;

    float det = e0 + e1 + e2;
    if (det == 0)   // on edge
        return false;

    // compute z value (i.e. t)
    p0_[2] *= Sz;
    p1_[2] *= Sz;
    p2_[2] *= Sz;

    float tScaled = e0*p0_.z() + e1*p1_.z() + e2*p2_.z();
    if (det < 0 && (tScaled >= 0 || tScaled < r._tmax*det) ) return false;
    else if (det > 0 && (tScaled <= 0 || tScaled > r._tmax*det) ) return false;

    return true;
}

Unlike pbrt, my implement use right-hand coordinate and row-major matrices, but I think this do not influence the correctness of this piece of code.

Most of the time, this function works fine. However, artifacts appears when I do occlusion test to the light sources, where r._tmax has a specific value (the distance between intersection point and light sample point).

I tested the algorithm with a plane and a point light above it, however just use the position of the light, radiance is not computed, and this is why I think the problem is caused by intersection implementation. If unoccluded, a value is returned as radiance, otherwise zero.

Vec3f PointEmitter::evalDirect(std::shared_ptr<Scene> scene, const IntersectInfo &isect_info, const Point2f &u) const {
    Vec3f wi;
    float light_pdf;
    VisibilityTester vt;
    Vec3f Li_sampled = sample_Li(isect_info, wi, light_pdf, vt, u);
    Vec3f f = isect_info._bsdf->f(isect_info._wo, wi)*std::abs(isect_info._shading._n.dot(wi));

    return (scene->intersectP(IntersectInfo(_p).spawnRayTo(isect_info))) *100;
}

And another open box-like polygon, with same light: There's a circle-like region obviously stands out on the surface.

I consider this comes from the float point error when comparing r_tmax and the distance. But the code already uses double-precision. Can anyone help?