Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorBrecht Van Lommel <brecht@blender.org>2021-10-24 15:19:19 +0300
committerBrecht Van Lommel <brecht@blender.org>2021-10-26 16:37:04 +0300
commitfd25e883e2807a151f673b87c152a59701a0df80 (patch)
tree9441933f32ba2672ca71c58842342a9c525e123e /intern/cycles/util/math.h
parentd7d40745fa09061a3117bd3669c5a46bbf611eae (diff)
Cycles: remove prefix from source code file names
Remove prefix of filenames that is the same as the folder name. This used to help when #includes were using individual files, but now they are always relative to the cycles root directory and so the prefixes are redundant. For patches and branches, git merge and rebase should be able to detect the renames and move over code to the right file.
Diffstat (limited to 'intern/cycles/util/math.h')
-rw-r--r--intern/cycles/util/math.h870
1 files changed, 870 insertions, 0 deletions
diff --git a/intern/cycles/util/math.h b/intern/cycles/util/math.h
new file mode 100644
index 00000000000..e7fc492733f
--- /dev/null
+++ b/intern/cycles/util/math.h
@@ -0,0 +1,870 @@
+/*
+ * Copyright 2011-2013 Blender Foundation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef __UTIL_MATH_H__
+#define __UTIL_MATH_H__
+
+/* Math
+ *
+ * Basic math functions on scalar and vector types. This header is used by
+ * both the kernel code when compiled as C++, and other C++ non-kernel code. */
+
+#ifndef __KERNEL_GPU__
+# include <cmath>
+#endif
+
+#ifdef __HIP__
+# include <hip/hip_vector_types.h>
+#endif
+
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "util/types.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* Float Pi variations */
+
+/* Division */
+#ifndef M_PI_F
+# define M_PI_F (3.1415926535897932f) /* pi */
+#endif
+#ifndef M_PI_2_F
+# define M_PI_2_F (1.5707963267948966f) /* pi/2 */
+#endif
+#ifndef M_PI_4_F
+# define M_PI_4_F (0.7853981633974830f) /* pi/4 */
+#endif
+#ifndef M_1_PI_F
+# define M_1_PI_F (0.3183098861837067f) /* 1/pi */
+#endif
+#ifndef M_2_PI_F
+# define M_2_PI_F (0.6366197723675813f) /* 2/pi */
+#endif
+#ifndef M_1_2PI_F
+# define M_1_2PI_F (0.1591549430918953f) /* 1/(2*pi) */
+#endif
+#ifndef M_SQRT_PI_8_F
+# define M_SQRT_PI_8_F (0.6266570686577501f) /* sqrt(pi/8) */
+#endif
+#ifndef M_LN_2PI_F
+# define M_LN_2PI_F (1.8378770664093454f) /* ln(2*pi) */
+#endif
+
+/* Multiplication */
+#ifndef M_2PI_F
+# define M_2PI_F (6.2831853071795864f) /* 2*pi */
+#endif
+#ifndef M_4PI_F
+# define M_4PI_F (12.566370614359172f) /* 4*pi */
+#endif
+
+/* Float sqrt variations */
+#ifndef M_SQRT2_F
+# define M_SQRT2_F (1.4142135623730950f) /* sqrt(2) */
+#endif
+#ifndef M_LN2_F
+# define M_LN2_F (0.6931471805599453f) /* ln(2) */
+#endif
+#ifndef M_LN10_F
+# define M_LN10_F (2.3025850929940457f) /* ln(10) */
+#endif
+
+/* Scalar */
+
+#ifndef __HIP__
+# ifdef _WIN32
+ccl_device_inline float fmaxf(float a, float b)
+{
+ return (a > b) ? a : b;
+}
+
+ccl_device_inline float fminf(float a, float b)
+{
+ return (a < b) ? a : b;
+}
+
+# endif /* _WIN32 */
+#endif /* __HIP__ */
+
+#ifndef __KERNEL_GPU__
+using std::isfinite;
+using std::isnan;
+using std::sqrt;
+
+ccl_device_inline int abs(int x)
+{
+ return (x > 0) ? x : -x;
+}
+
+ccl_device_inline int max(int a, int b)
+{
+ return (a > b) ? a : b;
+}
+
+ccl_device_inline int min(int a, int b)
+{
+ return (a < b) ? a : b;
+}
+
+ccl_device_inline uint min(uint a, uint b)
+{
+ return (a < b) ? a : b;
+}
+
+ccl_device_inline float max(float a, float b)
+{
+ return (a > b) ? a : b;
+}
+
+ccl_device_inline float min(float a, float b)
+{
+ return (a < b) ? a : b;
+}
+
+ccl_device_inline double max(double a, double b)
+{
+ return (a > b) ? a : b;
+}
+
+ccl_device_inline double min(double a, double b)
+{
+ return (a < b) ? a : b;
+}
+
+/* These 2 guys are templated for usage with registers data.
+ *
+ * NOTE: Since this is CPU-only functions it is ok to use references here.
+ * But for other devices we'll need to be careful about this.
+ */
+
+template<typename T> ccl_device_inline T min4(const T &a, const T &b, const T &c, const T &d)
+{
+ return min(min(a, b), min(c, d));
+}
+
+template<typename T> ccl_device_inline T max4(const T &a, const T &b, const T &c, const T &d)
+{
+ return max(max(a, b), max(c, d));
+}
+#endif /* __KERNEL_GPU__ */
+
+ccl_device_inline float min4(float a, float b, float c, float d)
+{
+ return min(min(a, b), min(c, d));
+}
+
+ccl_device_inline float max4(float a, float b, float c, float d)
+{
+ return max(max(a, b), max(c, d));
+}
+
+/* Int/Float conversion */
+
+ccl_device_inline int as_int(uint i)
+{
+ union {
+ uint ui;
+ int i;
+ } u;
+ u.ui = i;
+ return u.i;
+}
+
+ccl_device_inline uint as_uint(int i)
+{
+ union {
+ uint ui;
+ int i;
+ } u;
+ u.i = i;
+ return u.ui;
+}
+
+ccl_device_inline uint as_uint(float f)
+{
+ union {
+ uint i;
+ float f;
+ } u;
+ u.f = f;
+ return u.i;
+}
+
+#ifndef __HIP__
+ccl_device_inline int __float_as_int(float f)
+{
+ union {
+ int i;
+ float f;
+ } u;
+ u.f = f;
+ return u.i;
+}
+
+ccl_device_inline float __int_as_float(int i)
+{
+ union {
+ int i;
+ float f;
+ } u;
+ u.i = i;
+ return u.f;
+}
+
+ccl_device_inline uint __float_as_uint(float f)
+{
+ union {
+ uint i;
+ float f;
+ } u;
+ u.f = f;
+ return u.i;
+}
+
+ccl_device_inline float __uint_as_float(uint i)
+{
+ union {
+ uint i;
+ float f;
+ } u;
+ u.i = i;
+ return u.f;
+}
+#endif
+
+ccl_device_inline int4 __float4_as_int4(float4 f)
+{
+#ifdef __KERNEL_SSE__
+ return int4(_mm_castps_si128(f.m128));
+#else
+ return make_int4(
+ __float_as_int(f.x), __float_as_int(f.y), __float_as_int(f.z), __float_as_int(f.w));
+#endif
+}
+
+ccl_device_inline float4 __int4_as_float4(int4 i)
+{
+#ifdef __KERNEL_SSE__
+ return float4(_mm_castsi128_ps(i.m128));
+#else
+ return make_float4(
+ __int_as_float(i.x), __int_as_float(i.y), __int_as_float(i.z), __int_as_float(i.w));
+#endif
+}
+
+template<typename T> ccl_device_inline uint pointer_pack_to_uint_0(T *ptr)
+{
+ return ((uint64_t)ptr) & 0xFFFFFFFF;
+}
+
+template<typename T> ccl_device_inline uint pointer_pack_to_uint_1(T *ptr)
+{
+ return (((uint64_t)ptr) >> 32) & 0xFFFFFFFF;
+}
+
+template<typename T> ccl_device_inline T *pointer_unpack_from_uint(const uint a, const uint b)
+{
+ return (T *)(((uint64_t)b << 32) | a);
+}
+
+ccl_device_inline uint uint16_pack_to_uint(const uint a, const uint b)
+{
+ return (a << 16) | b;
+}
+
+ccl_device_inline uint uint16_unpack_from_uint_0(const uint i)
+{
+ return i >> 16;
+}
+
+ccl_device_inline uint uint16_unpack_from_uint_1(const uint i)
+{
+ return i & 0xFFFF;
+}
+
+/* Versions of functions which are safe for fast math. */
+ccl_device_inline bool isnan_safe(float f)
+{
+ unsigned int x = __float_as_uint(f);
+ return (x << 1) > 0xff000000u;
+}
+
+ccl_device_inline bool isfinite_safe(float f)
+{
+ /* By IEEE 754 rule, 2*Inf equals Inf */
+ unsigned int x = __float_as_uint(f);
+ return (f == f) && (x == 0 || x == (1u << 31) || (f != 2.0f * f)) && !((x << 1) > 0xff000000u);
+}
+
+ccl_device_inline float ensure_finite(float v)
+{
+ return isfinite_safe(v) ? v : 0.0f;
+}
+
+ccl_device_inline int clamp(int a, int mn, int mx)
+{
+ return min(max(a, mn), mx);
+}
+
+ccl_device_inline float clamp(float a, float mn, float mx)
+{
+ return min(max(a, mn), mx);
+}
+
+ccl_device_inline float mix(float a, float b, float t)
+{
+ return a + t * (b - a);
+}
+
+ccl_device_inline float smoothstep(float edge0, float edge1, float x)
+{
+ float result;
+ if (x < edge0)
+ result = 0.0f;
+ else if (x >= edge1)
+ result = 1.0f;
+ else {
+ float t = (x - edge0) / (edge1 - edge0);
+ result = (3.0f - 2.0f * t) * (t * t);
+ }
+ return result;
+}
+
+#ifndef __KERNEL_CUDA__
+ccl_device_inline float saturate(float a)
+{
+ return clamp(a, 0.0f, 1.0f);
+}
+#endif /* __KERNEL_CUDA__ */
+
+ccl_device_inline int float_to_int(float f)
+{
+ return (int)f;
+}
+
+ccl_device_inline int floor_to_int(float f)
+{
+ return float_to_int(floorf(f));
+}
+
+ccl_device_inline int quick_floor_to_int(float x)
+{
+ return float_to_int(x) - ((x < 0) ? 1 : 0);
+}
+
+ccl_device_inline float floorfrac(float x, ccl_private int *i)
+{
+ *i = quick_floor_to_int(x);
+ return x - *i;
+}
+
+ccl_device_inline int ceil_to_int(float f)
+{
+ return float_to_int(ceilf(f));
+}
+
+ccl_device_inline float fractf(float x)
+{
+ return x - floorf(x);
+}
+
+/* Adapted from godot-engine math_funcs.h. */
+ccl_device_inline float wrapf(float value, float max, float min)
+{
+ float range = max - min;
+ return (range != 0.0f) ? value - (range * floorf((value - min) / range)) : min;
+}
+
+ccl_device_inline float pingpongf(float a, float b)
+{
+ return (b != 0.0f) ? fabsf(fractf((a - b) / (b * 2.0f)) * b * 2.0f - b) : 0.0f;
+}
+
+ccl_device_inline float smoothminf(float a, float b, float k)
+{
+ if (k != 0.0f) {
+ float h = fmaxf(k - fabsf(a - b), 0.0f) / k;
+ return fminf(a, b) - h * h * h * k * (1.0f / 6.0f);
+ }
+ else {
+ return fminf(a, b);
+ }
+}
+
+ccl_device_inline float signf(float f)
+{
+ return (f < 0.0f) ? -1.0f : 1.0f;
+}
+
+ccl_device_inline float nonzerof(float f, float eps)
+{
+ if (fabsf(f) < eps)
+ return signf(f) * eps;
+ else
+ return f;
+}
+
+/* `signum` function testing for zero. Matches GLSL and OSL functions. */
+ccl_device_inline float compatible_signf(float f)
+{
+ if (f == 0.0f) {
+ return 0.0f;
+ }
+ else {
+ return signf(f);
+ }
+}
+
+ccl_device_inline float smoothstepf(float f)
+{
+ float ff = f * f;
+ return (3.0f * ff - 2.0f * ff * f);
+}
+
+ccl_device_inline int mod(int x, int m)
+{
+ return (x % m + m) % m;
+}
+
+ccl_device_inline float3 float2_to_float3(const float2 a)
+{
+ return make_float3(a.x, a.y, 0.0f);
+}
+
+ccl_device_inline float3 float4_to_float3(const float4 a)
+{
+ return make_float3(a.x, a.y, a.z);
+}
+
+ccl_device_inline float4 float3_to_float4(const float3 a)
+{
+ return make_float4(a.x, a.y, a.z, 1.0f);
+}
+
+ccl_device_inline float inverse_lerp(float a, float b, float x)
+{
+ return (x - a) / (b - a);
+}
+
+/* Cubic interpolation between b and c, a and d are the previous and next point. */
+ccl_device_inline float cubic_interp(float a, float b, float c, float d, float x)
+{
+ return 0.5f *
+ (((d + 3.0f * (b - c) - a) * x + (2.0f * a - 5.0f * b + 4.0f * c - d)) * x +
+ (c - a)) *
+ x +
+ b;
+}
+
+CCL_NAMESPACE_END
+
+#include "util/math_int2.h"
+#include "util/math_int3.h"
+#include "util/math_int4.h"
+
+#include "util/math_float2.h"
+#include "util/math_float3.h"
+#include "util/math_float4.h"
+
+#include "util/rect.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* Interpolation */
+
+template<class A, class B> A lerp(const A &a, const A &b, const B &t)
+{
+ return (A)(a * ((B)1 - t) + b * t);
+}
+
+/* Triangle */
+
+ccl_device_inline float triangle_area(ccl_private const float3 &v1,
+ ccl_private const float3 &v2,
+ ccl_private const float3 &v3)
+{
+ return len(cross(v3 - v2, v1 - v2)) * 0.5f;
+}
+
+/* Orthonormal vectors */
+
+ccl_device_inline void make_orthonormals(const float3 N,
+ ccl_private float3 *a,
+ ccl_private float3 *b)
+{
+#if 0
+ if (fabsf(N.y) >= 0.999f) {
+ *a = make_float3(1, 0, 0);
+ *b = make_float3(0, 0, 1);
+ return;
+ }
+ if (fabsf(N.z) >= 0.999f) {
+ *a = make_float3(1, 0, 0);
+ *b = make_float3(0, 1, 0);
+ return;
+ }
+#endif
+
+ if (N.x != N.y || N.x != N.z)
+ *a = make_float3(N.z - N.y, N.x - N.z, N.y - N.x); //(1,1,1)x N
+ else
+ *a = make_float3(N.z - N.y, N.x + N.z, -N.y - N.x); //(-1,1,1)x N
+
+ *a = normalize(*a);
+ *b = cross(N, *a);
+}
+
+/* Color division */
+
+ccl_device_inline float3 safe_invert_color(float3 a)
+{
+ float x, y, z;
+
+ x = (a.x != 0.0f) ? 1.0f / a.x : 0.0f;
+ y = (a.y != 0.0f) ? 1.0f / a.y : 0.0f;
+ z = (a.z != 0.0f) ? 1.0f / a.z : 0.0f;
+
+ return make_float3(x, y, z);
+}
+
+ccl_device_inline float3 safe_divide_color(float3 a, float3 b)
+{
+ float x, y, z;
+
+ x = (b.x != 0.0f) ? a.x / b.x : 0.0f;
+ y = (b.y != 0.0f) ? a.y / b.y : 0.0f;
+ z = (b.z != 0.0f) ? a.z / b.z : 0.0f;
+
+ return make_float3(x, y, z);
+}
+
+ccl_device_inline float3 safe_divide_even_color(float3 a, float3 b)
+{
+ float x, y, z;
+
+ x = (b.x != 0.0f) ? a.x / b.x : 0.0f;
+ y = (b.y != 0.0f) ? a.y / b.y : 0.0f;
+ z = (b.z != 0.0f) ? a.z / b.z : 0.0f;
+
+ /* try to get gray even if b is zero */
+ if (b.x == 0.0f) {
+ if (b.y == 0.0f) {
+ x = z;
+ y = z;
+ }
+ else if (b.z == 0.0f) {
+ x = y;
+ z = y;
+ }
+ else
+ x = 0.5f * (y + z);
+ }
+ else if (b.y == 0.0f) {
+ if (b.z == 0.0f) {
+ y = x;
+ z = x;
+ }
+ else
+ y = 0.5f * (x + z);
+ }
+ else if (b.z == 0.0f) {
+ z = 0.5f * (x + y);
+ }
+
+ return make_float3(x, y, z);
+}
+
+/* Rotation of point around axis and angle */
+
+ccl_device_inline float3 rotate_around_axis(float3 p, float3 axis, float angle)
+{
+ float costheta = cosf(angle);
+ float sintheta = sinf(angle);
+ float3 r;
+
+ r.x = ((costheta + (1 - costheta) * axis.x * axis.x) * p.x) +
+ (((1 - costheta) * axis.x * axis.y - axis.z * sintheta) * p.y) +
+ (((1 - costheta) * axis.x * axis.z + axis.y * sintheta) * p.z);
+
+ r.y = (((1 - costheta) * axis.x * axis.y + axis.z * sintheta) * p.x) +
+ ((costheta + (1 - costheta) * axis.y * axis.y) * p.y) +
+ (((1 - costheta) * axis.y * axis.z - axis.x * sintheta) * p.z);
+
+ r.z = (((1 - costheta) * axis.x * axis.z - axis.y * sintheta) * p.x) +
+ (((1 - costheta) * axis.y * axis.z + axis.x * sintheta) * p.y) +
+ ((costheta + (1 - costheta) * axis.z * axis.z) * p.z);
+
+ return r;
+}
+
+/* NaN-safe math ops */
+
+ccl_device_inline float safe_sqrtf(float f)
+{
+ return sqrtf(max(f, 0.0f));
+}
+
+ccl_device_inline float inversesqrtf(float f)
+{
+ return (f > 0.0f) ? 1.0f / sqrtf(f) : 0.0f;
+}
+
+ccl_device float safe_asinf(float a)
+{
+ return asinf(clamp(a, -1.0f, 1.0f));
+}
+
+ccl_device float safe_acosf(float a)
+{
+ return acosf(clamp(a, -1.0f, 1.0f));
+}
+
+ccl_device float compatible_powf(float x, float y)
+{
+#ifdef __KERNEL_GPU__
+ if (y == 0.0f) /* x^0 -> 1, including 0^0 */
+ return 1.0f;
+
+ /* GPU pow doesn't accept negative x, do manual checks here */
+ if (x < 0.0f) {
+ if (fmodf(-y, 2.0f) == 0.0f)
+ return powf(-x, y);
+ else
+ return -powf(-x, y);
+ }
+ else if (x == 0.0f)
+ return 0.0f;
+#endif
+ return powf(x, y);
+}
+
+ccl_device float safe_powf(float a, float b)
+{
+ if (UNLIKELY(a < 0.0f && b != float_to_int(b)))
+ return 0.0f;
+
+ return compatible_powf(a, b);
+}
+
+ccl_device float safe_divide(float a, float b)
+{
+ return (b != 0.0f) ? a / b : 0.0f;
+}
+
+ccl_device float safe_logf(float a, float b)
+{
+ if (UNLIKELY(a <= 0.0f || b <= 0.0f))
+ return 0.0f;
+
+ return safe_divide(logf(a), logf(b));
+}
+
+ccl_device float safe_modulo(float a, float b)
+{
+ return (b != 0.0f) ? fmodf(a, b) : 0.0f;
+}
+
+ccl_device_inline float sqr(float a)
+{
+ return a * a;
+}
+
+ccl_device_inline float pow20(float a)
+{
+ return sqr(sqr(sqr(sqr(a)) * a));
+}
+
+ccl_device_inline float pow22(float a)
+{
+ return sqr(a * sqr(sqr(sqr(a)) * a));
+}
+
+ccl_device_inline float beta(float x, float y)
+{
+ return expf(lgammaf(x) + lgammaf(y) - lgammaf(x + y));
+}
+
+ccl_device_inline float xor_signmask(float x, int y)
+{
+ return __int_as_float(__float_as_int(x) ^ y);
+}
+
+ccl_device float bits_to_01(uint bits)
+{
+ return bits * (1.0f / (float)0xFFFFFFFF);
+}
+
+ccl_device_inline uint count_leading_zeros(uint x)
+{
+#if defined(__KERNEL_CUDA__) || defined(__KERNEL_OPTIX__) || defined(__KERNEL_HIP__)
+ return __clz(x);
+#else
+ assert(x != 0);
+# ifdef _MSC_VER
+ unsigned long leading_zero = 0;
+ _BitScanReverse(&leading_zero, x);
+ return (31 - leading_zero);
+# else
+ return __builtin_clz(x);
+# endif
+#endif
+}
+
+ccl_device_inline uint count_trailing_zeros(uint x)
+{
+#if defined(__KERNEL_CUDA__) || defined(__KERNEL_OPTIX__) || defined(__KERNEL_HIP__)
+ return (__ffs(x) - 1);
+#else
+ assert(x != 0);
+# ifdef _MSC_VER
+ unsigned long ctz = 0;
+ _BitScanForward(&ctz, x);
+ return ctz;
+# else
+ return __builtin_ctz(x);
+# endif
+#endif
+}
+
+ccl_device_inline uint find_first_set(uint x)
+{
+#if defined(__KERNEL_CUDA__) || defined(__KERNEL_OPTIX__) || defined(__KERNEL_HIP__)
+ return __ffs(x);
+#else
+# ifdef _MSC_VER
+ return (x != 0) ? (32 - count_leading_zeros(x & (-x))) : 0;
+# else
+ return __builtin_ffs(x);
+# endif
+#endif
+}
+
+/* projections */
+ccl_device_inline float2 map_to_tube(const float3 co)
+{
+ float len, u, v;
+ len = sqrtf(co.x * co.x + co.y * co.y);
+ if (len > 0.0f) {
+ u = (1.0f - (atan2f(co.x / len, co.y / len) / M_PI_F)) * 0.5f;
+ v = (co.z + 1.0f) * 0.5f;
+ }
+ else {
+ u = v = 0.0f;
+ }
+ return make_float2(u, v);
+}
+
+ccl_device_inline float2 map_to_sphere(const float3 co)
+{
+ float l = len(co);
+ float u, v;
+ if (l > 0.0f) {
+ if (UNLIKELY(co.x == 0.0f && co.y == 0.0f)) {
+ u = 0.0f; /* Otherwise domain error. */
+ }
+ else {
+ u = (1.0f - atan2f(co.x, co.y) / M_PI_F) / 2.0f;
+ }
+ v = 1.0f - safe_acosf(co.z / l) / M_PI_F;
+ }
+ else {
+ u = v = 0.0f;
+ }
+ return make_float2(u, v);
+}
+
+/* Compares two floats.
+ * Returns true if their absolute difference is smaller than abs_diff (for numbers near zero)
+ * or their relative difference is less than ulp_diff ULPs.
+ * Based on
+ * https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+ */
+
+ccl_device_inline float compare_floats(float a, float b, float abs_diff, int ulp_diff)
+{
+ if (fabsf(a - b) < abs_diff) {
+ return true;
+ }
+
+ if ((a < 0.0f) != (b < 0.0f)) {
+ return false;
+ }
+
+ return (abs(__float_as_int(a) - __float_as_int(b)) < ulp_diff);
+}
+
+/* Calculate the angle between the two vectors a and b.
+ * The usual approach `acos(dot(a, b))` has severe precision issues for small angles,
+ * which are avoided by this method.
+ * Based on "Mangled Angles" from https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf
+ */
+ccl_device_inline float precise_angle(float3 a, float3 b)
+{
+ return 2.0f * atan2f(len(a - b), len(a + b));
+}
+
+/* Return value which is greater than the given one and is a power of two. */
+ccl_device_inline uint next_power_of_two(uint x)
+{
+ return x == 0 ? 1 : 1 << (32 - count_leading_zeros(x));
+}
+
+/* Return value which is lower than the given one and is a power of two. */
+ccl_device_inline uint prev_power_of_two(uint x)
+{
+ return x < 2 ? x : 1 << (31 - count_leading_zeros(x - 1));
+}
+
+#ifndef __has_builtin
+# define __has_builtin(v) 0
+#endif
+
+/* Reverses the bits of a 32 bit integer. */
+ccl_device_inline uint32_t reverse_integer_bits(uint32_t x)
+{
+ /* Use a native instruction if it exists. */
+#if defined(__arm__) || defined(__aarch64__)
+ __asm__("rbit %w0, %w1" : "=r"(x) : "r"(x));
+ return x;
+#elif defined(__KERNEL_CUDA__)
+ return __brev(x);
+#elif __has_builtin(__builtin_bitreverse32)
+ return __builtin_bitreverse32(x);
+#else
+ /* Flip pairwise. */
+ x = ((x & 0x55555555) << 1) | ((x & 0xAAAAAAAA) >> 1);
+ /* Flip pairs. */
+ x = ((x & 0x33333333) << 2) | ((x & 0xCCCCCCCC) >> 2);
+ /* Flip nibbles. */
+ x = ((x & 0x0F0F0F0F) << 4) | ((x & 0xF0F0F0F0) >> 4);
+ /* Flip bytes. CPUs have an instruction for that, pretty fast one. */
+# ifdef _MSC_VER
+ return _byteswap_ulong(x);
+# elif defined(__INTEL_COMPILER)
+ return (uint32_t)_bswap((int)x);
+# else
+ /* Assuming gcc or clang. */
+ return __builtin_bswap32(x);
+# endif
+#endif
+}
+
+CCL_NAMESPACE_END
+
+#endif /* __UTIL_MATH_H__ */