#if !defined HAMMON_DIFFUSE_LIB
#define HAMMON_DIFFUSE_LIB
    #ifdef FAST_HAMMON
        vec3 hammonDiffuse(in vec3 albedo, in vec3 n, in float nDotV, in float nDotL, in float nDotH, in float lDotV, in float roughness) {
            nDotH = abs(nDotH) + 1e-5;
            float facing = 0.5 + 0.5 * lDotV;
            float rough = nDotH <= 0.0 ? 0.0 : facing * (0.9 - 0.4 * facing) * ((0.5 + nDotH) / max(nDotH, 0.05));
            vec3 fresnel_v = 1.0 - fresnelFunction(nDotV, n, vec3(0.0), 0.0);
            vec3 fresnel_l = 1.0 - fresnelFunction(nDotL, n, vec3(0.0), 0.0);
            vec3 energyConservationFactor = vec3(1.0);
            {
                energyConservationFactor -= fresnelFunction(1.0 - randNextF(), n, vec3(0.0), 0.0) * (1.0 - randNextF());
                energyConservationFactor -= fresnelFunction(1.0 - randNextF(), n, vec3(0.0), 0.0) * (1.0 - randNextF());
                energyConservationFactor -= fresnelFunction(1.0 - randNextF(), n, vec3(0.0), 0.0) * (1.0 - randNextF());
                energyConservationFactor -= fresnelFunction(1.0 - randNextF(), n, vec3(0.0), 0.0) * (1.0 - randNextF());
            }
            vec3 smooth_v = (fresnel_l * fresnel_v) / energyConservationFactor;
            vec3 single = mix(smooth_v, vec3(rough), roughness) * rcp(pi);
            float multi = 0.1159 * roughness;

            return max((albedo * single + albedo * multi) * nDotL, 0.0);
        }
    #else
        vec3 hammonDiffuse(in vec3 albedo, in vec3 n, in float nDotV, in float nDotL, in float nDotH, in float lDotV, in float roughness) {
            nDotH = abs(nDotH) + 1e-5;
            float facing = 0.5 + 0.5 * lDotV;
            float rough = nDotH <= 0.0 ? 0.0 : facing * (0.9 - 0.4 * facing) * ((0.5 + nDotH) / max(nDotH, 0.05));
            vec3 fresnel_v = 1.0 - fresnelFunction(nDotV, n, vec3(0.0), 0.0);
            vec3 fresnel_l = 1.0 - fresnelFunction(nDotL, n, vec3(0.0), 0.0);
            vec3 energyConservationFactor = vec3(1.0);
            uint N = 64u;
            for(uint i = 0u; i < N; ++i) {
                vec2 hammersley = hammersley2d(i, N);
                energyConservationFactor -= fresnelFunction(1.0 - hammersley.x, n, vec3(0.0), 0.0) / N;
            }
            vec3 smooth_v = (fresnel_l * fresnel_v) / energyConservationFactor;
            vec3 single = mix(smooth_v, vec3(rough), roughness) * rcp(pi);
            float multi = 0.1159 * roughness;

            return max((albedo * single + albedo * multi) * nDotL, 0.0);
        }
    #endif
#endif