# a generic polynomial type in F#

In this short article I want to show how to implement a simple polynomial type in F# over a generic field. For this I will use the technique introduced in this great article: Writing generic numeric code

So you will have to get and reference the F#-Powerpack (Fsharp.PowerPack.dll).

I choose to represent the polynomial as a list of it’s coefficients so that

[a0; a1; a2; a3; ...]


will be mapped to the Polynomial $a_{0} + a_{1} X + a_{2} X^{2} + a_{3} X^{3} + ..$.

I implemented the basic operations (+,-,*) and division with remainder using the algorithm I’m sure we all learned in school (although I think I have seen some reports that there are debates of replacing this with more “calculator-work” in the US ).

I thinks it’s all so simple that the code tells all there is to say, but feel free to ask if there are some questions – so without further due: here is the code:

namespace Math.Polynoms

open System

[<AutoOpen>]
module Helpers =

/// zips-two list using an operator
/// if one list is longer it's elements
/// are pasted at the end
/// so [a1;a2;a3;a4] zipWith [b1;b2] becomse
/// [op a1 b1; op a2 b2; a3; a4]
let zipWith op l l' =
let rec inner l l' cont =
match l,l' with
| a::l, b::l' -> inner l l' (fun r -> op (a, b)::r |> cont)
| a::l, []    -> inner l l' (fun r -> a::r |> cont)
| [], b::l'   -> inner l l' (fun r -> b::r |> cont)
| [], []      -> cont []
inner l l' (fun r -> r)

/// removes pending zeroes at the end of a list
let removePendingZeros eqZero l =
let rec inner l zeroAcc cont =
match l with
| []   -> cont []
| a::l -> match eqZero a with
| false -> inner l (fun r -> r) (fun r -> (cont << zeroAcc) (a::r))
| true  -> inner l (fun r -> zeroAcc (a::r)) (fun r -> cont r)
inner l (fun r -> r) (fun r -> r)

/// polynom-type over a generic field 'a
type Polynom<'a>(coeffs : 'a list, ops : IFractional<'a>) =

// remove pending zeros, as ... aX^n + 0X^(n+1) = ... aX^n
let coeffizients = coeffs |> removePendingZeros (fun x -> ops.Equals(ops.Zero, x))

member p.Coeffizients = coeffizients
member p.Operations = ops

// create a runtime
static member CreateRuntime(coeffs) =
let ops = GlobalAssociations.GetNumericAssociation<'a>() :?> IFractional<'a>
new Polynom<_>(coeffs, ops)

// write the polynom nicely
override p.ToString() =
if coeffizients.Length = 0 then "0" else
let showMonom i c =
let x = if i = 0 then "" elif i = 1 then "X" else sprintf "X^%d" i
if ops.Equals (c, ops.Zero) then "" else sprintf "%A %s" c x
let conc s s' =
if System.String.IsNullOrEmpty s then s'
else s + " + " + s'
coeffizients
|> List.mapi (fun i c -> showMonom i c)
|> List.filter (not << System.String.IsNullOrEmpty)
|> List.rev
|> List.fold conc ""

/// Zero = 0*X^0
static member Zero : Polynom<'a> =
let ops = GlobalAssociations.GetNumericAssociation<'a>() :?> IFractional<'a>
new Polynom<_>([ ops.Zero ], ops)

/// One = 1*X^0
static member One : Polynom<'a> =
let ops = GlobalAssociations.GetNumericAssociation<'a>() :?> IFractional<'a>
new Polynom<_>([ ops.One ], ops)

