66 lines
1.9 KiB
Plaintext
66 lines
1.9 KiB
Plaintext
(* You will need
|
|
https://sourceforge.net/p/chemoelectric/ats2-xprelude/ *)
|
|
|
|
#include "share/atspre_staload.hats"
|
|
|
|
#include "xprelude/HATS/xprelude.hats"
|
|
staload "xprelude/SATS/exrat.sats"
|
|
staload _ = "xprelude/DATS/exrat.dats"
|
|
|
|
val a = exrat_make_string_exn "2988348162058574136915891421498819466320163312926952423791023078876139"
|
|
val b = exrat_make_string_exn "2351399303373464486466122544523690094744975233415544072992656881240319"
|
|
|
|
val modulus = exrat_make (10, 1) ** 40
|
|
|
|
(* xprelude/SATS/exrat.sats includes the "exrat_numerator_modular_pow"
|
|
function, based on GMP's mpz_powm. *)
|
|
val result1 = exrat_numerator_modular_pow (a, b, modulus)
|
|
|
|
(* But that was too easy. Here is the right-to-left binary method,
|
|
https://en.wikipedia.org/w/index.php?title=Modular_exponentiation&oldid=1136216610#Right-to-left_binary_method
|
|
*)
|
|
val result2 =
|
|
(lam (base : exrat,
|
|
exponent : exrat,
|
|
modulus : exrat) : exrat =>
|
|
let
|
|
val zero = exrat_make (0, 1)
|
|
and one = exrat_make (1, 1)
|
|
and two = exrat_make (2, 1)
|
|
macdef divrem = exrat_numerator_euclid_division
|
|
macdef rem = exrat_numerator_euclid_remainder
|
|
in
|
|
if modulus = one then
|
|
zero
|
|
else
|
|
let
|
|
fun
|
|
loop (result : exrat,
|
|
base : exrat,
|
|
exponent : exrat) : exrat =
|
|
if iseqz exponent then
|
|
result
|
|
else
|
|
let
|
|
val @(exponent, remainder) = exponent \divrem two
|
|
val result =
|
|
if remainder = one then
|
|
(result * base) \rem modulus
|
|
else
|
|
result
|
|
val base = (base * base) \rem modulus
|
|
in
|
|
loop (result, base, exponent)
|
|
end
|
|
in
|
|
loop (one, base \rem modulus, exponent)
|
|
end
|
|
end) (a, b, modulus)
|
|
|
|
implement
|
|
main0 () =
|
|
begin
|
|
println! result1;
|
|
println! result2
|
|
end
|