Created concepts for samplers, added quotient_and_pdf variants to satisfy the concepts

include/nbl/builtin/hlsl/algorithm.hlsl

@@ -92,12 +92,12 @@ NBL_CONSTEXPR_FUNC void swap(NBL_REF_ARG(T) lhs, NBL_REF_ARG(T) rhs)
 namespace impl
 {
 
-template<class Accessor, class Comparator>
+template<class Accessor, class Comparator, bool IsUpper = false>
 struct bound_t
 {
-    static bound_t<Accessor,Comparator> setup(uint32_t begin, const uint32_t end, const typename Accessor::value_type _value, const Comparator _comp)
+    static bound_t<Accessor,Comparator,IsUpper> setup(uint32_t begin, const uint32_t end, const typename Accessor::value_type _value, const Comparator _comp)
     {
-        bound_t<Accessor,Comparator> retval;
+        bound_t<Accessor,Comparator,IsUpper> retval;
         retval.comp = _comp;
         retval.value = _value;
         retval.it = begin;
@@ -110,7 +110,7 @@ struct bound_t
     {
         if (isNPoT(len))
         {
-            const uint32_t newLen = 0x1u<<nbl::hlsl::findMSB(len);
+            const uint32_t newLen = 0x1u<<_static_cast<uint32_t>(nbl::hlsl::findMSB(len));
             const uint32_t testPoint = it+(len-newLen);
             len = newLen;
             comp_step(accessor,testPoint);
@@ -132,9 +132,17 @@ struct bound_t
         comp_step(accessor,mid);
     }
 
+    bool compare(const typename Accessor::value_type lhs, const typename Accessor::value_type rhs)
+    {
+        if (IsUpper)
+            return !comp(rhs,lhs);
+        else
+            return comp(lhs,rhs);
+    }
+
     void comp_step(NBL_REF_ARG(Accessor) accessor, const uint32_t testPoint, const uint32_t rightBegin)
     {
-        if (comp(accessor[testPoint],value))
+        if (compare(accessor[testPoint],value))
             it = rightBegin;
     }
     void comp_step(NBL_REF_ARG(Accessor) accessor, const uint32_t testPoint)
@@ -148,36 +156,25 @@ struct bound_t
     uint32_t len;
 };
 
-template<class Accessor, class Comparator>
-struct lower_to_upper_comparator_transform_t
-{
-    bool operator()(const typename Accessor::value_type lhs, const typename Accessor::value_type rhs)
-    {
-        return !comp(rhs,lhs);
-    }
-
-    Comparator comp;
-};
-
 }
 
 
 template<class Accessor, class Comparator>
-uint32_t lower_bound(NBL_REF_ARG(Accessor) accessor, const uint32_t begin, const uint32_t end, const typename Accessor::value_type value, const Comparator comp)
+uint32_t lower_bound(NBL_REF_ARG(Accessor) accessor, const uint32_t begin, const uint32_t end, const typename Accessor::value_type value, NBL_REF_ARG(Comparator) comp)
 {
-    impl::bound_t<Accessor,Comparator> implementation = impl::bound_t<Accessor,Comparator>::setup(begin,end,value,comp);
-    return implementation(accessor);
+    impl::bound_t<Accessor,Comparator,false> implementation = impl::bound_t<Accessor,Comparator,false>::setup(begin,end,value,comp);
+    const uint32_t retval = implementation(accessor);
+    comp = implementation.comp;
+    return retval;
 }
 
 template<class Accessor, class Comparator>
-uint32_t upper_bound(NBL_REF_ARG(Accessor) accessor, const uint32_t begin, const uint32_t end, const typename Accessor::value_type value, const Comparator comp)
+uint32_t upper_bound(NBL_REF_ARG(Accessor) accessor, const uint32_t begin, const uint32_t end, const typename Accessor::value_type value, NBL_REF_ARG(Comparator) comp)
 {
-    //using TransformedComparator = impl::lower_to_upper_comparator_transform_t<Accessor,Comparator>;
-    //TransformedComparator transformedComparator;
-
-    impl::lower_to_upper_comparator_transform_t<Accessor,Comparator> transformedComparator;
-    transformedComparator.comp = comp;
-    return lower_bound<Accessor,impl::lower_to_upper_comparator_transform_t<Accessor,Comparator> >(accessor,begin,end,value,transformedComparator);
+    impl::bound_t<Accessor,Comparator,true> implementation = impl::bound_t<Accessor,Comparator,true>::setup(begin,end,value,comp);
+    const uint32_t retval = implementation(accessor);
+    comp = implementation.comp;
+    return retval;
 }
 
 

include/nbl/builtin/hlsl/bxdf/base/lambertian.hlsl

@@ -37,16 +37,18 @@ struct SLambertianBase
     template<typename C=bool_constant<!IsBSDF> >
     enable_if_t<C::value && !IsBSDF, sample_type> generate(NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction, const vector2_type u) NBL_CONST_MEMBER_FUNC
     {
+        typename sampling::ProjectedHemisphere<scalar_type>::cache_type cache;
         ray_dir_info_type L;
-        L.setDirection(sampling::ProjectedHemisphere<scalar_type>::generate(u));
+        L.setDirection(sampling::ProjectedHemisphere<scalar_type>::generate(u, cache));
         return sample_type::createFromTangentSpace(L, interaction.getFromTangentSpace());
     }
     template<typename C=bool_constant<IsBSDF> >
     enable_if_t<C::value && IsBSDF, sample_type> generate(NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction, const vector3_type u) NBL_CONST_MEMBER_FUNC
     {
+        typename sampling::ProjectedSphere<scalar_type>::cache_type cache;
         vector3_type _u = u;
         ray_dir_info_type L;
-        L.setDirection(sampling::ProjectedSphere<scalar_type>::generate(_u));
+        L.setDirection(sampling::ProjectedSphere<scalar_type>::generate(_u, cache));
         return sample_type::createFromTangentSpace(L, interaction.getFromTangentSpace());
     }
     template<typename C=bool_constant<!IsBSDF> >
@@ -63,9 +65,9 @@ struct SLambertianBase
     scalar_type pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(isotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC
     {
         NBL_IF_CONSTEXPR (IsBSDF)
-            return sampling::ProjectedSphere<scalar_type>::pdf(_sample.getNdotL(_clamp));
+            return sampling::ProjectedSphere<scalar_type>::backwardPdf(vector<scalar_type, 3>(0, 0, _sample.getNdotL(_clamp)));
         else
-            return sampling::ProjectedHemisphere<scalar_type>::pdf(_sample.getNdotL(_clamp));
+            return sampling::ProjectedHemisphere<scalar_type>::backwardPdf(vector<scalar_type, 3>(0, 0, _sample.getNdotL(_clamp)));
     }
     scalar_type pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC
     {
@@ -74,12 +76,15 @@ struct SLambertianBase
 
     quotient_pdf_type quotient_and_pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(isotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC
     {
-        sampling::quotient_and_pdf<monochrome_type, scalar_type> qp;
+        // only z component matters: ProjectedHemisphere::forwardPdf reads v.z
+        const vector3_type L = vector3_type(0, 0, _sample.getNdotL(_clamp));
+        typename sampling::ProjectedHemisphere<scalar_type>::cache_type cache;
+        scalar_type p;
         NBL_IF_CONSTEXPR (IsBSDF)
-            qp = sampling::ProjectedSphere<scalar_type>::template quotient_and_pdf(_sample.getNdotL(_clamp));
+            p = sampling::ProjectedSphere<scalar_type>::forwardPdf(L, cache);
         else
-            qp = sampling::ProjectedHemisphere<scalar_type>::template quotient_and_pdf(_sample.getNdotL(_clamp));
-        return quotient_pdf_type::create(qp.quotient[0], qp.pdf);
+            p = sampling::ProjectedHemisphere<scalar_type>::forwardPdf(L, cache);
+        return quotient_pdf_type::create(scalar_type(1.0), p);
     }
     quotient_pdf_type quotient_and_pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC
     {

include/nbl/builtin/hlsl/bxdf/base/oren_nayar.hlsl

@@ -72,16 +72,18 @@ struct SOrenNayarBase
     template<typename C=bool_constant<!IsBSDF> >
     enable_if_t<C::value && !IsBSDF, sample_type> generate(NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction, const vector2_type u) NBL_CONST_MEMBER_FUNC
     {
+        typename sampling::ProjectedHemisphere<scalar_type>::cache_type cache;
         ray_dir_info_type L;
-        L.setDirection(sampling::ProjectedHemisphere<scalar_type>::generate(u));
+        L.setDirection(sampling::ProjectedHemisphere<scalar_type>::generate(u, cache));
         return sample_type::createFromTangentSpace(L, interaction.getFromTangentSpace());
     }
     template<typename C=bool_constant<IsBSDF> >
     enable_if_t<C::value && IsBSDF, sample_type> generate(NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction, const vector3_type u) NBL_CONST_MEMBER_FUNC
     {
+        typename sampling::ProjectedSphere<scalar_type>::cache_type cache;
         vector3_type _u = u;
         ray_dir_info_type L;
-        L.setDirection(sampling::ProjectedSphere<scalar_type>::generate(_u));
+        L.setDirection(sampling::ProjectedSphere<scalar_type>::generate(_u, cache));
         return sample_type::createFromTangentSpace(L, interaction.getFromTangentSpace());
     }
     template<typename C=bool_constant<!IsBSDF> >
@@ -98,9 +100,9 @@ struct SOrenNayarBase
     scalar_type pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(isotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC
     {
         if (IsBSDF)
-            return sampling::ProjectedSphere<scalar_type>::pdf(_sample.getNdotL(_clamp));
+            return sampling::ProjectedSphere<scalar_type>::backwardPdf(vector<scalar_type, 3>(0, 0, _sample.getNdotL(_clamp)));
         else
-            return sampling::ProjectedHemisphere<scalar_type>::pdf(_sample.getNdotL(_clamp));
+            return sampling::ProjectedHemisphere<scalar_type>::backwardPdf(vector<scalar_type, 3>(0, 0, _sample.getNdotL(_clamp)));
     }
     scalar_type pdf(NBL_CONST_REF_ARG(sample_type) _sample, NBL_CONST_REF_ARG(anisotropic_interaction_type) interaction) NBL_CONST_MEMBER_FUNC
     {

include/nbl/builtin/hlsl/ieee754.hlsl

@@ -291,6 +291,41 @@ NBL_CONSTEXPR_FUNC bool isZero(T val)
 	return !(ieee754::impl::bitCastToUintType(val) & exponentAndMantissaMask);
 }
 
+// Returns the largest representable value less than `val`.
+// For positive values this decrements the bit representation; for negative values it increments.
+// Caller must guarantee val is finite and non-zero.
+template <bool CanBeNeg = true, typename T NBL_FUNC_REQUIRES(hlsl::is_floating_point_v<T>)
+NBL_CONSTEXPR_FUNC T nextDown(T val)
+{
+	using AsUint = typename unsigned_integer_of_size<sizeof(T)>::type;
+	using traits_t = traits<T>;
+
+	const AsUint bits = ieee754::impl::bitCastToUintType(val);
+
+	// positive: decrement; negative: increment
+	AsUint result;
+	if (CanBeNeg)
+	{
+		const bool isNegative = (bits & traits_t::signMask) != AsUint(0);
+		result = isNegative ? (bits + AsUint(1)) : (bits - AsUint(1));
+	}
+	else
+		result = bits - AsUint(1);
+	return impl::castBackToFloatType<T>(result);
+}
+
+// Returns the representable value nearest to `val` in the direction of zero.
+// For positive values this decrements the bit representation; for negative values it decrements (moving toward zero).
+// Caller must guarantee val is finite and non-zero.
+template <typename T NBL_FUNC_REQUIRES(hlsl::is_floating_point_v<T>)
+NBL_CONSTEXPR_FUNC T nextTowardZero(T val)
+{
+	using AsUint = typename unsigned_integer_of_size<sizeof(T)>::type;
+
+	const AsUint bits = ieee754::impl::bitCastToUintType(val);
+	return impl::castBackToFloatType<T>(bits - AsUint(1));
+}
+
 }
 }
 }

include/nbl/builtin/hlsl/math/angle_adding.hlsl

@@ -16,10 +16,10 @@ namespace hlsl
 namespace math
 {
 
-template <typename T>
+template <typename T, template<typename NBL_STRUCT_CONSTRAINABLE> class AcosHelper = tgmath_impl::acos_helper>
 struct sincos_accumulator
 {
-    using this_t = sincos_accumulator<T>;
+    using this_t = sincos_accumulator<T, AcosHelper>;
 
     static this_t create(T cosA)
     {
@@ -66,10 +66,29 @@ struct sincos_accumulator
         addAngle(cosA, sqrt<T>(T(1.0) - cosA * cosA));
     }
 
-    // TODO: rename to `getSumOfArccos`
-    T getSumofArccos()
+    T getSumOfArccos()
     {
-        return acos<T>(runningSum.real()) + wraparound;
+        return AcosHelper<T>::__call(runningSum.real()) + wraparound;
+    }
+
+    // Returns (sum of angles) - pi using the identity acos(x) - pi = -acos(-x).
+    // Avoids catastrophic cancellation when wraparound == 0 and sum ~ pi,
+    // where the naive `acos(~-1) - pi` subtracts two values both near pi.
+    // Note: when wraparound == pi (typical for spherical triangles with small
+    // excess e), the naive path `acos(cos(e))` is actually slightly more precise
+    // since it returns a small value directly without a near-pi cancellation.
+    T getSumOfArccosMinusPi()
+    {
+        return -AcosHelper<T>::__call(-runningSum.real()) + wraparound;
+    }
+
+    // Same as getSumOfArccosMinusPi but clamps the accumulated cosine to [-1,1]
+    // before calling acos. The angle addition in addAngle() uses sum-of-products
+    // which can push the accumulated cosine slightly outside the unit circle due
+    // to FP rounding, causing acos to return NaN on GPU where max(NaN,0)=0.
+    T getClampedSumOfArccosMinusPi()
+    {
+        return -AcosHelper<T>::__call(clamp<T>(-runningSum.real(), T(-1), T(1))) + wraparound;
     }
 
     complex_t<T> runningSum;

include/nbl/builtin/hlsl/math/fast_acos.hlsl

@@ -0,0 +1,87 @@
+// Copyright (C) 2018-2023 - DevSH Graphics Programming Sp. z O.O.
+// This file is part of the "Nabla Engine".
+// For conditions of distribution and use, see copyright notice in nabla.h
+#ifndef _NBL_BUILTIN_HLSL_MATH_FAST_ACOS_INCLUDED_
+#define _NBL_BUILTIN_HLSL_MATH_FAST_ACOS_INCLUDED_
+
+#include "nbl/builtin/hlsl/cpp_compat.hlsl"
+#include "nbl/builtin/hlsl/numbers.hlsl"
+#include "nbl/builtin/hlsl/tgmath.hlsl"
+#include "nbl/builtin/hlsl/concepts.hlsl"
+
+namespace nbl
+{
+namespace hlsl
+{
+namespace math
+{
+// https://www.desmos.com/calculator/a59pbwgwof
+// a comparison between fast acos methods
+
+// Attempt 1: Polynomial approximation of acos(x) for x in [-1,1]
+// Based on an odd cubic fit: acos(x) ~ (a * x^2 + b) * x + pi/2
+// Max absolute error: ~9.8e-2 rad
+// Very fast: no sqrt, no branches
+template<typename T NBL_STRUCT_CONSTRAINABLE>
+struct fast_acos_cubic
+{
+	static T __call(const T val)
+	{
+		return (T(-0.621565443625098) * val * val - T(0.837561827808302)) * val + numbers::pi<T> / T(2.0);
+	}
+};
+
+// Attempt 2: Abramowitz & Stegun formula 4.4.45 (adapted for full [-1,1] range)
+// acos(x) ~ sqrt(1-x) * (a0 + a1*x + a2*x^2 + a3*x^3) for x in [0,1]
+// For x in [-1,0]: acos(x) = pi - acos(-x)
+// Max absolute error: ~6.3e-5 rad
+template<typename T NBL_STRUCT_CONSTRAINABLE>
+struct fast_acos_stegun_poly3
+{
+	static T __call(const T val)
+	{
+		const T ax = abs<T>(val);
+		const T s = sqrt<T>(T(1.0) - ax);
+		const T poly = ((T(-0.019771840941) * ax + T(0.075701735421)) * ax - T(0.212644584569)) * ax + T(1.570771931669); // not pi/2, free constant for better fit
+		const T result = s * poly;
+		return hlsl::select(val < T(0.0), numbers::pi<T> - result, result);
+	}
+};
+
+// Attempt 3: Degree-4 polynomial, good accuracy/cost tradeoff
+// Fitted as: acos(x) ~ sqrt(1-x) * (a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4) for x in [0,1]
+// Max absolute error: ~8.6e-6 rad
+template<typename T NBL_STRUCT_CONSTRAINABLE>
+struct fast_acos_stegun_poly4
+{
+	static T __call(const T val)
+	{
+		const T ax = abs<T>(val);
+		const T s = sqrt<T>(T(1.0) - ax);
+		const T poly = (((T(0.0102831457) * ax - T(0.0387865708)) * ax + T(0.0864014439)) * ax - T(0.2144371342)) * ax + numbers::pi<T> / T(2.0);
+		const T result = s * poly;
+		return hlsl::select(val < T(0.0), numbers::pi<T> - result, result);
+	}
+};
+
+// Attempt 4: Degree-5 polynomial for best accuracy
+// Fitted as: acos(x) ~ sqrt(1-x) * (a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5) for x in [0,1]
+// Max absolute error: ~1.1e-6 rad
+template<typename T NBL_STRUCT_CONSTRAINABLE>
+struct fast_acos_stegun_poly5
+{
+	static T __call(const T val)
+	{
+		const T ax = abs<T>(val);
+		const T s = sqrt<T>(T(1.0) - ax);
+		const T poly = ((((T(-0.0051108043) * ax + T(0.0211860706)) * ax - T(0.0464796661)) * ax + T(0.0883884911)) * ax - T(0.2145725267)) * ax + numbers::pi<T> / T(2.0);
+		const T result = s * poly;
+		return hlsl::select(val < T(0.0), numbers::pi<T> - result, result);
+	}
+};
+
+}
+}
+}
+
+#endif

include/nbl/builtin/hlsl/math/functions.hlsl

@@ -96,9 +96,10 @@ scalar_type_t<T> lpNorm(NBL_CONST_REF_ARG(T) v)
 template <typename T NBL_FUNC_REQUIRES(concepts::FloatingPointLikeScalar<T>)
 void sincos(T theta, NBL_REF_ARG(T) s, NBL_REF_ARG(T) c)
 {
+    s = sin<T>(theta);
     c = cos<T>(theta);
-    s = sqrt<T>(T(NBL_FP64_LITERAL(1.0))-c*c);
-    s = ieee754::flipSign(s, theta < T(NBL_FP64_LITERAL(0.0)));
+    // s = sqrt<T>(T(NBL_FP64_LITERAL(1.0))-c*c);
+    // s = ieee754::flipSign(s, theta < T(NBL_FP64_LITERAL(0.0)));
 }
 
 template <typename T NBL_FUNC_REQUIRES(vector_traits<T>::Dimension == 3)
@@ -136,7 +137,7 @@ struct conditionalAbsOrMax_helper<T NBL_PARTIAL_REQ_BOT(concepts::FloatingPointL
         using UintOfTSize = unsigned_integer_of_size_t<sizeof(T)>;
         const T condAbs = bit_cast<T>(bit_cast<UintOfTSize>(x) & (cond ? (numeric_limits<UintOfTSize>::max >> 1) : numeric_limits<UintOfTSize>::max));
 
-        return max<T>(condAbs, limit);
+        return nbl::hlsl::max<T>(condAbs, limit);
     }
 };
 
@@ -156,7 +157,7 @@ struct conditionalAbsOrMax_helper<T NBL_PARTIAL_REQ_BOT(concepts::FloatingPointL
         const Uint32VectorWithDimensionOfT condAbsAsUint = xAsUintVec & mask;
         T condAbs = bit_cast<T, Uint32VectorWithDimensionOfT>(condAbsAsUint);
 
-        return max<T>(condAbs, limit);
+        return nbl::hlsl::max<T>(condAbs, limit);
     }
 };
 }