Skip to content

Commit

Permalink
fix main conflicts
Browse files Browse the repository at this point in the history
  • Loading branch information
Rinzii committed Apr 6, 2024
2 parents 41f7f15 + b9dae17 commit 20409c4
Show file tree
Hide file tree
Showing 5 changed files with 316 additions and 2 deletions.
124 changes: 124 additions & 0 deletions include/ccmath/detail/basic/impl/remquo_double_impl.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
/*
* Copyright (c) 2024-Present Ian Pike
* Copyright (c) 2024-Present ccmath contributors
*
* This library is provided under the MIT License.
* See LICENSE for more information.
*/

#pragma once

#include <cstdint>
#include <limits>

#include "ccmath/internal/helpers/bits.hpp"
#include "ccmath/internal/predef/unlikely.hpp"
#include "ccmath/detail/basic/abs.hpp"
#include "ccmath/detail/basic/fmod.hpp"

namespace ccm::internal
{
namespace
{
namespace impl
{
inline constexpr double remquo_double_impl(double x, double y, int * quo) noexcept
{
std::int64_t x_i64{};
std::int64_t y_i64{};
std::uint64_t x_sign{};
std::uint64_t quotient_sign{};
int computed_quotient{};

x_i64 = ccm::helpers::double_to_int64(x);
y_i64 = ccm::helpers::double_to_int64(y);

// Determine the signs of x and the quotient.
x_sign = static_cast<std::uint64_t>(x_i64) & 0x8000000000000000ULL;
quotient_sign = x_sign ^ (static_cast<std::uint64_t>(y_i64) & 0x8000000000000000ULL);

// Clear the sign bits from the int64_t representations of x and y.
x_i64 &= 0x7fffffffffffffffULL;
y_i64 &= 0x7fffffffffffffffULL;

// If y is zero.
if (CCM_UNLIKELY(y_i64 == 0)) { return (x * y) / (x * y); }

// If x is not finite or y is NaN.
if (CCM_UNLIKELY(x_i64 >= 0x7ff0000000000000ULL || y_i64 > 0x7ff0000000000000ULL)) { return (x * y) / (x * y); }

// b (or bit 54) represents the highest bit we can compare against for both signed and unsigned integers without causing an overflow.
// Here we are checking that if y_i64 is within the range of signed 64-bit integers that can be represented without setting the MSB (i.e.,
// positive or zero).
if (y_i64 <= 0x7fbfffffffffffffULL)
{
x = ccm::fmod(x, 8 * y); // now x < (8 * y)
}

if (CCM_UNLIKELY(x_i64 == y_i64))
{
*quo = quotient_sign ? -1 : 1;
return 0.0 * x;
}

x = ccm::fabs(x);
y = ccm::helpers::int64_to_double(y_i64);
computed_quotient = 0;

if (y_i64 <= 0x7fcfffffffffffffULL && x >= 4 * y)
{
x -= 4 * y;
computed_quotient += 4;
}

if (y_i64 <= 0x7fdfffffffffffffULL && x >= 2 * y)
{
x -= 2 * y;
computed_quotient += 2;
}

if (y_i64 < 0x0020000000000000ULL)
{
if (x + x > y)
{
x -= y;
++computed_quotient;
if (x + x >= y)
{
x -= y;
++computed_quotient;
}
}
}
else
{
double y_half = 0.5 * y;
if (x > y_half)
{
x -= y;
++computed_quotient;
if (x >= y_half)
{
x -= y;
++computed_quotient;
}
}
}

*quo = quotient_sign ? -computed_quotient : computed_quotient;

// Make sure that the correct sign of zero results in round down mode.
if (x == 0.0) { x = 0.0; }
if (x_sign) { x = -x; }

return x;
}
} // namespace impl
} // namespace

template <typename T, typename = std::enable_if_t<std::is_floating_point_v<T>>>
inline constexpr T remquo_double(T x, T y, int * quo) noexcept
{
return static_cast<T>(impl::remquo_double_impl(x, y, quo));
}
} // namespace ccm::internal
118 changes: 118 additions & 0 deletions include/ccmath/detail/basic/impl/remquo_float_impl.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,118 @@
/*
* Copyright (c) 2024-Present Ian Pike
* Copyright (c) 2024-Present ccmath contributors
*
* This library is provided under the MIT License.
* See LICENSE for more information.
*/

#pragma once

#include <cstdint>
#include <limits>

#include "ccmath/internal/helpers/bits.hpp"
#include "ccmath/internal/predef/unlikely.hpp"
#include "ccmath/detail/basic/abs.hpp"
#include "ccmath/detail/basic/fmod.hpp"

