/* 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; }