Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions libc/config/linux/riscv/entrypoints.txt
Original file line number Diff line number Diff line change
Expand Up @@ -986,6 +986,7 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
libc.src.stdfix.idivulr
libc.src.stdfix.idivuk
libc.src.stdfix.idivulk
libc.src.stdfix.rdivi
)
endif()

Expand Down
1 change: 1 addition & 0 deletions libc/config/linux/x86_64/entrypoints.txt
Original file line number Diff line number Diff line change
Expand Up @@ -1019,6 +1019,7 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
libc.src.stdfix.idivulr
libc.src.stdfix.idivuk
libc.src.stdfix.idivulk
libc.src.stdfix.rdivi
)
endif()

Expand Down
8 changes: 8 additions & 0 deletions libc/include/stdfix.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -544,3 +544,11 @@ functions:
arguments:
- type: unsigned long accum
guard: LIBC_COMPILER_HAS_FIXED_POINT
- name: rdivi
standards:
- stdc_ext
return_type: fract
arguments:
- type: int
- type: int
guard: LIBC_COMPILER_HAS_FIXED_POINT
64 changes: 64 additions & 0 deletions libc/src/__support/fixed_point/fx_bits.h
Original file line number Diff line number Diff line change
Expand Up @@ -224,6 +224,70 @@ idiv(T x, T y) {
return static_cast<XType>(result);
}

LIBC_INLINE long accum nrstep(long accum d, long accum x0) {
auto v = x0 * (2.lk - (d * x0));
return v;
}

/* Divide the two integers and return a fixed_point value
*
* For reference, see:
* https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division
* https://stackoverflow.com/a/9231996
*/
template <typename XType> LIBC_INLINE constexpr XType divi(int n, int d) {
// If the value of the second operand of the / operator is zero, the
// behavior is undefined. Ref: ISO/IEC TR 18037:2008(E) p.g. 16
LIBC_CRASH_ON_VALUE(d, 0);

if (LIBC_UNLIKELY(n == 0)) {
return FXRep<XType>::ZERO();
}
bool d_is_signed = false;
if (d < 0) {
d = (d * -1);
d_is_signed = true;
}

unsigned int nv = static_cast<unsigned int>(n);
unsigned int dv = static_cast<unsigned int>(d);
unsigned int clz = cpp::countl_zero<unsigned int>(dv) - 1;
unsigned long int scaled_val = dv << clz;
/* Scale denominator to be in the range of [0.5,1] */
FXBits<long accum> d_scaled{scaled_val};
unsigned long int scaled_val_n = nv << clz;
/* Scale the numerator as much as the denominator to maintain correctness of
* the original equation
*/
FXBits<long accum> n_scaled{scaled_val_n};
long accum n_scaled_val = n_scaled.get_val();
long accum d_scaled_val = d_scaled.get_val();
/* x0 = (48/17) - (32/17) * d_n */
long accum a = 2.8235lk; /* 48/17 */
long accum b = 1.8823lk; /* 32/17 */
/* Error of the initial approximation, as derived
* from the wikipedia article is
* E0 = 1/17 = 0.059 (5.9%)
*/
long accum initial_approx = a - (b * d_scaled_val);
/* Each newton-raphson iteration will square the error, due
* to quadratic convergence. So,
* E1 = (0.059)^2 = 0.0034
*/
long accum val = nrstep(d_scaled_val, initial_approx);
/* E2 = 0.0000121 */
val = nrstep(d_scaled_val, val);
/* E3 = 1.468e−10 */
val = nrstep(d_scaled_val, val);
/* E4 = 2.155e−20 */
val = nrstep(d_scaled_val, val);
long accum res = n_scaled_val * val;
if (d_is_signed) {
res *= static_cast<XType>(-1);
}
return static_cast<XType>(res);
}

} // namespace fixed_point
} // namespace LIBC_NAMESPACE_DECL

Expand Down
14 changes: 14 additions & 0 deletions libc/src/stdfix/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,20 @@ foreach(suffix IN ITEMS r lr k lk ur ulr uk ulk)
)
endforeach()

foreach(suffix IN ITEMS r)
add_entrypoint_object(
${suffix}divi
HDRS
${suffix}divi.h
SRCS
${suffix}divi.cpp
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.fixed_point.fx_bits
)
endforeach()

add_entrypoint_object(
uhksqrtus
HDRS
Expand Down
21 changes: 21 additions & 0 deletions libc/src/stdfix/rdivi.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
//===-- Implementation of rdivi function ---------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//

#include "rdivi.h"
#include "include/llvm-libc-macros/stdfix-macros.h" // fract
#include "src/__support/common.h" // LLVM_LIBC_FUNCTION
#include "src/__support/fixed_point/fx_bits.h" // fixed_point
#include "src/__support/macros/config.h" // LIBC_NAMESPACE_DECL

namespace LIBC_NAMESPACE_DECL {

LLVM_LIBC_FUNCTION(fract, rdivi, (int a, int b)) {
return fixed_point::divi<fract>(a,b);
}

} // namespace LIBC_NAMESPACE_DECL
21 changes: 21 additions & 0 deletions libc/src/stdfix/rdivi.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
//===-- Implementation header for rdivi ------------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC_STDFIX_RDIVI_H
#define LLVM_LIBC_SRC_STDFIX_RDIVI_H

#include "include/llvm-libc-macros/stdfix-macros.h"
#include "src/__support/macros/config.h"

namespace LIBC_NAMESPACE_DECL {

fract rdivi(int a, int b);

} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC_STDFIX_RDIVI_H
16 changes: 16 additions & 0 deletions libc/test/src/stdfix/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,22 @@ foreach(suffix IN ITEMS r lr k lk ur ulr uk ulk)
)
endforeach()

