FXP: add fxp_mul_frac.

This can multiply size_ts by a fraction without the risk of overflow.
This commit is contained in:
David Goldblatt 2020-12-08 13:22:59 -08:00 committed by David Goldblatt
parent 56e85c0e47
commit caef4c2868
2 changed files with 50 additions and 0 deletions

View File

@ -90,6 +90,31 @@ fxp_round_nearest(fxp_t a) {
return (a >> 16) + increment;
}
/*
* Approximately computes x * frac, without the size limitations that would be
* imposed by converting u to an fxp_t.
*/
static inline size_t
fxp_mul_frac(size_t x_orig, fxp_t frac) {
assert(frac <= (1U << 16));
/*
* Work around an over-enthusiastic warning about type limits below (on
* 32-bit platforms, a size_t is always less than 1ULL << 48).
*/
uint64_t x = (uint64_t)x_orig;
/*
* If we can guarantee no overflow, multiply first before shifting, to
* preserve some precision. Otherwise, shift first and then multiply.
* In the latter case, we only lose the low 16 bits of a 48-bit number,
* so we're still accurate to within 1/2**32.
*/
if (x < (1ULL << 48)) {
return (size_t)((x * frac) >> 16);
} else {
return (size_t)((x >> 16) * (uint64_t)frac);
}
}
/*
* Returns true on error. Otherwise, returns false and updates *ptr to point to
* the first character not parsed (because it wasn't a digit).

View File

@ -222,6 +222,30 @@ TEST_BEGIN(test_round_simple) {
}
TEST_END
static void
expect_mul_frac(size_t a, const char *fracstr, size_t expected) {
fxp_t frac = xparse_fxp(fracstr);
size_t result = fxp_mul_frac(a, frac);
expect_true(double_close(expected, result),
"Expected %zu * %s == %zu (fracmul); got %zu", a, fracstr,
expected, result);
}
TEST_BEGIN(test_mul_frac_simple) {
expect_mul_frac(SIZE_MAX, "1.0", SIZE_MAX);
expect_mul_frac(SIZE_MAX, ".75", SIZE_MAX / 4 * 3);
expect_mul_frac(SIZE_MAX, ".5", SIZE_MAX / 2);
expect_mul_frac(SIZE_MAX, ".25", SIZE_MAX / 4);
expect_mul_frac(1U << 16, "1.0", 1U << 16);
expect_mul_frac(1U << 30, "0.5", 1U << 29);
expect_mul_frac(1U << 30, "0.25", 1U << 28);
expect_mul_frac(1U << 30, "0.125", 1U << 27);
expect_mul_frac((1U << 30) + 1, "0.125", 1U << 27);
expect_mul_frac(100, "0.25", 25);
expect_mul_frac(1000 * 1000, "0.001", 1000);
}
TEST_END
static void
expect_print(const char *str) {
fxp_t fxp = xparse_fxp(str);
@ -339,6 +363,7 @@ main(void) {
test_mul_simple,
test_div_simple,
test_round_simple,
test_mul_frac_simple,
test_print_simple,
test_stress);
}