#include "int64x64-128.h"
#include "ns3/abort.h"
#include "ns3/assert.h"

namespace ns3 {

#define OUTPUT_SIGN(sa,sb,ua,ub)					\
  ({bool negA, negB;							\
    negA = sa < 0;							\
    negB = sb < 0;							\
    ua = negA?-sa:sa;							\
    ub = negB?-sb:sb;							\
    (negA && !negB) || (!negA && negB);})


#define MASK_LO ((((int128_t)1)<<64)-1)
#define MASK_HI (~MASK_LO)

void
int64x64_t::Mul (int64x64_t const &o)
{
  bool negResult;
  uint128_t a, b;
  negResult = OUTPUT_SIGN (_v, o._v, a, b);
  int128_t result = Umul (a, b);
  // add the sign to the result
  result = negResult ? -result : result;
  _v = result;
}

uint128_t
int64x64_t::Umul (uint128_t a, uint128_t b)
{
  uint128_t aL = a & MASK_LO;
  uint128_t bL = b & MASK_LO;
  uint128_t aH = (a >> 64) & MASK_LO;
  uint128_t bH = (b >> 64) & MASK_LO;

  uint128_t result;
  uint128_t hiPart,loPart,midPart;

  // Multiplying (a.h 2^64 + a.l) x (b.h 2^64 + b.l) =
  //			2^128 a.h b.h + 2^64*(a.h b.l+b.h a.l) + a.l b.l
  // get the low part a.l b.l
  // multiply the fractional part
  loPart = aL * bL;
  // compute the middle part 2^64*(a.h b.l+b.h a.l)
  midPart = aL * bH + aH * bL;
  // truncate the low part
  result = (loPart >> 64) + (midPart & MASK_LO);
  // compute the high part 2^128 a.h b.h
  hiPart = aH * bH;
  // truncate the high part and only use the low part
  result |= ((hiPart & MASK_LO) << 64) + (midPart & MASK_HI);
  // if the high part is not zero, put a warning
  NS_ABORT_MSG_IF ((hiPart & MASK_HI) != 0,
		   "High precision 128 bits multiplication error: multiplication overflow.");
  return result;
}
void
int64x64_t::Div (int64x64_t const &o)
{
  bool negResult;
  uint128_t a, b;
  negResult = OUTPUT_SIGN (_v, o._v, a, b);  
  int128_t result = Divu (a, b);
  result = negResult ? -result:result;
  _v = result;
}

uint128_t
int64x64_t::Divu (uint128_t a, uint128_t b)
{
  uint128_t quo = a / b;
  uint128_t rem = (a % b);
  uint128_t result = quo << 64;
  // Now, manage the remainder
  uint128_t tmp = rem >> 64;
  uint128_t div;
  if (tmp == 0)
    {
      rem = rem << 64;
      div = b;
    }
  else
    {
      rem = rem;
      div = b >> 64;
    }
  quo = rem / div;
  result = result + quo;
  return result;
}

void 
int64x64_t::MulByInvert (const int64x64_t &o)
{
  bool negResult = _v < 0;
  uint128_t a = negResult?-_v:_v;
  uint128_t result = UmulByInvert (a, o._v);

  _v = negResult?-result:result;
}
uint128_t
int64x64_t::UmulByInvert (uint128_t a, uint128_t b)
{
  uint128_t result, ah, bh, al, bl;
  uint128_t hi, mid;
  ah = a >> 64;
  bh = b >> 64;
  al = a & MASK_LO;
  bl = b & MASK_LO;
  hi = ah * bh;
  mid = ah * bl + al * bh;
  mid >>= 64;
  result = ah * bh + mid;
  return result;
}
int64x64_t 
int64x64_t::Invert (uint64_t v)
{
  NS_ASSERT (v > 1);
  uint128_t a;
  a = 1;
  a <<= 64;
  int64x64_t result;
  result._v = Divu (a, v);
  int64x64_t tmp = int64x64_t (v, false);
  tmp.MulByInvert (result);
  if (tmp.GetHigh () != 1)
    {
      result._v += 1;
    }
  return result;
}

} // namespace ns3