multmod.ml

(* a + b (mod m) *)
let addmod a b m =
  assert (0 <= a && 0 <= b);
  let a = a mod m in
  let b = b mod m in
  if a < m - b then
    (* Case A:  a + b < m *)
    a + b
  else
    (* Case B:  a + b >= m

            a >= m - b

            a + b ≡ a - (m - b)
      *)
    a - (m - b)

(* a * b (mod m) *)
let multmod a b m =
  let a = a mod m in
  let b = b mod m in
  let rec loop a b =
    if a > b then loop b a
    else if a = 0 then 0
    else if a mod 2 = 1 then addmod (loop (a-1) b) b m
    else
      let res = loop (a/2) b in
      addmod res res m
  in
  loop a b

let () =
  Random.self_init ();
  for i = 0 to 100_000_000 do
    if i mod 1000000 = 0 then Format.eprintf "%d@." i;
    let a = Random.int (1 lsl 30 - 1) in
    let b = Random.int (1 lsl 30 - 1) in
    let c = Random.int (1 lsl 30 - 1) in
    let r = (a * b) mod c in
    assert (r = multmod a b c)
  done