How do I make sure my cosine-weighted sampling works fine? I'm working on a pathtracer and while it's working mostly fine, I'm getting different outputs with cosine-weighted sampling and uniform sampling : 



Uniform sampling is on the left, cosine-weighted is on the right. The difference is somewhat subtle, but it does not go away with more samples so it's not some kind of start up bias.

I've checked my code against PBRT and Wikipedia but I cannot figure out what is wrong, because all the factors seem to be there (brdf, pdf, Li, and cos(theta)). I'm hoping someone here could take a fresh look at it :

// Sampling functions :
Vec3 sample_uniform_hemisphere(Sampler* sampler, float* pdf) {
    float cosTheta, v;
    sample_unit_square(sampler, &amp;cosTheta, &amp;v, NULL);
    float sinTheta = sqrtf(1 - cosTheta * cosTheta);
    float cosPhi = cosf(2 * M_PI * v);
    float sinPhi = sinf(2 * M_PI * v);

    if (pdf)
        *pdf = .5f / M_PI;
    return (Vec3) {
        sinTheta * cosPhi,
        sinTheta * sinPhi,
        cosTheta
    };
}

Vec3 sample_cosine_weighted_hemisphere(Sampler* sampler, float* pdf) {
    Vec3 sample;
    sample_unit_disc(sampler, &amp;sample.x, &amp;sample.y, NULL);

    float z2 = 1.0f - sample.x * sample.x - sample.y * sample.y;
    sample.z = sqrtf(z2 &gt; 0.0f ? z2 : 0.0f);
    if (pdf)
        *pdf = sample.z / M_PI; // p(theta, phi) = cos(theta) / pi
    return sample;
}

// Snippet that does the actual shading :
        float pdf;
        Vec3 bounce_sample = it.material-&gt;bsdf-&gt;sampler(it.sampler, local_out, &amp;pdf);
        Ray bounce = {
        .o = surface_point,
        .d = basis2world(bounce_sample, u, v, it.normal)
        };

        Intersect bounce_it = trace_ray(sc, bounce, it.depth + 1, it.sampler);
        Color f = it.material-&gt;bsdf-&gt;f(it.material-&gt;params,
                       bounce_sample, local_out);
        Color dc = cmul(
            bounce_it.outgoing_radiance,
            f);
        float inv_pdf = pdf &gt; EPSILON ? (1.0f / pdf) : 1.0f;
        rval = cadd(rval,
            cscale(dc, inv_pdf * bounce_sample.z));


Thanks in advance!
The code you've shown here looks right to me. In order to verify the implementations, you could try writing some code that just integrates a few test functions over the hemisphere using both methods and verify that you get the same result (within some sampling error). Also if you haven't already, verify that the outputs from sample_unit_square and sample_unit_disc are doing what they're supposed to.