ldexp implementation

Rinzii · Mar 25, 2024 · 6790252 · 6790252
1 parent 536c1e0
commit 6790252
Show file tree

Hide file tree

Showing 2 changed files with 83 additions and 0 deletions.
diff --git a/include/ccmath/internal/support/bits.hpp b/include/ccmath/internal/support/bits.hpp
@@ -56,6 +56,26 @@ namespace ccm::helpers
 		return x && !(x & (x - 1));
 	}
 
+	template <typename T, std::enable_if_t<!std::is_integral_v<T>, bool> = true>
+	inline constexpr std::uint32_t get_exponent_of_floating_point(T x) noexcept
+	{
+		const auto bits				= bit_cast<float_bits_t<T>>(x);
+		const auto shifted_exponent = bits >> ccm::helpers::floating_point_traits<T>::exponent_shift;
+		const auto masked_exponent	= shifted_exponent & ccm::helpers::floating_point_traits<T>::exponent_mask;
+		return masked_exponent - ccm::helpers::floating_point_traits<T>::exponent_bias;
+	}
+
+	template <typename T, std::enable_if_t<!std::is_integral_v<T>, bool> = true>
+	inline constexpr void set_exponent_of_floating_point(T & x, int exp) noexcept
+	{
+		const auto bit_casted			  = bit_cast<float_bits_t<T>>(x);
+		const auto inverted_exponent_mask = ~ccm::helpers::floating_point_traits<T>::shifted_exponent_mask;
+		const auto masked_exponent = (exp + ccm::helpers::floating_point_traits<T>::exponent_bias) & ccm::helpers::floating_point_traits<T>::exponent_mask;
+		const auto shifted_masked_exponent = masked_exponent << ccm::helpers::floating_point_traits<T>::exponent_shift;
+		const auto final_bits			   = (bit_casted & inverted_exponent_mask) | shifted_masked_exponent;
+
+		x = bit_cast<T>(final_bits);
+	}
 
 
 

diff --git a/include/ccmath/math/fmanip/ldexp.hpp b/include/ccmath/math/fmanip/ldexp.hpp
@@ -10,5 +10,68 @@
 
 namespace ccm
 {
+	// this function multiples parameter x by two to the power of paramter powerOf2
+	// x * 2^powerOf2
+	template <typename T, std::enable_if_t<!std::is_integral_v<T>, bool> = true>
+	inline constexpr T ldexp(T x, int32_t powerOf2) noexcept
+	{
+		if (!ccm::isfinite(x)) { return x; }
+		int32_t oldexp = ccm::helpers::get_exponent_of_floating_point<T>(x);
 
+		// if the mantissa is 0 and the original exponent is 0
+		if ((oldexp == 0) && ((ccm::helpers::bit_cast<ccm::helpers::floating_point_traits<T>::uint_type>(x) &
+							   ccm::helpers::floating_point_traits<T>::normal_mantissa_mask) == 0))
+		{
+			return x;
+		}
+
+		if (powerOf2 > ccm::helpers::floating_point_traits<T>::maximum_binary_exponent)
+		{
+			// error == hugeval
+			return std::numeric_limits<T>::infinity();
+		}
+		// the reference source says -2 * exp_max
+		else if (powerOf2 < ccm::helpers::floating_point_traits<T>::minimum_binary_exponent)
+		{
+			// error == range
+			return 0;
+		}
+		T factor;
+		// normalizes an abnormal floating point
+		if (oldexp == 0)
+		{
+			factor = 1 << (sizeof(T) * CHAR_BIT);
+			x *= factor;
+			powerOf2 = -sizeof(T) * CHAR_BIT;
+			oldexp	 = ccm::helpers::get_exponent_of_floating_point<T>(x);
+		}
+
+		powerOf2 = oldexp + powerOf2;
+		if (powerOf2 >= ccm::helpers::floating_point_traits<T>::maximum_binary_exponent)
+		{
+			// overflow
+			return std::numeric_limits<T>::infinity();
+		}
+		if (powerOf2 > 0)
+		{
+			ccm::helpers::set_exponent_of_floating_point<T>(x, powerOf2);
+			return x;
+		}
+		// denormal, or underflow
+		powerOf2 += sizeof(T) * CHAR_BIT;
+		ccm::helpers::set_exponent_of_floating_point<T>(x, powerOf2);
+		factor = 1 << (sizeof(T) * CHAR_BIT);
+		x /= factor;
+		if (x == static_cast<T>(0))
+		{
+			// underflow report
+		}
+		return x;
+	}
+
+	template <typename T, std::enable_if_t<std::is_integral_v<T>, bool> = true>
+	inline constexpr double ldexp(T x, int32_t powerOf2)
+	{
+		return ldexp(static_cast<double>(x), powerOf2);
+	}
 } // namespace ccm