void ModExp(uint1024_t M, uint1024_t e, uint1024_t n,
            uint1024_t Mbar, uint1024_t xbar, uint1024_t* out)
{
    uint1024_t xbar_temp = 0;  // temporary intermediate variable

    for (int i=NUM_BITS-1; i>=0; i--)
    {
        // square
        montMult(xbar,xbar,n,&xbar);

        // multiply if odd
        if (e.test(i)) // if (e.bit(i) == 1), then e is odd
        {
            montMult(Mbar,xbar,n, &xbar_temp);
            xbar = xbar_temp; // update temporary variable
        }
    }

    // undo montgomery residue transformation
    montMult(xbar_temp,1,n, out);
}