I would like to try to implement Yun's algorithm for square-free factorization of polynomials. From Wikipedia (f is the polynomial):
a0 = gcd(f, f'); b1 = f/a0; c1 = f'/a0; d1 = c1 - b1'; i = 1
repeat
ai = gcd(bi, di); bi+1 = bi/ai; ci+1 = di/ai; i = i + 1; di = ci - bi'
until b = 1
However, I'm not sure about the second step. I would like to use it for polynomials with integer coefficients (not necessary monic or primitive). Is it possible to realize the division b1 = f/a0 using just integers?
I found the code for synthetic division:
def extended_synthetic_division(dividend, divisor):
'''Fast polynomial division by using Extended Synthetic Division. Also works with non-monic polynomials.'''
# dividend and divisor are both polynomials, which are here simply lists of coefficients. Eg: x^2 + 3x + 5 will be represented as [1, 3, 5]
out = list(dividend) # Copy the dividend
normalizer = divisor[0]
for i in xrange(len(dividend)-(len(divisor)-1)):
out[i] /= normalizer # for general polynomial division (when polynomials are non-monic),
# we need to normalize by dividing the coefficient with the divisor's first coefficient
coef = out[i]
if coef != 0: # useless to multiply if coef is 0
for j in xrange(1, len(divisor)): # in synthetic division, we always skip the first coefficient of the divisor,
# because it is only used to normalize the dividend coefficients
out[i + j] += -divisor[j] * coef
# The resulting out contains both the quotient and the remainder, the remainder being the size of the divisor (the remainder
# has necessarily the same degree as the divisor since it is what we couldn't divide from the dividend), so we compute the index
# where this separation is, and return the quotient and remainder.
separator = -(len(divisor)-1)
return out[:separator], out[separator:] # return quotient, remainder.
The problem for me is that out[i] /= normalizer. Would it always work with integer (floor) division for Yun's b1 = f/a0? Is it so that it is always possible to divide f/gcd(f, f')? Is the out[separator:] (remainder) always going to zero?