Skip to content

Conversation

bojle
Copy link
Contributor

@bojle bojle commented Aug 22, 2025

This PR includes only one of the fxdivi functions (rdivi). It uses a polynomial function for initial approximation followed by 4 newton-raphson iterations to calculate the reciprocal and finally multiplies the numerator with it to get the result.

@lntue @PiJoules

Copy link

Thank you for submitting a Pull Request (PR) to the LLVM Project!

This PR will be automatically labeled and the relevant teams will be notified.

If you wish to, you can add reviewers by using the "Reviewers" section on this page.

If this is not working for you, it is probably because you do not have write permissions for the repository. In which case you can instead tag reviewers by name in a comment by using @ followed by their GitHub username.

If you have received no comments on your PR for a week, you can request a review by "ping"ing the PR by adding a comment “Ping”. The common courtesy "ping" rate is once a week. Please remember that you are asking for valuable time from other developers.

If you have further questions, they may be answered by the LLVM GitHub User Guide.

You can also ask questions in a comment on this PR, on the LLVM Discord or on the forums.

@llvmbot
Copy link
Member

llvmbot commented Aug 22, 2025

@llvm/pr-subscribers-libc

Author: Shreeyash Pandey (bojle)

Changes

This PR includes only one of the fxdivi functions (rdivi). It uses a polynomial function for initial approximation followed by 4 newton-raphson iterations to calculate the reciprocal and finally multiplies the numerator with it to get the result.


Full diff: https://github.com/llvm/llvm-project/pull/154914.diff

6 Files Affected:

  • (modified) libc/config/linux/riscv/entrypoints.txt (+1)
  • (modified) libc/config/linux/x86_64/entrypoints.txt (+1)
  • (modified) libc/include/stdfix.yaml (+8)
  • (modified) libc/src/__support/fixed_point/fx_bits.h (+41)
  • (modified) libc/src/stdfix/CMakeLists.txt (+14)
  • (modified) libc/test/src/stdfix/CMakeLists.txt (+16)
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 1a03683d72e61..b783dd4b04c74 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -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()
 
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5590f1a15ac57..f6beccc7229dd 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -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()
 
diff --git a/libc/include/stdfix.yaml b/libc/include/stdfix.yaml
index 5b385124eb63d..451330c3478d2 100644
--- a/libc/include/stdfix.yaml
+++ b/libc/include/stdfix.yaml
@@ -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
diff --git a/libc/src/__support/fixed_point/fx_bits.h b/libc/src/__support/fixed_point/fx_bits.h
index 00c6119b4f353..13eecba981664 100644
--- a/libc/src/__support/fixed_point/fx_bits.h
+++ b/libc/src/__support/fixed_point/fx_bits.h
@@ -224,6 +224,47 @@ idiv(T x, T y) {
   return static_cast<XType>(result);
 }
 
+inline long accum nrstep(long accum d, long accum x0) {
+  auto v = x0 * (2 - (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);
+
+  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 */
+  long accum initial_approx = a - (b * d_scaled_val);
+  long accum val = nrstep(d_scaled_val, initial_approx);
+  val = nrstep(d_scaled_val, val);
+  val = nrstep(d_scaled_val, val);
+  val = nrstep(d_scaled_val, val);
+  long accum res = n_scaled_val * val;
+  return static_cast<XType>(res);
+}
+
 } // namespace fixed_point
 } // namespace LIBC_NAMESPACE_DECL
 
diff --git a/libc/src/stdfix/CMakeLists.txt b/libc/src/stdfix/CMakeLists.txt
index 843111e3f80b1..3cbabd17ad34c 100644
--- a/libc/src/stdfix/CMakeLists.txt
+++ b/libc/src/stdfix/CMakeLists.txt
@@ -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
diff --git a/libc/test/src/stdfix/CMakeLists.txt b/libc/test/src/stdfix/CMakeLists.txt
index e2b4bc1805f7c..741522276feaa 100644
--- a/libc/test/src/stdfix/CMakeLists.txt
+++ b/libc/test/src/stdfix/CMakeLists.txt
@@ -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

