Fix 128 bit integer math rounding (#9057)

The previous code only used the 1st multiplicand was use to
determine the direction of rounding, breaking commutative property

`muldiv_round (1, 3, 4) != muldiv_round (3, 1, 4)`
This commit is contained in:
Robin Gareus 2023-01-12 16:08:53 +01:00
parent 1ff5592731
commit 117cfc844b
Signed by: rgareus
GPG Key ID: A090BCE02CF57F04
1 changed files with 5 additions and 11 deletions

View File

@ -52,20 +52,11 @@ namespace PBD {
inline
int64_t muldiv_round (int64_t v, int64_t n, int64_t d)
{
/* either n or d or both could be negative but for now we assume that
only d could be (that is, n and d represent negative rational numbers of the
form 1/-2 rather than -1/2). This follows the behavior of the
ratio_t type in the temporal library.
Consequently, we only use d in the rounding-signdness expression.
*/
const int64_t hd = PBD_IDIV_ROUNDING (v, d);
#ifndef COMPILER_INT128_SUPPORT
boost::multiprecision::int512_t bignum = v;
bignum *= n;
bignum += hd;
bignum += PBD_IDIV_ROUNDING (bignum, d);
bignum /= d;
try {
@ -82,6 +73,9 @@ int64_t muldiv_round (int64_t v, int64_t n, int64_t d)
__int128 _n (n);
__int128 _d (d);
__int128 _v (v);
__int128 vn (_v * _n);
const int64_t hd = PBD_IDIV_ROUNDING (vn, d);
/* this could overflow, but will not do so merely because we are
* multiplying two int64_t together and storing the result in an
@ -89,7 +83,7 @@ int64_t muldiv_round (int64_t v, int64_t n, int64_t d)
* limit) or where v * n / d > INT64_T (i.e. n > d)
*/
return(int64_t) (((_v * _n) + hd) / _d);
return(int64_t) ((vn + hd) / _d);
#endif
}