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 Trauriges Smiley ).

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

    /// addition
    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 add a b = ops.Add(a,b)
        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)
            ] |> List.fold add ops.Zero
        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