We can easily calculate the greatest common divisor of two numbers using the following classical method:
let rec gcd a b =
if b = 0 then a else gcd b (a % b)
This method can be extended to univariate polynomials if we give a way to obtain the remainder of the division of two polynomials, using long division. Assuming rem a b is the remainder of polynomials a and b, then we find their GCD in the same way. We can represent a polynomial of degree n as its coefficient list, a list of size n. Then we have the following:
let rec polynomialGCD a b =
if b = zeroPolynomial then a else polynomialGCD b (rem a b)
Where zeroPolynomial therefore represents an empty list. To define rem, we need the long division. An algorithm proposed on the RosettaCode website gives the following code:
let private longDivision (xs: float list) (ys: float list) =
let rec shift n l = if n <= 0 then l else shift (n-1) (l @ [0.0])
let rec pad n l = if n <= 0 then l else pad (n-1) (0.0 :: l)
let rec norm = function | 0.0 :: tl -> norm tl | x -> x
let deg l = List.length (norm l) - 1
let zip op p q =
let d = (List.length p) - (List.length q) in
List.map2 op (pad (-d) p) (pad d q)
let polydiv f g =
let rec aux f s q =
let ddif = (deg f) - (deg s) in
if ddif < 0 then (q, f) else
let k = (List.head f) / (List.head s) in
let ks = List.map (( * ) k) (shift ddif s) in
let q' = zip (+) q (shift ddif [k])
let f' = norm (List.tail (zip (-) f ks)) in
aux f' s q' in
aux (norm f) (norm g) []
polydiv xs ys
The function longDivision returns the tuple (quotient, remainder), so we can define our rem a b as follows:
let rem a b = snd (longDivision a b)
And finally, I recall the GDC code of two polynomials, we have:
let rec polynomialGCD a b =
if b = [] then a else polynomialGCD b (rem a b)
However, it doesn't work. I tried with some examples and the result is always bad, see below:
let poly1 = [ 1.; 0.; 6.; 7. ] // x³ + 6x + 7
let poly2 = [ 1.; 3.; 2. ] // x² + 3x + 2
let gcd = polynomialGCD poly1 poly2 // Normally, the GCD is x + 1, i.e. [1; 1]
let poly1' = [ 1.; 2.; 0.; -1.;0 ; 1. ] // x⁵ + 2x⁴ - x² + 1
let poly2' = [ 1.; 0.; 0.; 0.; -1. ] // x⁴ - 1
let gcd' = polynomialGCD poly1' poly2' // Normally, the GCD is 1, i.e. [1]
I obtain gcd = [13.0; 13.0] and gdc' = [-1.040816327], I could conceive that 13x + 13 still factorizes to 13(x+1), although the answer to be obtained would normally be x+1 (so [1; 1]) but the result of gcd' is completely wrong. I don't really understand what is wrong with my code, what do I need to change to make it work?
I think your code is actually correct, but it's missing a final step. If you want the GCD to be monic (that is, to have 1 as its coefficient of the highest degree), you have to multiply by a constant at the end.
So, as you noticed, in your first example, 13x + 13 can be multiplied by 1/13 to get x + 1, which is monic (i.e. 1x^1 + 1x^0).
Similarly, in your second example, -51/49 can be multiplied by -49/51 to get 1, which is the monic result you expect.
Here's a simple fix that implements this step:
let polynomialGCD a b =
let rec loop a b =
if b = [] then a
else loop b (rem a b)
let raw = loop a b
let q, r = polydiv raw [List.head raw] // force a monic result
assert(List.head q = 1.0)
assert(r = [])
q
Testing your examples:
let gcd = polynomialGCD poly1 poly2
let gcd' = polynomialGCD poly1' poly2'
printfn "%A" gcd // output: [1.0; 1.0]
printfn "%A" gcd' // output: [1.0]