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 (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?

CodePudding user response：

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]

CodePudding user response：

The Rosetta code assumes that [a, b, c, ...] represents a bx cx^2 .... But you are assuming it should represent ax^n bx^[n-1] cx^(n-2) ....

There may be other bugs, but that one leapt out at me.