program ModuloMultiplication;
{$mode objfpc}{$H+}
uses
SysUtils;
(*
================================================================
Modulo Multiplication in C (Slow and Fast Versions)
================================================================
Purpose:
Compute (a * b) % mod safely for large 64‑bit values
without using big‑integer libraries.
Why two versions?
1. Slow version:
- Adds 'b' to the result 'a' times.
- Always correct.
- Very slow for large numbers.
- Useful as a correctness reference.
2. Fast version:
- Uses the "double‑and‑add" technique.
- Runs in O(log b).
- Uses unsigned __int128 internally to avoid overflow.
- Produces the same result as the slow version.
================================================================
*)
(* ---------------------------------------------------------------
SLOW VERSION (simple, correct, but extremely slow)
--------------------------------------------------------------- *)
function MulModSlow(a, b, mod_val: QWord): QWord;
var
tmp, result_val, i: QWord;
begin
(*
Computes:
(b + b + b + ... a times) % mod
This avoids overflow because:
- result stays below mod
- b fits in uint64_t
- addition is safe
But it is O(a), which is too slow for large inputs.
*)
if b < a then
begin
tmp := a;
a := b;
b := tmp;
end; // reduce loop count
result_val := 0;
for i := 0 to a - 1 do
begin
result_val := result_val + b;
result_val := result_val mod mod_val;
end;
Exit(result_val);
end;
(* ---------------------------------------------------------------
FAST VERSION (Safe 64-bit double-and-add without UInt128)
--------------------------------------------------------------- *)
function MulModFast(a, b, mod_val: QWord): QWord;
var
result_val: QWord;
begin
(*
Uses the classic "double‑and‑add" method:
- If the lowest bit of b is set, add 'a' to result.
- Double 'a' each step.
- Shift b right each step.
The key detail:
We use double-and-add on 64-bit numbers directly.
To avoid (a + a) overflow during doubling, we take
(a mod mod_val) and subtract mod_val when (a + a) exceeds it.
*)
result_val := 0;
a := a mod mod_val;
b := b mod mod_val;
while b > 0 do
begin
if (b and 1) <> 0 then
begin
if mod_val - result_val > a then
result_val := result_val + a
else
result_val := result_val + a - mod_val;
end;
if mod_val - a > a then
a := a + a
else
a := a + a - mod_val;
b := b shr 1;
end;
Exit(result_val);
end;
(* ---------------------------------------------------------------
MAIN PROGRAM
--------------------------------------------------------------- *)
var
x, y, mod_val: QWord;
slow_result, fast_result: QWord;
begin
x := 798345;
y := 20289473612815;
mod_val := 100000000000003;
slow_result := MulModSlow(x, y, mod_val);
fast_result := MulModFast(x, y, mod_val);
WriteLn(Format('Slow result: %d', [slow_result]));
WriteLn(Format('Fast result: %d', [fast_result]));
end.
(*
run:
Slow result: 99811422305238
Fast result: 99811422305238
*)