@llvmbot
Copy link
Member

llvmbot commented Aug 22, 2025

@llvm/pr-subscribers-backend-risc-v

Author: Shreeyash Pandey (bojle)

Changes

This PR includes only one of the fxdivi functions (rdivi). It uses a polynomial function for initial approximation followed by 4 newton-raphson iterations to calculate the reciprocal and finally multiplies the numerator with it to get the result.


Full diff: https://github.com/llvm/llvm-project/pull/154914.diff

6 Files Affected:

  • (modified) libc/config/linux/riscv/entrypoints.txt (+1)
  • (modified) libc/config/linux/x86_64/entrypoints.txt (+1)
  • (modified) libc/include/stdfix.yaml (+8)
  • (modified) libc/src/__support/fixed_point/fx_bits.h (+41)
  • (modified) libc/src/stdfix/CMakeLists.txt (+14)
  • (modified) libc/test/src/stdfix/CMakeLists.txt (+16)
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 1a03683d72e61..b783dd4b04c74 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -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()
 
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5590f1a15ac57..f6beccc7229dd 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -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()
 
diff --git a/libc/include/stdfix.yaml b/libc/include/stdfix.yaml
index 5b385124eb63d..451330c3478d2 100644
--- a/libc/include/stdfix.yaml
+++ b/libc/include/stdfix.yaml
@@ -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
diff --git a/libc/src/__support/fixed_point/fx_bits.h b/libc/src/__support/fixed_point/fx_bits.h
index 00c6119b4f353..13eecba981664 100644
--- a/libc/src/__support/fixed_point/fx_bits.h
+++ b/libc/src/__support/fixed_point/fx_bits.h
@@ -224,6 +224,47 @@ idiv(T x, T y) {
   return static_cast<XType>(result);
 }
 
+inline long accum nrstep(long accum d, long accum x0) {
+  auto v = x0 * (2 - (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);
+
+  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 */
+  long accum initial_approx = a - (b * d_scaled_val);
+  long accum val = nrstep(d_scaled_val, initial_approx);
+  val = nrstep(d_scaled_val, val);
+  val = nrstep(d_scaled_val, val);
+  val = nrstep(d_scaled_val, val);
+  long accum res = n_scaled_val * val;
+  return static_cast<XType>(res);
+}
+
 } // namespace fixed_point
 } // namespace LIBC_NAMESPACE_DECL
 
diff --git a/libc/src/stdfix/CMakeLists.txt b/libc/src/stdfix/CMakeLists.txt
index 843111e3f80b1..3cbabd17ad34c 100644
--- a/libc/src/stdfix/CMakeLists.txt
+++ b/libc/src/stdfix/CMakeLists.txt
@@ -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
diff --git a/libc/test/src/stdfix/CMakeLists.txt b/libc/test/src/stdfix/CMakeLists.txt
index e2b4bc1805f7c..741522276feaa 100644
--- a/libc/test/src/stdfix/CMakeLists.txt
+++ b/libc/test/src/stdfix/CMakeLists.txt
@@ -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

@PiJoules
Copy link
Contributor

Thanks. Could you add some tests for this?

@bojle
Copy link
Contributor Author

bojle commented Aug 22, 2025

Oops, forgot to add them. I've pushed the tests too. Please check now. @PiJoules

Copy link

github-actions bot commented Aug 26, 2025

✅ With the latest revision this PR passed the C/C++ code formatter.

Copy link
Contributor

@PiJoules PiJoules left a comment

Choose a reason for hiding this comment

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

I think we're almost there. Just a couple more changes to test I'd like to see for more confidence.

Copy link
Contributor

@PiJoules PiJoules left a comment

Choose a reason for hiding this comment

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

LGTM. Lets wait for @lntue to take another look before submitting.

@bojle
Copy link
Contributor Author

bojle commented Sep 12, 2025

Hi @PiJoules,
If @lntue isn't available, can we have someone else from libc take a look at this?
Once this is merged, I plan to start and finish other fxdivi implementations.

