From 12e5042ece6bc78f00facb39bdbe47799cecc533 Mon Sep 17 00:00:00 2001
From: Paul Davis <paul@linuxaudiosystems.com>
Date: Fri, 27 May 2022 12:45:47 -0600
Subject: [PATCH] libpbd: add muldiv() to compute v * (n/d) without overflow

---
 libs/pbd/pbd/integer_division.h | 58 +++++++++++++++++++++++++++++++++
 1 file changed, 58 insertions(+)

diff --git a/libs/pbd/pbd/integer_division.h b/libs/pbd/pbd/integer_division.h
index 323fbf03c8..e2fe056397 100644
--- a/libs/pbd/pbd/integer_division.h
+++ b/libs/pbd/pbd/integer_division.h
@@ -19,6 +19,13 @@
 #ifndef __libpbd_integer_division_h__
 #define __libpbd_integer_division_h__
 
+#include <cstdint>
+
+#ifndef COMPILER_INT128_SUPPORT
+#include <boost/multiprecision/cpp_int.hpp>
+#include "pbd/error.h"
+#endif
+
 #define PBD_IDIV_ASR(x) ((x) < 0 ? -1 : 0)  // Compiles into a (N-1)-bit arithmetic shift right
 
 /* The value of PBD_IDIV_ROUNDING will have the same sign as the dividend (x) and half
@@ -36,4 +43,55 @@ T int_div_round (T x, T y)
 	return (x + PBD_IDIV_ROUNDING(x,y)) / y ;
 }
 
+namespace PBD {
+
+/* this computes v * (n/d) where v, n and d are all 64 bit integers, without
+ * overflow, and with appropriate rounding given that this is integer division.
+ */
+
+inline
+int64_t muldiv (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 /= d;
+
+	try {
+
+		return bignum.convert_to<int64_t> ();
+
+	} catch (...) {
+		fatal << "arithmetic overflow in timeline math\n" << endmsg;
+		/* NOTREACHED */
+		return 0;
+	}
+
+#else
+	__int128 _n (n);
+	__int128 _d (d);
+	__int128 _v (v);
+
+	/* this could overflow, but will not do so merely because we are
+	 * multiplying two int64_t together and storing the result in an
+	 * int64_t. Overflow will occur where (v*n)+hd > INT128_MAX (hard
+	 * limit) or where v * n / d > INT64_T (i.e. n > d)
+	 */
+
+	return(int64_t) (((_v * _n) + hd) / _d);
+#endif
+}
+} /* namespace */
+
 #endif /* __libpbd_integer_division_h___ */