poly_gcd(p,q)

In the longhand polynomial division given as
P(k) = P(k-2) - P(k-1)*Q(k)
The quotient Q(k) and the remainder P(k) are obtained from dividing the dividend P(k-2) by the divisor P(k-1). If we can make Q(k) = 1, by converting P(k-2) and P(k-1) into equal degree and monic, then the longhand polynomial division becomes simply the "monic polynomial subtraction" (MPS):
P(k) = P(k-2) - P(k-1)
For a pair of given polynomials p(x) and q(x) of degree n and m, n>m, we set
P(1) = p(x)/p_0
P(2) = q(x)*x^(n-m)/q_0
Applying the MPS repeatedly starting from k=3, until k=K+1, such that
P(K+1) = P(k-1) - P(k) = 0
then we get our desired polynomial GCD as
gcd(p,q) = P(K).
The source code uses only basic MATLAB built-in functions. Its listing is only 17 lines total !
Amazingly, this simple routine gives the expected results for the test polynomials and their derivatives of very high degree, such as
p(x) = (x + 1)^1000
p(x) = (x + 123456789)^30
p(x) = (1234x + 56789)^60
p(x) = (x^4-2x^3+3x^2-4x +5)^50
p(x) = (x^4 - 1)^25
*************** UPDATE (10/05/09): **************
The approach "Leading-coefficient Elinimation" is revised from the original "Monic Polynomial Subtraction".
It also reduces almost half of the total arithematic operations.
The total source code listing is only 12 lines!
*************** UPDATE (01/22/2018): **************
The source code function g = poly_gcd(p,q) is revised and updated. It greatly reduces the overall operation procedures.
Please see the typical examples in the comment section.