Modular arithmetic: Difference between revisions
→{{header|C}}
m (missing article) |
|||
Line 80:
1
</pre>
=={{header|ATS}}==
The following program uses the <code>unsigned __int128</code> type of GNU C. It could easily be modified, however, not to use that. I use the compiler extension to implement multiplication that, for all the ordinary integer types on AMD64, does not overflow.
<syntaxhighlight lang="ATS">
(* The program is compiled to C and the integer types have C
semantics. This means ordinary unsigned arithmetic is already
modular! However, the modulus is fixed at 2**n, where n is the
number of bits in the unsigned integer type.
Below, I let a "modulus" of zero mean to use 2**n as the
modulus. *)
(*------------------------------------------------------------------*)
#include "share/atspre_staload.hats"
staload UN = "prelude/SATS/unsafe.sats"
(*------------------------------------------------------------------*)
(* An abstract type, the size of @(g0uint tk, g1uint (tk, modulus)) *)
abst@ype modular_g0uint (tk : tkind, modulus : int) =
@(g0uint tk, g1uint (tk, modulus))
(* Because the type is abstract, we need a constructor: *)
extern fn {tk : tkind}
modular_g0uint_make
{m : int} (* "For any integer m" *)
(a : g0uint tk,
m : g1uint (tk, m))
:<> modular_g0uint (tk, m)
(* A deconstructor: *)
extern fn {tk : tkind}
modular_g0uint_unmake
{m : int}
(a : modular_g0uint (tk, m))
:<> @(g0uint tk, g1uint (tk, m))
extern fn {tk : tkind}
modular_g0uint_succ (* "Successor" *)
{m : int}
(a : modular_g0uint (tk, m))
:<> modular_g0uint (tk, m)
extern fn {tk : tkind} (* This won't be used, but let us write it. *)
modular_g0uint_pred (* "Predecessor" *)
{m : int}
(a : modular_g0uint (tk, m))
:<> modular_g0uint (tk, m)
extern fn {tk : tkind} (* This won't be used, but let us write it.*)
modular_g0uint_neg
{m : int}
(a : modular_g0uint (tk, m))
:<> modular_g0uint (tk, m)
extern fn {tk : tkind}
modular_g0uint_add
{m : int}
(a : modular_g0uint (tk, m),
b : modular_g0uint (tk, m))
:<> modular_g0uint (tk, m)
extern fn {tk : tkind} (* This won't be used, but let us write it.*)
modular_g0uint_sub
{m : int}
(a : modular_g0uint (tk, m),
b : modular_g0uint (tk, m))
:<> modular_g0uint (tk, m)
extern fn {tk : tkind}
modular_g0uint_mul
{m : int}
(a : modular_g0uint (tk, m),
b : modular_g0uint (tk, m))
:<> modular_g0uint (tk, m)
extern fn {tk : tkind}
modular_g0uint_npow
{m : int}
(a : modular_g0uint (tk, m),
i : intGte 0)
:<> modular_g0uint (tk, m)
overload succ with modular_g0uint_succ
overload pred with modular_g0uint_pred
overload ~ with modular_g0uint_neg
overload + with modular_g0uint_add
overload - with modular_g0uint_sub
overload * with modular_g0uint_mul
overload ** with modular_g0uint_npow
(*------------------------------------------------------------------*)
local
(* We make the type be @(g0uint tk, g1uint (tk, modulus)).
The first element is the least residue, the second is the
modulus. A modulus of 0 indicates that the modulus is 2**n, where
n is the number of bits in the typekind. *)
typedef _modular_g0uint (tk : tkind, modulus : int) =
@(g0uint tk, g1uint (tk, modulus))
in (* local *)
assume modular_g0uint (tk, modulus) = _modular_g0uint (tk, modulus)
implement {tk}
modular_g0uint_make (a, m) =
if m = g1i2u 0 then
@(a, m)
else
@(a mod m, m)
implement {tk}
modular_g0uint_unmake a =
a
implement {tk}
modular_g0uint_succ a =
let
val @(a, m) = a
in
if (m = g1i2u 0) || (succ a <> m) then
@(succ a, m)
else
@(g1i2u 0, m)
end
implement {tk}
modular_g0uint_pred a =
let
val @(a, m) = a
prval () = lemma_g1uint_param m
in
(* An exercise for the advanced reader: how come in
modular_g0uint_succ I could use "||", but here I have to use
"+" instead? *)
if (m = g1i2u 0) + (a <> g1i2u 0) then
@(pred a, m)
else
@(pred m, m)
end
implement {tk}
modular_g0uint_neg a =
let
val @(a, m) = a
in
if m = g1i2u 0 then
@(succ (lnot a), m) (* Two's complement. *)
else if a = g0i2u 0 then
@(a, m)
else
@(m - a, m)
end
implement {tk}
modular_g0uint_add (a, b) =
let
(* The modulus of b WILL be same as that of a. The type system
guarantees this at compile time. *)
val @(a, m) = a
and @(b, _) = b
in
if m = g1i2u 0 then
@(a + b, m)
else
@((a + b) mod m, m)
end
implement {tk}
modular_g0uint_mul (a, b) =
(* For multiplication there is a complication, which is that the
product might overflow the register and so end up reduced
modulo the 2**(wordsize). Approaches to that problem are
discussed here:
https://en.wikipedia.org/w/index.php?title=Modular_arithmetic&oldid=1145603919#Example_implementations
However, what I will do is inline some C, and use a GNU C
extension for an integer type that (on AMD64, at least) is
twice as large as uintmax_t.
In so doing, perhaps I help demonstrate how suitable ATS is for
low-level systems programming. Inlining the C is very easy to
do. *)
let
val @(a, m) = a
and @(b, _) = b
in
if m = g1i2u 0 then
@(a * b, m)
else
let
typedef big = $extype"unsigned __int128"
(* A call to _modular_g0uint_mul will actually be a call to
a C function or macro, which happens also to be named
_modular_g0uint_mul. *)
extern fn
_modular_g0uint_mul
(a : big,
b : big,
m : big)
:<> big = "mac#_modular_g0uint_mul"
in
(* The following will work only as long as the C compiler
itself knows how to cast the integer types. There are
safer methods of casting, but, for this task, let us
ignore that. *)
@($UN.cast (_modular_g0uint_mul ($UN.cast a,
$UN.cast b,
$UN.cast m)),
m)
end
end
(* The following puts a static inline function _modular_g0uint_mul
near the top of the C source file. *)
%{^
ATSinline() unsigned __int128
_modular_g0uint_mul (unsigned __int128 a,
unsigned __int128 b,
unsigned __int128 m)
{
return ((a * b) % m);
}
%}
end (* local *)
implement {tk}
modular_g0uint_sub (a, b) =
a + (~b)
implement {tk}
modular_g0uint_npow {m} (a, i) =
(* To compute a power, the multiplication implementation devised
above can be used. The algorithm here is simply the squaring
method:
https://en.wikipedia.org/w/index.php?title=Exponentiation_by_squaring&oldid=1144956501 *)
let
fun
repeat {i : nat} (* <-- This number consistently shrinks. *)
.<i>. (* <-- Proof the recursion will terminate. *)
(accum : modular_g0uint (tk, m), (* "Accumulator" *)
base : modular_g0uint (tk, m),
i : int i)
:<> modular_g0uint (tk, m) =
if i = 0 then
accum
else
let
val i_halved = half i (* Integer division. *)
and base_squared = base * base
in
if i_halved + i_halved = i then
repeat (accum, base_squared, i_halved)
else
repeat (base * accum, base_squared, i_halved)
end
val @(_, m) = modular_g0uint_unmake<tk> a
in
repeat (modular_g0uint_make<tk> (g0i2u 1, m), a, i)
end
(*------------------------------------------------------------------*)
extern fn {tk : tkind}
f : {m : int} modular_g0uint (tk, m) -<> modular_g0uint (tk, m)
(* Using the "successor" function below means that, to add 1, we do
not need to know the modulus. That is why I added "succ". *)
implement {tk}
f(x) = succ (x**100 + x)
implement
main0 () =
let
val x = modular_g0uint_make (10U, 13U)
in
println! ((modular_g0uint_unmake (f(x))).0)
end
(*------------------------------------------------------------------*)
</syntaxhighlight>
{{out}}
<pre>$ patscc -std=gnu2x -g -O2 modular_arithmetic_task.dats && ./a.out
1</pre>
=={{header|C}}==
| |||