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