Barrett Reduction - Possible without overflow and floating point arithmetic? I have two variables, $x$ and $y$. I'm looking to compute the product of $x$ and $y$ modulo $p$, where $p$ is prime.
For small values of $x$ and $y$ (less than 43 bits each), I was able to perform the computation using Barrett Reduction.
However, for large values of $x$ and $y$, overflow exceeding 128 bits occurs, and the computed value is incorrect. This is because $x$ and $y$ must be multiplied together, which must be stored in a 128-bit datatype. That product must then be multiplied by a precomputed constant, $r$, resulting in a value that can be stored in 192 bits.
Is it possible to compute xy mod p using Barrett Reduction where overflow over 128 bits does not occur, also while not relying on floating point arithmetic?
If it helps, I have the ability to perform the multiplication of $x$ and $y$ and store the product in two variables: $a$ and $b$, where $a$ holds the high 64 bits of the product, and $b$ holds the low 64 bits of the product.
 A: Yes, but in addition to $r = \left\lfloor \dfrac{2^{64}}p\right\rfloor$, you will also need $$s =\left\lfloor \dfrac{2^{128}}p\right\rfloor - 2^{64}r$$
and $$t = 2^{64} - rp$$
Note that $s,t < 2^{64}$.
So you have $xy = a2^{64} + b$. Then 
$$\begin{align}xy \bmod p &= (a2^{64} + b) \bmod p\\&= \big[(a2^{64}\bmod p) + (b\bmod p)\big]\bmod p\end{align}$$
Per Barrett's algorithm, to calculate $a2^{64}\bmod p$, you first calculate
$$q = \left\lfloor\dfrac {a2^{64}(2^{24}r + s)}{2^{128}}\right\rfloor = ar + \left\lfloor\frac{as}{2^{64}}\right\rfloor$$
The approximate modulus is then
$$2^{64}a - qp = a(2^{64} - rp) - \left\lfloor\frac{as}{2^{64}}\right\rfloor p= at - \left\lfloor\frac{as}{2^{64}}\right\rfloor p$$
In pseudocode, to find $x \bmod p$ for a 128-bit $x$
function Mod128(x, p)
   a = x >> 64
   b = x - (a << 64)
   qa = a * s >> 64
   qb = b * r >> 64
   a1 = a * t - qa * p
   if a1 >= p then a1 = a1 - p
   b1 = b - qb * p
   if b1 >= p then b1 = b1 - p
   x1 = a1 + b1
   if x1 >= p then x1 = x1 - p
   return x1

Variation assuming two functions hi64(x,y) and lo64(x,y) that return the high and low 64 bits of the product of 64 bit integers x and y. If $p <2^{63}$, this version is guaranteed to never overflow using unsigned 64 bit arithmetic, since qb * p can be at most b:
function ModProd(x,y,p)
   a = hi64(x,y)
   b = lo64(x,y)
   qa = hi64(a,s)
   qb = hi64(b,r)
   at = lo64(a,t)
   qap = lo64(qa,p)
   if at < qap then 
      a1 = (2^63 - qap) + 2^63 + at
   else
      a1 = at - qap
   if a1 >= p then a1 = a1 - p
   b1 = b - qb * p
   if b1 >= p then b1 = b1 - p
   prod = a1 + b1
   if prod >= p then prod = prod - p
   return prod

