/*
    Copyright (C) 2020 Paul Davis

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

#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
 * the magnitude of the divisor (y). Adding ROUNDING to the dividend thus
 * increases its magnitude before the integer division truncates the resulting
 * quotient.
 */

#define PBD_IDIV_ROUNDING(x,y) ( (y)/2 - (PBD_IDIV_ASR((x)^(y)) & (y)))

template<typename T>
T int_div_round (T x, T y)
{
	/* essentially ((x + (y/2)) / y) but handles signed/negative values correcvtly */
	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_round (int64_t v, int64_t n, int64_t d)
{
#ifndef COMPILER_INT128_SUPPORT
	boost::multiprecision::int512_t bignum = v;

	bignum *= n;
	bignum += PBD_IDIV_ROUNDING (bignum, d);
	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);
	__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
	 * 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) ((vn + hd) / _d);
#endif
}

inline
int64_t muldiv_floor (int64_t v, int64_t n, int64_t d)
{
#ifndef COMPILER_INT128_SUPPORT
	boost::multiprecision::int512_t bignum = v;

	bignum *= n;
	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) / _d);
#endif
}
} /* namespace */

#endif /* __libpbd_integer_division_h___ */