namespace ccm::internal
{
namespace
{
namespace impl
{
inline constexpr float remquo_float_impl(float x, float y, int * quo) noexcept
{
std::int32_t x_i32{};
std::int32_t y_i32{};
std::uint32_t x_sign{};
int quotient_sign{};
int computed_quotient{};

x_i32 = ccm::helpers::float_to_int32(x);
y_i32 = ccm::helpers::float_to_int32(y);

// Determine the signs of x and the quotient.
x_sign = static_cast<std::uint32_t>(x_i32) & 0x80000000;
quotient_sign = x_sign ^ (static_cast<std::uint32_t>(y_i32) & 0x80000000);

// Clear the sign bits from the int32_t representations of x and y.
x_i32 &= 0x7fffffff;
y_i32 &= 0x7fffffff;

// If y is zero.
if (CCM_UNLIKELY(y_i32 == 0)) { return (x * y) / (x * y); }

// If x is not finite or y is NaN.
if (CCM_UNLIKELY(x_i32 >= 0x7f800000 || y_i32 > 0x7f800000)) { return (x * y) / (x * y); }

if (y_i32 <= 0x7dffffff)
{
x = ccm::fmod(x, 8 * y); // now x < (8 * y)
}

if ((x_i32 - y_i32) == 0)
{
*quo = quotient_sign ? -1 : 1;
return 0.0f * x;
}

x = ccm::fabsf(x);
y = ccm::fabsf(y);
computed_quotient = 0;

if (y_i32 <= 0x7e7fffff && x >= 4 * y)
{
x -= 4 * y;
computed_quotient += 4;
}
if (y_i32 <= 0x7effffff && x >= 2 * y)
{
x -= 2 * y;
computed_quotient += 2;
}
if (y_i32 < 0x01000000)
{
if (x + x > y)
{
x -= y;
++computed_quotient;
if (x + x >= y)
{
x -= y;
++computed_quotient;
}
}
}
else
{
float y_half = 0.5 * y;
if (x > y_half)
{
x -= y;
++computed_quotient;
if (x >= y_half)
{
x -= y;
++computed_quotient;
}
}
}
*quo = quotient_sign ? -computed_quotient : computed_quotient;

// Make sure that the correct sign of zero results in round down mode.
if (x == 0.0f) { x = 0.0f; }
if (x_sign) { x = -x; }

return x;
}
} // namespace impl
} // namespace

template <typename T, typename = std::enable_if_t<std::is_floating_point_v<T>>>
inline constexpr T remquo_float(T x, T y, int * quo) noexcept
{
return static_cast<T>(impl::remquo_float_impl(x, y, quo));
}
} // namespace ccm::internal
55 changes: 55 additions & 0 deletions include/ccmath/detail/basic/impl/remquo_ldouble_impl.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
/*
* Copyright (c) 2024-Present Ian Pike
* Copyright (c) 2024-Present ccmath contributors
*
* This library is provided under the MIT License.
* See LICENSE for more information.
*/

#pragma once

#include <cstdint>
#include <limits>

#include "ccmath/internal/helpers/bits.hpp"
#include "ccmath/internal/predef/unlikely.hpp"
#include "ccmath/detail/basic/abs.hpp"
#include "ccmath/detail/basic/fmod.hpp"

namespace ccm::internal
{
namespace
{
namespace impl
{
inline constexpr long double remquo_ldouble_impl(long double x, long double y, int * quo) noexcept
{
std::int32_t x_exponent{};
std::int32_t y_exponent{};
std::int32_t x_int32{};
std::int32_t y_int32{};

std::uint32_t x_sign{};
std::uint32_t y_msb{};
std::uint32_t x_lsb{};

int quotient_sign{};
int computed_quotient{};


// TODO: For the time being this is mega on hold until we have access to
// necessary implementations of int128_t. Without them extracting the
// needed bits to perform the necessary operations would be extremely challenging
// due to the fact that long double is extremely platform dependent.

return 0;
}
} // namespace impl
} // namespace

template <typename T, typename = std::enable_if_t<std::is_floating_point_v<T>>>
inline constexpr T remquo_ldouble(T x, T y, int * quo) noexcept
{
return static_cast<T>(impl::remquo_ldouble_impl(x, y, quo));
}
} // namespace ccm::internal
19 changes: 19 additions & 0 deletions include/ccmath/internal/type/int128.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
/*
* Copyright (c) 2024-Present Ian Pike
* Copyright (c) 2024-Present ccmath contributors
*
* This library is provided under the MIT License.
* See LICENSE for more information.
*/

#pragma once

namespace ccm
{
// TODO: Implement int128_t at some point.
// REASONING: Having access to int128_t would be extremely helpful in many cases when implementing specific
// algorithms such as remquo and many trigonometric functions.
//
// For the time being this is mega on hold, but I may return to this at some point.
// This file is really just here to keep a reminder of the idea.
}
2 changes: 0 additions & 2 deletions include/ccmath/math/power/sqrt.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,6 @@

#pragma once

#include <cfenv> // To get rounding mode

namespace ccm
{

Expand Down

0 comments on commit 20409c4

Please sign in to comment.