foreach(suffix IN ITEMS r)
add_libc_test(
${suffix}divi_test
SUITE
libc-stdfix-tests
HDRS
DivITest.h
SRCS
${suffix}divi_test.cpp
DEPENDS
libc.src.stdfix.${suffix}divi
libc.src.__support.fixed_point.fx_bits
libc.hdr.signal_macros
)
endforeach()

add_libc_test(
uhksqrtus_test
SUITE
Expand Down
65 changes: 65 additions & 0 deletions libc/test/src/stdfix/DivITest.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
//===-- Utility class to test fxdivi functions ------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//

#include "src/__support/CPP/type_traits.h"
#include "src/__support/fixed_point/fx_bits.h"
#include "src/__support/fixed_point/fx_rep.h"
#include "test/UnitTest/Test.h"

template <typename XType> XType get_epsilon() = delete;
template <> fract get_epsilon() { return FRACT_EPSILON; }
template <> unsigned fract get_epsilon() { return UFRACT_EPSILON; }
template <> long fract get_epsilon() { return LFRACT_EPSILON; }

template <typename XType>
class DivITest : public LIBC_NAMESPACE::testing::Test {
using FXRep = LIBC_NAMESPACE::fixed_point::FXRep<XType>;
using FXBits = LIBC_NAMESPACE::fixed_point::FXBits<XType>;

public:
typedef XType (*DivIFunc)(int, int);

void testBasic(DivIFunc func) {
XType epsilon = get_epsilon<XType>();
EXPECT_LT((func(2, 3) - 0.666656494140625r), epsilon);
EXPECT_LT((func(3, 4) - 0.75r), epsilon);
EXPECT_LT((func(1043, 2764) - 0.3773516643r), epsilon);
EXPECT_LT((func(60000, 720293) - 0.08329943509r), epsilon);

EXPECT_EQ(func(128, 256), 0.5r);
EXPECT_EQ(func(1, 2), 0.5r);
EXPECT_EQ(func(1, 4), 0.25r);
EXPECT_EQ(func(1, 8), 0.125r);
EXPECT_EQ(func(1, 16), 0.0625r);

EXPECT_EQ(func(-1, 2), -0.5r);
EXPECT_EQ(func(1, -4), -0.25r);
EXPECT_EQ(func(-1, 8), -0.125r);
EXPECT_EQ(func(1, -16), -0.0625r);
}

void testSpecial(DivIFunc func) {
EXPECT_EQ(func(0,10), 0.r);
EXPECT_EQ(func(0,-10), 0.r);
EXPECT_EQ(func(-32768,32768), FRACT_MIN);
EXPECT_EQ(func(32767,32768), FRACT_MAX);
Comment on lines +49 to +50
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this will only be true if the number of fractional bits for a fract type is 16, but it can vary depending on the platform, so let's maybe instead derive these number of fractional bits. Something like

EXPECT_EQ(func(-1 << FRACT_FBIT, 1 << FRACT_FBIT), FRACT_MIN);
EXPECT_EQ(func((1 << FRACT_FBIT) - 1, 1 << FRACT_FBIT), FRACT_MAX);

EXPECT_EQ(func(INT_MAX,INT_MAX), 1.0r);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is correct but it had me second guessing a little because of some nuances with the spec I had to look up.

Per 7.18a.6.1

For functions returning a fixed-point value, the return value is saturated on overflow. 

so rdivi(INT_MAX,INT_MAX) should definitely return FRACT_MAX, and per Clause 6.4.4:

The value of a constant shall be in the range of representable values for its type, with exception for constants of a fract type with a value of exactly 1; such a constant shall denote the maximal value for the type.

so 1.0r should also resolve to FRACT_MAX hence the comparison is correct. Just to make things clearer, could you add these as comments above this check, or alternatively make the RHS FRACT_MAX while keeping the rdivi saturation comment?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tangentially, it also might be good to add some tests which would trigger the saturated overflow case. So maybe some integral values which would resolve to either less than zero or greater than or equal to one which would saturate to zero and FRACT_MAX.

EXPECT_EQ(func(INT_MAX-1,INT_MAX), 0.99999999r);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since 0.99999999r also can't fit precisely, let's use the EXPECT_LT epsilon approach used above.

EXPECT_EQ(func(INT_MIN,INT_MAX), FRACT_MIN);
/* Expecting 0 here as fract is not precise enough to
* handle 1/INT_MAX
*/
EXPECT_EQ(func(1, INT_MAX), 0.r);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similar to before, rounding can be either up or down, so I think we also want to use the epsilon approach here.

}
};

#define LIST_DIVI_TESTS(Name, XType, func) \
using LlvmLibc##Name##diviTest = DivITest<XType>; \
TEST_F(LlvmLibc##Name##diviTest, Basic) { testBasic(&func); } \
TEST_F(LlvmLibc##Name##diviTest, Special) { testSpecial(&func); } \
static_assert(true, "Require semicolon.")
14 changes: 14 additions & 0 deletions libc/test/src/stdfix/rdivi_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
//===-- Unittests for rdivi -----------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//

#include "DivITest.h"

#include "llvm-libc-macros/stdfix-macros.h" // fract
#include "src/stdfix/rdivi.h"

LIST_DIVI_TESTS(r, fract, LIBC_NAMESPACE::rdivi);
Loading