Comment on lines 262 to 263
long accum a = 2.8235lk; /* 48/17 */
long accum b = 1.8823lk; /* 32/17 */
Copy link
Contributor

Choose a reason for hiding this comment

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

Use hexadecimal to be precise. 5-decimal-digit precision is not enough to specified 32-bit fractional part of long accum

@lntue
Copy link
Contributor

lntue commented Sep 12, 2025

From clang-format CI:

error: clang-format exited with code 1
diff --git a/libc/src/stdfix/rdivi.cpp b/libc/src/stdfix/rdivi.cpp
index 56de3eb3a..7021a8b1c 100644
--- a/libc/src/stdfix/rdivi.cpp
+++ b/libc/src/stdfix/rdivi.cpp
@@ -6,11 +6,11 @@
 //
 //===----------------------------------------------------------------------===//
 
-#include "include/llvm-libc-macros/stdfix-macros.h" // fract
 #include "rdivi.h"
-#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
+#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 {

Signed-off-by: Shreeyash Pandey <[email protected]>
Signed-off-by: Shreeyash Pandey <[email protected]>
Signed-off-by: Shreeyash Pandey <[email protected]>
Comment on lines 284 to 288
// E2 = 0.0000121
val = nrstep(d_scaled_val, val);
// E3 = 1.468e−10
val = nrstep(d_scaled_val, val);

Copy link
Contributor

Choose a reason for hiding this comment

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

We can skip further Newton-Raphson steps for smaller types:

if constexpr (FXRep<XType>::FRACTION_LENGTH > 8) {
  val = nrstep(...);
  if constexpr(FXRep<XType>::FRACTION_LENGTH > 16)
    val = nrstep(...);
}

Copy link
Contributor Author

@bojle bojle Sep 27, 2025

Choose a reason for hiding this comment

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

Doing this, I end up losing precision for finite result tests. For example, 128/256 results in 0.499969. With 3 iterations, the result is precisely 0.5. If the results, even for finite tests is acceptable to be +- epsilon, I can re-write the tests to do EXPECT_LT instead of EXPECT_EQ

Copy link
Contributor

Choose a reason for hiding this comment

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

Theoretically, to get correct rounding for division, we will need twice the precision. So with 32-bit inputs, we will need 64-bit precision before rounding to the targeted format. In here, we are actually operating on a bit more than 32-bit precision, so there will be results that are not correctly rounded.

Nonetheless, dividing a power of 2 not correctly might be a surprise to people, so maybe we can add a special branch for the case d is a power of 2. WDYT?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Since this is a performance tradeoff, having performance numbers would be a
good idea. The difference is 3 inline function calls back to back vs a for
loop with a bound only known at runtime. The former has a clear dataflow with
no branches (pipeline friendly?), and latter has a cmp and branch with unknown
loop count. We could measure both, then make this decision?

Copy link
Contributor Author

@bojle bojle Oct 4, 2025

Choose a reason for hiding this comment

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

I made these mods to the code:

  long accum val = nrstep(d_scaled_val, initial_approx);
  if (isPowerOfTwo(d)) {
    // E2 = 0.0000121
    val = nrstep(d_scaled_val, val);
    // E3 = 1.468e−10
    val = nrstep(d_scaled_val, val);
  }

and set up a naive time measurement script.

The results:

With contiguos 3 iterations (dividing 76/128 for 10k iterations): 650us
With if statements for po2 (76/128 for 10k iterations): 526us

So, I think I should add checks for po2.

Signed-off-by: Shreeyash Pandey <[email protected]>
// Division with a power of 2 would generally be expected to be
// exact, we handle this by specially treating po2 cases and having
// extra iterations for them.
if (FXRep<XType>::FRACTION_LEN > 8 || isPowerOfTwo(cpp::abs(d))) {
Copy link
Contributor

Choose a reason for hiding this comment

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

when the denominator is a power of two, it's better to just perform right shifts (be careful with signed-ness and overflow). Then the other conditions on the fraction lengths can be constexpr, i.e. no runtime cost.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants