This post provides a workout in generic programming using modules & functors.
The program presented here models univariate polynomials over rings based on an exercise in "The Module Language" chapter, of Didier Rémy's book, Using, Understanding and Unraveling the OCaml Lanaguage.
Arithmetics and rings
We begin with a type for modules implementing arithmetic.
module type ARITH = sig
type t
val of_int : int -> t val to_int : t -> int
val of_string : string -> t val to_string : t -> string
val zero : t val one : t
val add : t -> t -> t val sub : t -> t -> t
val mul : t -> t -> t val div : t -> t -> t
val compare : t -> t -> int val equal : t -> t -> bool
end;;
A ring is a set equipped with two binary operations that
generalize the arithmetic operations of addition and
multiplication.
module type RING = sig
type t
type extern_t
val print : t -> unit
val make : extern_t -> t val show : t -> extern_t
val zero : t val one : t
val add : t -> t -> t val mul : t -> t -> t
val equal : t -> t -> bool
end;;
We can build rings over arithmetics with functors. This particular
one fixes the external representation of the elements of the ring
to int.
module Ring (A : ARITH) :
RING with type t = A.t and type extern_t = int =
struct
include A
type extern_t = int
let make = of_int let show = to_int
let print x = print_int (show x)
end;;
Thus, here for example are rings over various specific arithmetic
modules.
module Ring_int32 = Ring (Int32);;
module Ring_int64 = Ring (Int64);;
module Ring_nativeint = Ring (Nativeint);;
module Ring_int = Ring (
struct
type t = int
let of_int x = x let to_int x = x
let of_string = int_of_string let to_string = string_of_int
let zero = 0 let one = 1
let add = ( + ) let sub = ( - )
let mul = ( * ) let div = ( / )
let compare = Pervasives.compare let equal = ( = )
end
);;
Polynomials
We define now the type of polynomials.
module type POLYNOMIAL = sig
type coeff (*Type of coefficients*)
type coeff_extern_t (*Type of coeff. external rep*)
(*Polynomials satisfy the ring interface*)
include RING (*Declares a type [t] and [extern_t]*)
(*Function to evaluate a polynomial at a point*)
val eval : t -> coeff -> coeff
end;;
Given a module implementing a ring, we can generate a module
implementing polynomials with coefficients drawn from the ring.
module Polynomial (R : RING) :
POLYNOMIAL with type coeff = R.t
and type coeff_extern_t = R.extern_t
and type extern_t = (R.extern_t * int) list =
struct
type coeff = R.t (*Coefficient type*)
type coeff_extern_t = R.extern_t (*External coeff. rep*)
type extern_t = (coeff_extern_t * int) list (*External polynomial rep*)
(*List of coefficients and their powers*)
type t = (coeff * int) list (*Invariant : Ordered by powers,
lower order terms at the front*)
(* ... *)
end;;
As the comments indicate, the polynomial data structure is a list
of pairs of coefficients and powers, ordered so that lower powers
come before higher ones. Here's a simple printing utility to aid
visualization.
let print p =
List.iter
(fun (c, k) -> Printf.printf "+ (";
R.print c;
Printf.printf ")X^%d " k)
p
In order that we get a canonical representation, null
coefficients are eliminated. In particular, the null monomial is
simply the empty list.
let zero = []
The multiplicative identity one is not really
necessary as it is just a particular monomial however, its
presence makes polynomials themselves satisfy the interface of
rings.
let one = [R.one, 0]
This helper function constructs monomials.
let monomial (a : coeff) (k : int) =
if k < 0 then
failwith "monomial : negative powers not supported"
else if R.equal a R.zero then [] else [a, k]
Next up, we define addition of polynomials by the following
function. Care is taken to ensure the representation invariant is
respected.
let rec add u v =
match u, v with
| [], _ -> v
| _, [] -> u
| ((c1, k1) :: r1 as p1), ((c2, k2) :: r2 as p2) ->
if k1 < k2 then
(c1, k1) :: (add r1 p2)
else if k1 = k2 then
let c = R.add c1 c2 in
if R.equal c R.zero then add r1 r2
else (c, k1) :: (add r1 r2)
else (c2, k2) :: (add p1 r2)
With monomial and add avaialable, we can
now write make that computes a polynomial from an
external representation. We also give the inverse
function show here too.
let make l =
List.fold_left (fun acc (c, k) ->
add (monomial (R.make c) k) acc) zero l
let show p =
List.fold_right (fun (c, k) acc -> (R.show c, k) :: acc) p []
The module private function times left-multiplies a
polynomial by a monomial.
let rec times (c, k) = function
| [] -> []
| (c1, k1) :: q ->
let c2 = R.mul c c1 in
if R.equal c2 R.zero then times (c, k) q
else (c2, k + k1) :: times (c, k) q
Given the existence of times, polynomial
multiplication can be expressed in a "one-liner".
let mul p = List.fold_left (fun r m -> add r (times m p)) zero
Comparing two polynomials for equality is achieved with the
following predicate.
let rec equal p1 p2 =
match p1, p2 with
| [], [] -> true
| (c1, k1) :: q1, (c2, k2) :: q2 ->
k1 = k2 && R.equal c1 c2 && equal q1 q2
| _ -> false
In the course of evaluating polynomials for a specific value of
their indeterminate, we'll require a function for computing
powers. The following routine uses
the exponentiation
by squaring technique.
let rec pow c = function
| 0 -> R.one
| 1 -> c
| k ->
let l = pow c (k lsr 1) in
let l2 = R.mul l l in
if k land 1 = 0 then l2 else R.mul c l2
Finally, the function eval for evaluation of a
polynomial at a specific point. The implementation
uses Horner's
rule for computationally efficient evaluation.
let eval p c = match List.rev p with
| [] -> R.zero
| (h :: t) ->
let reduce (a, k) (b, l) =
let n = pow c (k - l) in
let t = R.add (R.mul a n) b in
(t, l) in
let a, k = List.fold_left reduce h t in
R.mul (pow c k) a
Testing and example usage
The following interactive session creates a polynomial with integer coefficients module and uses it to confirm the equivalence of $(1 + x)(1 - x)$ with $(1 - x^{2})$.
# #use "polynomial.ml";;
# module R = Ring_int;;
# module P = Polynomial (R);;
# let p = P.mul (P.make [(1, 0); (1, 1)]) (P.make [(1, 0); (-1, 1)]);;
# let q = P.make [(1, 0); (-1, 2)];;
# P.equal p q;;
- : bool = true
Polynomials in two variables, can be treated as univariate
polynomials with polynomial coefficients. For example, the
polynomial $(X + Y)$ can be regarded as $(1 \times X^{1})Y^{0} +
(1 \times X^{0})Y^{1}$. Similarly we can write $X - Y$ as $(1
\times X^{1})Y^{0} + (-1 \times X^{0}) Y^1$ and now check the
equivalence $(X + Y)(X - Y) = (X^{2} - Y^{2})$.
#module Y = Polynomial (P);;
#(* (X + Y) *)
#let r = Y.make [
# ([1, 1], 0); (*(1 X^1) Y^0*)
# ([1, 0], 1) (*(1 X^0) Y^1*)
#];;
#(* (X - Y) *)
#let s = Y.make [
# ([1, 1], 0); (*(1 X^1) Y^0*)
# ([-1, 0], 1) (*((-1) X^0) Y^1*)
#];;
#(* (X^2 - Y^2) *)
#let t = Y.make [
# ([1, 2], 0); (*(1 X^2) Y^0*)
# ([-1, 0], 2) (* (-1 X^0) Y^2*)
#];;
#Y.equal (Y.mul r s) t;;
- : bool = true