/// equals
member a.Equals (b : Polynom<'a>) : bool =
if a.Coeffizients.Length <> b.Coeffizients.Length then false else
let ops = a.Operations
List.zip a.Coeffizients b.Coeffizients
|> List.forall ops.Equals

override a.Equals (b : obj) : bool =
match b with
| 😕 Polynom<'a> as b -> a.Equals b
| _                   -> false

override p.GetHashCode() = (coeffizients :> obj).GetHashCode()

static member (==) (a,b) = a.Equals b
static member (!=) (a,b) = a = b |> not

static member (+) (a : Polynom<'a>, b : Polynom<'a>) : Polynom<'a> =
let coeffs = zipWith (a.Operations.Add) a.Coeffizients b.Coeffizients
new Polynom<_>(coeffs, a.Operations)

/// substraction
static member (-) (a : Polynom<'a>, b : Polynom<'a>) : Polynom<'a> =
let coeffs = zipWith (a.Operations.Subtract) a.Coeffizients b.Coeffizients
new Polynom<_>(coeffs, a.Operations)

/// multiplication
static member (*) (a : Polynom<'a>, b : Polynom<'a>) : Polynom<'a> =
let ops = a.Operations
let maxPot = a.Coeffizients.Length + b.Coeffizients.Length - 2
let getCoeff i (p : Polynom<'a>) =
if i >= p.Coeffizients.Length then ops.Zero else p.Coeffizients.[i]
let getProdCoeff k =
[ for i in 0..k do
let j = k-i
yield ops.Multiply (a |> getCoeff i, b |> getCoeff j)
let coeffs = [0 .. maxPot] |> List.map getProdCoeff
new Polynom<_>(coeffs, ops)

/// DivideWithReminder (P1,P2) -> (D,R) so that P1 = D*P2 + R
static member DivideWithReminder (a : Polynom<'a>, b : Polynom<'a>) =
let ops = a.Operations
let divMonoms (a, i) (b, j) =
let k = i-j
if k < 0 then None
else (ops.Divide(a,b), k) |> Some
let getMonom (p : Polynom<_>) =
let i = p.Coeffizients.Length - 1
let a = p.Coeffizients.[i]
(a,i)
let multByMonom (p : Polynom<_>) (c,k) =
let rec prepZeroes k acc =
if k <= 0 then acc
else prepZeroes (k-1) (ops.Zero::acc)
let coeffs =
p.Coeffizients
|> List.map (fun a -> ops.Multiply(a,c))
|> prepZeroes k
new Polynom<_>(coeffs, ops)
let one = new Polynom<_>([ ops.One ], ops)
let zero = new Polynom<_>([ ops.Zero ], ops)
let rec div (a : Polynom<_>) (b : Polynom<_>) acc =
if a.Coeffizients.Length = 0 then acc, zero else
match divMonoms (getMonom a) (getMonom b) with
| Some (c, k) -> let rest = a - (multByMonom b (c, k))
let part = multByMonom one (c,k)
div rest b (acc + part)
| None        -> acc, a
let d, rest = div a b zero
d, rest


Here are some samples demonstrating the type – please not that you will have to reference FSharp.Powerpack.dll if you use F#-interactive first with:

#r FSharp.Powerpack.dll;;

module Sample =
let p1 = Polynom.CreateRuntime([1.0; 2.0; 0.5; 1.0]);
let p2 = Polynom.CreateRuntime([0.5; 1.0]);

let d,r = Polynom.DivideWithReminder(p1,p2);
let p1' = d*p2 + r;
System.Diagnostics.Debug.Assert( (p1 = p1') )

let p3 = Polynom.CreateRuntime([0.0; 0.5; 2.0]);
let d',r' = Polynom.DivideWithReminder(p1,p3);
let p1'' = d'*p3 + r';
System.Diagnostics.Debug.Assert( (p1 = p1'') )


Here is the output (sorry I missed this somehow):

   val p1 : Polynom<float> = 1.0 X^3 + 0.5 X^2 + 2.0 X + 1.0
val p2 : Polynom<float> = 1.0 X + 0.5
val r : Polynom<float> = 0
val d : Polynom<float> = 1.0 X^2 + 2.0
val p1' : Polynom<float> = 1.0 X^3 + 0.5 X^2 + 2.0 X + 1.0


and

   val p3 : Polynom<float> = 2.0 X^2 + 0.5 X
val r' : Polynom<float> = 1.9375 X + 1.0
val d' : Polynom<float> = 0.5 X + 0.125
val p1'' : Polynom<float> = 1.0 X^3 + 0.5 X^2 + 2.0 X + 1.0