## Tuesday, March 21, 2017

### Polynomials over rings

Polynomials over rings

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