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 
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
/// 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
