codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
/* 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; }
Private
[
?
]
Run code
Submit