/* montgomery multiplication */

typedef unsigned long long ull;
typedef signed long long sll;

void inverse(ull r, ull m, ull *r_inv, ull *m_prime) {
    ull r_saved, m_saved, u, v;
    r_saved = r; m_saved = m; u = 1; v = 0;
    while (r > 0) {
        r = r >> 1;
        if ((u & 1) == 0) {
            u = u >> 1; v = v >> 1; }
        else {
            u = ((u ^ m_saved) >> 1) + (u & m_saved);
            v = (v >> 1) + r_saved; } }
    *r_inv = u; *m_prime = v;
    return; }

void mul_ull(ull x, ull y, ull *xy_hi, ull *xy_lo) {
    ull x0, x1, y0, y1, t, p0, p1, p2;
    x1 = x >> 32; x0 = x & 0xFFFFFFFF;
    y1 = y >> 32; y0 = y & 0xFFFFFFFF;

    t = x0 * y0;
    p0 = t & 0xFFFFFFFF;

    t = x1 * y0 + (t >> 32);
    p1 = t & 0xFFFFFFFF;
    p2 = t >> 32;

    t = x0 * y1 + p1;

    *xy_lo = (t << 32) + p0;
    *xy_hi = x1 * y1 + p2 + (t >> 32);

    return; }

ull mod_ull(ull x, ull y, ull m) {
    sll i, t;
    for (i = 1; i <= 64; i++) {
        t = (sll) x >> 63;
        x = (x << 1) | (y >> 63);
        y = y << 1;
        if ((x | t) >= m) {
            x = x - m; y = y + 1; } }
    return x; }

ull mont_mul(ull a_bar, ull b_bar, ull m, ull m_prime) {
    ull t_hi, t_lo, tm, tmm_hi, tmm_lo, u_hi, u_lo, ov;
    mul_ull(a_bar, b_bar, &t_hi, &t_lo);
    tm = t_lo * m_prime;
    mul_ull(tm, m, &tmm_hi, &tmm_lo);
    u_lo = t_lo + tmm_lo; u_hi = t_hi + tmm_hi;
    if (u_lo < t_lo) u_hi = u_hi + 1;
    ov = (u_hi < t_hi) | ((u_hi == t_hi) & (u_lo < t_lo));
    u_lo = u_hi; u_hi = 0;
    u_lo = u_lo - (m & -(ov | (u_lo >= m)));
    return u_lo; }

int mulmod_ull(ull a, ull b, ull m, ull *p) {
    ull half_r = 0x8000000000000000LL;
    ull r_inv, m_prime, a_bar, b_bar, p_hi, p_lo;

    if ((m & 1) == 0) { return 0; } /* must be odd */
    if (m <= a) { return 0; } /* must be less than m */
    if (m <= b) { return 0; } /* must be less than m */

    inverse(half_r, m, &r_inv, &m_prime);
    a_bar = mod_ull(a, 0, m); b_bar = mod_ull(b, 0, m);
    *p = mont_mul(a_bar, b_bar, m, m_prime);
    mul_ull(*p, r_inv, &p_hi, &p_lo);
    *p = mod_ull(p_hi, p_lo, m);
    return 1; }

int powmod_ull(ull x, ull y, ull m, ull *x_to_y) {
    ull half_r = 0x8000000000000000LL;
    ull r_inv, m_prime, x_bar, hi, lo;

    if ((m & 1) == 0) { return 0; } /* must be odd */
    if (m <= x) { return 0; } /* must be less than m */

    inverse(half_r, m, &r_inv, &m_prime);
    x_bar = mod_ull(x, 0, m);
    *x_to_y = mod_ull(1, 0, m);
    while (y > 0) {
        if ((y & 1) == 1) {
            *x_to_y = mont_mul(*x_to_y, x_bar, m, m_prime);
            y = y - 1; }
        else {
            x_bar = mont_mul(x_bar, x_bar, m, m_prime);
            y = y / 2; } }
    mul_ull(*x_to_y, r_inv, &hi, &lo);
    *x_to_y = mod_ull(hi, lo, m);
    return 1; }

#include <stdio.h>

int main(void) {
    ull z; int t;
    ull x = 34721908534901LL;
    ull y = 72193687003295LL;
    ull m = 9412345678901731LL;

    t = mulmod_ull(x, y, m, &z);
    if (t) { printf("%lld\n", z); }
    else { printf("ERROR"); }

    t = powmod_ull(x, y, m, &z);
    if (t) { printf("%lld\n", z); }
    else { printf("ERROR"); }

    return 0; }

3751384291706939
7001634529421238