From b29a62d87cc0af3e9d134e9e0863b2cb053070b8 Mon Sep 17 00:00:00 2001 From: Nicolas Pitre Date: Sun, 7 Jul 2024 15:05:19 -0400 Subject: mul_u64_u64_div_u64: make it precise always MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Patch series "mul_u64_u64_div_u64: new implementation", v3. This provides an implementation for mul_u64_u64_div_u64() that always produces exact results. This patch (of 2): Library facilities must always return exact results. If the caller may be contented with approximations then it should do the approximation on its own. In this particular case the comment in the code says "the algorithm ... below might lose some precision". Well, if you try it with e.g.: a = 18446462598732840960 b = 18446462598732840960 c = 18446462598732840961 then the produced answer is 0 whereas the exact answer should be 18446462598732840959. This is _some_ precision lost indeed! Let's reimplement this function so it always produces the exact result regardless of its inputs while preserving existing fast paths when possible. Uwe said: : My personal interest is to get the calculations in pwm drivers right. : This function is used in several drivers below drivers/pwm/ . With the : errors in mul_u64_u64_div_u64(), pwm consumers might not get the : settings they request. Although I have to admit that I'm not aware it : breaks real use cases (because typically the periods used are too short : to make the involved multiplications overflow), but I pretty sure am : not aware of all usages and it breaks testing. : : Another justification is commits like : https://git.kernel.org/tip/77baa5bafcbe1b2a15ef9c37232c21279c95481c, : where people start to work around the precision shortcomings of : mul_u64_u64_div_u64(). Link: https://lkml.kernel.org/r/20240707190648.1982714-1-nico@fluxnic.net Link: https://lkml.kernel.org/r/20240707190648.1982714-2-nico@fluxnic.net Signed-off-by: Nicolas Pitre Tested-by: Uwe Kleine-König Reviewed-by: Uwe Kleine-König Tested-by: Biju Das Signed-off-by: Andrew Morton --- lib/math/div64.c | 108 +++++++++++++++++++++++++++++++++---------------------- 1 file changed, 65 insertions(+), 43 deletions(-) (limited to 'lib/math/div64.c') diff --git a/lib/math/div64.c b/lib/math/div64.c index 191761b1b623..b7fc75246399 100644 --- a/lib/math/div64.c +++ b/lib/math/div64.c @@ -186,55 +186,77 @@ EXPORT_SYMBOL(iter_div_u64_rem); #ifndef mul_u64_u64_div_u64 u64 mul_u64_u64_div_u64(u64 a, u64 b, u64 c) { - u64 res = 0, div, rem; - int shift; + if (ilog2(a) + ilog2(b) <= 62) + return div64_u64(a * b, c); - /* can a * b overflow ? */ - if (ilog2(a) + ilog2(b) > 62) { - /* - * Note that the algorithm after the if block below might lose - * some precision and the result is more exact for b > a. So - * exchange a and b if a is bigger than b. - * - * For example with a = 43980465100800, b = 100000000, c = 1000000000 - * the below calculation doesn't modify b at all because div == 0 - * and then shift becomes 45 + 26 - 62 = 9 and so the result - * becomes 4398035251080. However with a and b swapped the exact - * result is calculated (i.e. 4398046510080). - */ - if (a > b) - swap(a, b); +#if defined(__SIZEOF_INT128__) + + /* native 64x64=128 bits multiplication */ + u128 prod = (u128)a * b; + u64 n_lo = prod, n_hi = prod >> 64; + +#else + + /* perform a 64x64=128 bits multiplication manually */ + u32 a_lo = a, a_hi = a >> 32, b_lo = b, b_hi = b >> 32; + u64 x, y, z; + + x = (u64)a_lo * b_lo; + y = (u64)a_lo * b_hi + (u32)(x >> 32); + z = (u64)a_hi * b_hi + (u32)(y >> 32); + y = (u64)a_hi * b_lo + (u32)y; + z += (u32)(y >> 32); + x = (y << 32) + (u32)x; + + u64 n_lo = x, n_hi = z; + +#endif + + int shift = __builtin_ctzll(c); + /* try reducing the fraction in case the dividend becomes <= 64 bits */ + if ((n_hi >> shift) == 0) { + u64 n = (n_lo >> shift) | (n_hi << (64 - shift)); + + return div64_u64(n, c >> shift); /* - * (b * a) / c is equal to - * - * (b / c) * a + - * (b % c) * a / c - * - * if nothing overflows. Can the 1st multiplication - * overflow? Yes, but we do not care: this can only - * happen if the end result can't fit in u64 anyway. - * - * So the code below does - * - * res = (b / c) * a; - * b = b % c; + * The remainder value if needed would be: + * res = div64_u64_rem(n, c >> shift, &rem); + * rem = (rem << shift) + (n_lo - (n << shift)); */ - div = div64_u64_rem(b, c, &rem); - res = div * a; - b = rem; - - shift = ilog2(a) + ilog2(b) - 62; - if (shift > 0) { - /* drop precision */ - b >>= shift; - c >>= shift; - if (!c) - return res; - } } - return res + div64_u64(a * b, c); + if (n_hi >= c) { + /* overflow: result is unrepresentable in a u64 */ + return -1; + } + + /* Do the full 128 by 64 bits division */ + + shift = __builtin_clzll(c); + c <<= shift; + + int p = 64 + shift; + u64 res = 0; + bool carry; + + do { + carry = n_hi >> 63; + shift = carry ? 1 : __builtin_clzll(n_hi); + if (p < shift) + break; + p -= shift; + n_hi <<= shift; + n_hi |= n_lo >> (64 - shift); + n_lo <<= shift; + if (carry || (n_hi >= c)) { + n_hi -= c; + res |= 1ULL << p; + } + } while (n_hi); + /* The remainder value if needed would be n_hi << p */ + + return res; } EXPORT_SYMBOL(mul_u64_u64_div_u64); #endif -- cgit v1.2.3-73-gaa49b From 1635e62e75a7bbb1c6274f6b43911cedfe0da60a Mon Sep 17 00:00:00 2001 From: Nicolas Pitre Date: Sun, 7 Jul 2024 15:05:20 -0400 Subject: mul_u64_u64_div_u64: basic sanity test MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Verify that edge cases produce proper results, and some more. [npitre@baylibre.com: avoid undefined shift value] Link: https://lkml.kernel.org/r/7rrs9pn1-n266-3013-9q6n-1osp8r8s0rrn@syhkavp.arg Link: https://lkml.kernel.org/r/20240707190648.1982714-3-nico@fluxnic.net Signed-off-by: Nicolas Pitre Reviewed-by: Uwe Kleine-König Cc: Biju Das Signed-off-by: Andrew Morton --- lib/Kconfig.debug | 10 ++++ lib/math/Makefile | 1 + lib/math/div64.c | 9 +++- lib/math/test_mul_u64_u64_div_u64.c | 99 +++++++++++++++++++++++++++++++++++++ 4 files changed, 118 insertions(+), 1 deletion(-) create mode 100644 lib/math/test_mul_u64_u64_div_u64.c (limited to 'lib/math/div64.c') diff --git a/lib/Kconfig.debug b/lib/Kconfig.debug index a30c03a66172..bf0995d328b3 100644 --- a/lib/Kconfig.debug +++ b/lib/Kconfig.debug @@ -2280,6 +2280,16 @@ config TEST_DIV64 If unsure, say N. +config TEST_MULDIV64 + tristate "mul_u64_u64_div_u64() test" + depends on DEBUG_KERNEL || m + help + Enable this to turn on 'mul_u64_u64_div_u64()' function test. + This test is executed only once during system boot (so affects + only boot time), or at module load time. + + If unsure, say N. + config TEST_IOV_ITER tristate "Test iov_iter operation" if !KUNIT_ALL_TESTS depends on KUNIT diff --git a/lib/math/Makefile b/lib/math/Makefile index 91fcdb0c9efe..981a26127e08 100644 --- a/lib/math/Makefile +++ b/lib/math/Makefile @@ -6,4 +6,5 @@ obj-$(CONFIG_PRIME_NUMBERS) += prime_numbers.o obj-$(CONFIG_RATIONAL) += rational.o obj-$(CONFIG_TEST_DIV64) += test_div64.o +obj-$(CONFIG_TEST_MULDIV64) += test_mul_u64_u64_div_u64.o obj-$(CONFIG_RATIONAL_KUNIT_TEST) += rational-test.o diff --git a/lib/math/div64.c b/lib/math/div64.c index b7fc75246399..5faa29208bdb 100644 --- a/lib/math/div64.c +++ b/lib/math/div64.c @@ -212,11 +212,18 @@ u64 mul_u64_u64_div_u64(u64 a, u64 b, u64 c) #endif + /* make sure c is not zero, trigger exception otherwise */ +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wdiv-by-zero" + if (unlikely(c == 0)) + return 1/0; +#pragma GCC diagnostic pop + int shift = __builtin_ctzll(c); /* try reducing the fraction in case the dividend becomes <= 64 bits */ if ((n_hi >> shift) == 0) { - u64 n = (n_lo >> shift) | (n_hi << (64 - shift)); + u64 n = shift ? (n_lo >> shift) | (n_hi << (64 - shift)) : n_lo; return div64_u64(n, c >> shift); /* diff --git a/lib/math/test_mul_u64_u64_div_u64.c b/lib/math/test_mul_u64_u64_div_u64.c new file mode 100644 index 000000000000..58d058de4e73 --- /dev/null +++ b/lib/math/test_mul_u64_u64_div_u64.c @@ -0,0 +1,99 @@ +// SPDX-License-Identifier: GPL-2.0 +/* + * Copyright (C) 2024 BayLibre SAS + */ + +#define pr_fmt(fmt) KBUILD_MODNAME ": " fmt + +#include +#include +#include +#include + +typedef struct { u64 a; u64 b; u64 c; u64 result; } test_params; + +static test_params test_values[] = { +/* this contains many edge values followed by a couple random values */ +{ 0xb, 0x7, 0x3, 0x19 }, +{ 0xffff0000, 0xffff0000, 0xf, 0x1110eeef00000000 }, +{ 0xffffffff, 0xffffffff, 0x1, 0xfffffffe00000001 }, +{ 0xffffffff, 0xffffffff, 0x2, 0x7fffffff00000000 }, +{ 0x1ffffffff, 0xffffffff, 0x2, 0xfffffffe80000000 }, +{ 0x1ffffffff, 0xffffffff, 0x3, 0xaaaaaaa9aaaaaaab }, +{ 0x1ffffffff, 0x1ffffffff, 0x4, 0xffffffff00000000 }, +{ 0xffff000000000000, 0xffff000000000000, 0xffff000000000001, 0xfffeffffffffffff }, +{ 0x3333333333333333, 0x3333333333333333, 0x5555555555555555, 0x1eb851eb851eb851 }, +{ 0x7fffffffffffffff, 0x2, 0x3, 0x5555555555555554 }, +{ 0xffffffffffffffff, 0x2, 0x8000000000000000, 0x3 }, +{ 0xffffffffffffffff, 0x2, 0xc000000000000000, 0x2 }, +{ 0xffffffffffffffff, 0x4000000000000004, 0x8000000000000000, 0x8000000000000007 }, +{ 0xffffffffffffffff, 0x4000000000000001, 0x8000000000000000, 0x8000000000000001 }, +{ 0xffffffffffffffff, 0x8000000000000001, 0xffffffffffffffff, 0x8000000000000001 }, +{ 0xfffffffffffffffe, 0x8000000000000001, 0xffffffffffffffff, 0x8000000000000000 }, +{ 0xffffffffffffffff, 0x8000000000000001, 0xfffffffffffffffe, 0x8000000000000001 }, +{ 0xffffffffffffffff, 0x8000000000000001, 0xfffffffffffffffd, 0x8000000000000002 }, +{ 0x7fffffffffffffff, 0xffffffffffffffff, 0xc000000000000000, 0xaaaaaaaaaaaaaaa8 }, +{ 0xffffffffffffffff, 0x7fffffffffffffff, 0xa000000000000000, 0xccccccccccccccca }, +{ 0xffffffffffffffff, 0x7fffffffffffffff, 0x9000000000000000, 0xe38e38e38e38e38b }, +{ 0x7fffffffffffffff, 0x7fffffffffffffff, 0x5000000000000000, 0xccccccccccccccc9 }, +{ 0xffffffffffffffff, 0xfffffffffffffffe, 0xffffffffffffffff, 0xfffffffffffffffe }, +{ 0xe6102d256d7ea3ae, 0x70a77d0be4c31201, 0xd63ec35ab3220357, 0x78f8bf8cc86c6e18 }, +{ 0xf53bae05cb86c6e1, 0x3847b32d2f8d32e0, 0xcfd4f55a647f403c, 0x42687f79d8998d35 }, +{ 0x9951c5498f941092, 0x1f8c8bfdf287a251, 0xa3c8dc5f81ea3fe2, 0x1d887cb25900091f }, +{ 0x374fee9daa1bb2bb, 0x0d0bfbff7b8ae3ef, 0xc169337bd42d5179, 0x03bb2dbaffcbb961 }, +{ 0xeac0d03ac10eeaf0, 0x89be05dfa162ed9b, 0x92bb1679a41f0e4b, 0xdc5f5cc9e270d216 }, +}; + +/* + * The above table can be verified with the following shell script: + * + * #!/bin/sh + * sed -ne 's/^{ \+\(.*\), \+\(.*\), \+\(.*\), \+\(.*\) },$/\1 \2 \3 \4/p' \ + * lib/math/test_mul_u64_u64_div_u64.c | + * while read a b c r; do + * expected=$( printf "obase=16; ibase=16; %X * %X / %X\n" $a $b $c | bc ) + * given=$( printf "%X\n" $r ) + * if [ "$expected" = "$given" ]; then + * echo "$a * $b / $c = $r OK" + * else + * echo "$a * $b / $c = $r is wrong" >&2 + * echo "should be equivalent to 0x$expected" >&2 + * exit 1 + * fi + * done + */ + +static int __init test_init(void) +{ + int i; + + pr_info("Starting mul_u64_u64_div_u64() test\n"); + + for (i = 0; i < ARRAY_SIZE(test_values); i++) { + u64 a = test_values[i].a; + u64 b = test_values[i].b; + u64 c = test_values[i].c; + u64 expected_result = test_values[i].result; + u64 result = mul_u64_u64_div_u64(a, b, c); + + if (result != expected_result) { + pr_err("ERROR: 0x%016llx * 0x%016llx / 0x%016llx\n", a, b, c); + pr_err("ERROR: expected result: %016llx\n", expected_result); + pr_err("ERROR: obtained result: %016llx\n", result); + } + } + + pr_info("Completed mul_u64_u64_div_u64() test\n"); + return 0; +} + +static void __exit test_exit(void) +{ +} + +module_init(test_init); +module_exit(test_exit); + +MODULE_AUTHOR("Nicolas Pitre"); +MODULE_LICENSE("GPL"); +MODULE_DESCRIPTION("mul_u64_u64_div_u64() test module"); -- cgit v1.2.3-73-gaa49b