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