-
Notifications
You must be signed in to change notification settings - Fork 15.2k
[libc][stdfix] Implement fxdivi functions (rdivi) #154914
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Conversation
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 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. |
@llvm/pr-subscribers-libc Author: Shreeyash Pandey (bojle) ChangesThis 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:
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
|
@llvm/pr-subscribers-backend-risc-v Author: Shreeyash Pandey (bojle) ChangesThis 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:
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
|
Thanks. Could you add some tests for this? |
Oops, forgot to add them. I've pushed the tests too. Please check now. @PiJoules |
✅ With the latest revision this PR passed the C/C++ code formatter. |
Signed-off-by: Shreeyash Pandey <[email protected]>
Signed-off-by: Shreeyash Pandey <[email protected]>
There was a problem hiding this 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.
Signed-off-by: Shreeyash Pandey <[email protected]>
Signed-off-by: Shreeyash Pandey <[email protected]>
There was a problem hiding this 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.
long accum a = 2.8235lk; /* 48/17 */ | ||
long accum b = 1.8823lk; /* 32/17 */ |
There was a problem hiding this comment.
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
From clang-format CI:
|
Signed-off-by: Shreeyash Pandey <[email protected]>
Signed-off-by: Shreeyash Pandey <[email protected]>
Signed-off-by: Shreeyash Pandey <[email protected]>
Signed-off-by: Shreeyash Pandey <[email protected]>
Signed-off-by: Shreeyash Pandey <[email protected]>
// E2 = 0.0000121 | ||
val = nrstep(d_scaled_val, val); | ||
// E3 = 1.468e−10 | ||
val = nrstep(d_scaled_val, val); | ||
|
There was a problem hiding this comment.
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(...);
}
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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))) { |
There was a problem hiding this comment.
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.
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