def test_polydiv_type(self):
     # Make polydiv work for complex types
     msg = "Wrong type, should be complex"
     x = np.ones(3, dtype=complex)
     q, r = np.polydiv(x, x)
     assert_(q.dtype == complex, msg)
     msg = "Wrong type, should be float"
     x = np.ones(3, dtype=int)
     q, r = np.polydiv(x, x)
     assert_(q.dtype == float, msg)

 def test_polydiv(self):
     b = np.poly1d([2, 6, 6, 1])
     a = np.poly1d([-1j, (1+2j), -(2+1j), 1])
     q, r = np.polydiv(b, a)
     assert_equal(q.coeffs.dtype, np.complex128)
     assert_equal(r.coeffs.dtype, np.complex128)
     assert_equal(q*a + r, b)

    def __init__(self, poly, word_size, offset_words):
        self.word_size = word_size
        self.poly = numpy.array(poly, dtype=int)
        self.offset_words = offset_words
        
        # calculate the P matrix by polynomial division
        # each row is: e(i)*x^10 mod rds_poly
        # where e(i) is the i-th base vector in the canonical orthogonal base
        self.check_size = self.poly.size - 1
        self.matP = numpy.empty([0, self.check_size], dtype=int)
        for i in range(word_size):
            (q, r) = numpy.polydiv(numpy.identity(self.word_size+self.check_size, dtype=int)[i], self.poly)
            #print q, r
            # r may be "left-trimmed" => add missing zeros
            if self.check_size - r.size > 0:
                #print r
                #print numpy.zeros(check_size - r.size)
                r = numpy.append(numpy.zeros(self.check_size - r.size, dtype=int), r)

            rr = numpy.mod(numpy.array([r], dtype=int), 2)
            self.matP = numpy.append(self.matP, rr, axis=0)
            
        self.matG = numpy.append(numpy.identity(self.word_size, dtype=int), self.matP, axis=1)
        self.matH = numpy.append(self.matP, numpy.identity(self.check_size, dtype=int), axis=0)
        
        #self.offset_words = numpy.array(offset_words, dtype=int)
        self.syndromes = {}
        for ow_name, ow in offset_words.items():
            # actually it's useless to call syndrome here, because of the way
            # our H is constructed. Do the block-wise matrix multiplication
            # to be convinced of this.
            self.syndromes[ow_name] = self.syndrome(numpy.append(numpy.zeros(self.word_size, dtype=int), numpy.array(ow, dtype=int)))

 def __truediv__(self, other):
     #return self.__div(other)
     if sc.isscalar(other):
         return poly1d(self.coeffs/other, variable='s')
     else:
         other = sc.poly1d(other, variable='s')
         tplAt =  sc.polydiv(self, other)
         return (poly1d(tplAt[0]), poly1d(tplAt[1]))

 def adjust(self, u):
     assert self.poly[0] == 1
     u = np.array(u) % self.p
     div,mod = np.polydiv(u, self.poly)
     result = np.array(mod).astype(int) % self.p
     if len(result) == 0:
         result = np.array([0])
     return result

 def __rtruediv__(self, other):
     #return "XXX"
     #return self.__div__(a)
     if sc.isscalar(other):
         return sc.poly1d(other/self.coeffs, variable='s')
     else:
         other = sc.poly1d(other, variable='s')
         tplAt =  sc.polydiv(other, self)
         return (poly1d(tplAt[0]), poly1d(tplAt[1]))

def residue(b,a,tol=1e-3,rtype='avg'):
    """Compute partial-fraction expansion of b(s) / a(s).

    If M = len(b) and N = len(a)

            b(s)     b[0] s**(M-1) + b[1] s**(M-2) + ... + b[M-1]
    H(s) = ------ = ----------------------------------------------
            a(s)     a[0] s**(N-1) + a[1] s**(N-2) + ... + a[N-1]

             r[0]       r[1]             r[-1]
         = -------- + -------- + ... + --------- + k(s)
           (s-p[0])   (s-p[1])         (s-p[-1])

    If there are any repeated roots (closer than tol), then the partial
    fraction expansion has terms like

            r[i]      r[i+1]              r[i+n-1]
          -------- + ----------- + ... + -----------
          (s-p[i])  (s-p[i])**2          (s-p[i])**n

    See also:  invres, poly, polyval, unique_roots
    """

    b,a = map(asarray,(b,a))
    k,b = polydiv(b,a)
    p = roots(a)
    r = p*0.0
    pout, mult = unique_roots(p,tol=tol,rtype=rtype)
    p = []
    for n in range(len(pout)):
        p.extend([pout[n]]*mult[n])
    p = asarray(p)
    # Compute the residue from the general formula
    indx = 0
    for n in range(len(pout)):
        bn = b.copy()
        pn = []
        for l in range(len(pout)):
            if l != n:
                pn.extend([pout[l]]*mult[l])
        an = atleast_1d(poly(pn))
        # bn(s) / an(s) is (s-po[n])**Nn * b(s) / a(s) where Nn is
        # multiplicity of pole at po[n]
        sig = mult[n]
        for m in range(sig,0,-1):
            if sig > m:
                # compute next derivative of bn(s) / an(s)
                term1 = polymul(polyder(bn,1),an)
                term2 = polymul(bn,polyder(an,1))
                bn = polysub(term1,term2)
                an = polymul(an,an)
            r[indx+m-1] = polyval(bn,pout[n]) / polyval(an,pout[n]) \
                          / factorial(sig-m)
        indx += sig
    return r, p, k

def div_inc_pow(num, den, order):
    rnum =  np.zeros(len(den))
    for i in range(0,len(num)): rnum[i] = num[-i-1]
    rden = den[::-1]
    res = np.zeros(order)
    for i in range(0, order):
        quot, rem = np.polydiv(rnum, rden)
        res[i], rnum = quot, np.zeros(len(den))
        for i in range(0,len(rem)):
            rnum[i] = rem[i]
    return res[::-1]

def Bairstow(P):
    k = 0
    # compteur d'iteration;
    A = P.coeffs
    B = 1
    C = -8
    # B = A[1]/A[0];
    # C = A[2]/A[0];
    epsilon = 10 ** -10
    V = np.zeros(2)
    while (abs(P(V[1])) > epsilon) & (abs(P(V[0])) > epsilon):
        U = np.array([B, C])
        T = np.poly1d([1, B, C])

        Div = np.polydiv(P, T)  # 1ere div
        Reste = Div[1]
        Q = Div[0]
        R = Reste[0]
        S = Reste[1]

        Div = np.polydiv(Q, T)
        # 2nde div
        Reste = Div[1]
        G = Div[0]
        Rc = -Reste[0]
        Sc = -Reste[1]

        Rb = -B * Rc + Sc
        Sb = -C * Rc

        dv = 1.0 / (Sb * Rc - Sc * Rb)
        delb = (R * Sc - S * Rc) * dv
        delc = (-R * Sb + S * Rb) * dv
        diff = np.array([delb, delc])

        B = B + diff[0]
        C = C + diff[1]
        T = np.poly1d([1, B, C])
        V = T.roots
        k = k + 1
    return V, k

def polygcd(a,b, eps=1e-6):
    '''return monic GCD of polynomials a and b'''
    pa = a
    pb = b
    M = lambda x: x/x[0]
    # find gcd of a and b
    while len(pb) > 1 or pb[0] != 0:
        # Danger Will Robinson! requires numerical equality
        q,r = np.polydiv(pa,pb)
        pa = pb
        pb = r
    return M(pa)

def Bairstow(P,epsilon):
    k = 0; # compteur d'iteration;
    A = P.coeffs;
    #B = 1; C = -2;
    B = A[1]/A[0];
    C = A[2]/A[0];
    V = np.zeros(2);
    while ((abs(P(V[1])) > epsilon) & (abs(P(V[0])) > epsilon)):
        U = np.array([B, C]);
        T = np.poly1d([1, B, C]);
        
        Div = np.polydiv(P,T) #1ere div
        Reste = Div[1];
        Q = Div[0];       
        R = Reste[0];
        S = Reste[1];
        
        Div = np.polydiv(Q,T); #2nde div
        Reste = Div[1];
        G = Div[0];
        Rc = -Reste[0];
        Sc = -Reste[1];
        
        Rb = -B*Rc + Sc;
        Sb = -C*Rc;
        
        dv = 1.0/(Sb*Rc -Sc*Rb);
        delb = (R*Sc - S*Rc)*dv;
        delc = (-R*Sb + S*Rb)*dv;
        diff = np.array([delb, delc])
        
        
        B = B + diff[0];
        C = C + diff[1];        
        T = np.poly1d([1, B, C]);
        V = T.roots;
        k = k+ 1;
    return V,k;

def Bairstow2(P):  # Ici la resolution se fait avec Newton_Raphson
    k = 0
    B = 100
    C = -110
    epsilon = 10 ** -10
    V = np.zeros(2)
    while abs(P(V[1])) > epsilon and abs(P(V[0])) > epsilon:
        U = np.array([B, C])
        T = np.poly1d([1, B, C])

        Div = np.polydiv(P, T)  # 1ere div
        Reste = Div[1]
        Q = Div[0]
        R = Reste[0]
        S = Reste[1]

        Div = np.polydiv(Q, T)
        # 2nde div
        Reste = Div[1]
        G = Div[0]
        Rc = -Reste[0]
        Sc = -Reste[1]

        Rb = -B * Rc + Sc
        Sb = -C * Rc

        f = lambda x, y: np.array([R, S])
        J = np.array([[Rb, Rc], [Sb, Sc]])
        diff = Newton_Raphson(f, J, U, 100, 10 ** -3)

        B = B + diff[0]
        C = C + diff[1]
        T = np.poly1d([1, B, C])
        V = T.roots
        k = k + 1
    return V, k

    def test_poly1d_math(self):
        # here we use some simple coeffs to make calculations easier
        p = np.poly1d([1., 2, 4])
        q = np.poly1d([4., 2, 1])
        assert_equal(p/q, (np.poly1d([0.25]), np.poly1d([1.5, 3.75])))
        assert_equal(p.integ(), np.poly1d([1/3, 1., 4., 0.]))
        assert_equal(p.integ(1), np.poly1d([1/3, 1., 4., 0.]))

        p = np.poly1d([1., 2, 3])
        q = np.poly1d([3., 2, 1])
        assert_equal(p * q, np.poly1d([3., 8., 14., 8., 3.]))
        assert_equal(p + q, np.poly1d([4., 4., 4.]))
        assert_equal(p - q, np.poly1d([-2., 0., 2.]))
        assert_equal(p ** 4, np.poly1d([1., 8., 36., 104., 214., 312., 324., 216., 81.]))
        assert_equal(p(q), np.poly1d([9., 12., 16., 8., 6.]))
        assert_equal(q(p), np.poly1d([3., 12., 32., 40., 34.]))
        assert_equal(p.deriv(), np.poly1d([2., 2.]))
        assert_equal(p.deriv(2), np.poly1d([2.]))
        assert_equal(np.polydiv(np.poly1d([1, 0, -1]), np.poly1d([1, 1])),
                     (np.poly1d([1., -1.]), np.poly1d([0.])))

def polyInverse(num,field,modValue = 2,printOption = False):
	a = num
	ansOld = np.poly1d(a)

	b = field
	ans = np.poly1d(b)

	first = np.poly1d(0)
	sec = np.poly1d(1)

	otherSol = np.poly1d(1)
	secotherSol = np.poly1d(0)

	while  ansOld(0) != 0:
		array = np.polydiv(ans,ansOld)[0]
		array = modPolynomial(array)

		inv = compute(sec,first,array)
		first = sec
		sec = inv


		othInv = compute(secotherSol,otherSol,array)
		otherSol = secotherSol
		secotherSol = othInv
		
		newAns = compute(ansOld,ans,array)
		ans = ansOld
		ansOld = newAns

	list = (inv, othInv)

	if printOption:
		print 'the inverse of \n ',np.poly1d(field),' is \n',inv
		print 'and the other inverse is \n',othInv
	
	return list

def GetPolyEvalData(degrees, batch_size=128):
    """ Generate one batch of data
    """
    assert type(degrees) is list
    if FLAGS.task == "root_eval":
        return GetSmallestRootData(degrees, batch_size)
    degree = np.random.choice(degrees)
    X = np.zeros((degree, batch_size, 2), dtype=np.float32)
    if FLAGS.task.endswith("eval"):
        Y = np.zeros((batch_size, 1), dtype=np.float32)
    else:
        Y = np.zeros((degree, batch_size, 1), dtype=np.float32)
    for i in xrange(batch_size):
        roots = np.random.uniform(low=-1, high=1, size=degree - 1)
        roots.sort()
        coefs_0_n = np.polynomial.polynomial.polyfromroots(roots)
        f = np.poly1d(roots, True)
        coefs_0_n = np.random.uniform(low=-1, high=1, size=degree)
        f = np.poly1d(coefs_0_n[::-1])
        a = np.random.uniform(low=-1, high=1)
        if FLAGS.task == "poly_eval":
            Y[i, 0] = f(a)
        elif FLAGS.task == "poly_der_eval":
            Y[i, 0] = np.polyder(f)(a)
        elif FLAGS.task == "poly_div":
            Y[:, i, 0] = np.concatenate(np.polydiv(f, [1, -a]))
        elif FLAGS.task == "newton_eval":
            Y[i, 0] = a - f(a) / np.polyder(f)(a)
        else:
            raise ValueError("Task must be either poly_eval, poly_div, " +
                             "poly_der_eval, root_eval, or newton_eval")
        X[:, i, 0] = coefs_0_n[::-1]
        X[:, i, 1] = a
    if FLAGS.task.endswith("eval"):
        return list(X), Y
    else:
        return list(X), list(Y)

matriz = []

matriz.append(eval(input("Informe os coeficientes do polinômio: ")))

a = eval(input("Entre com o limite inferior: "))
b = eval(input("Entre com o limite superior: "))

size = len(matriz[0])
ordem = size -2

matriz.append(np.polyder(matriz[0]))

i = 0
while ordem:
  matriz.append(np.polydiv(matriz[i], matriz[i+1])[1]*-1)
  i += 1
  ordem -= 1

i = 0
mul = size - 1

resultadoA = []
resultadoB = []

while size:
  resultadoA.append(np.polyval(matriz[i], a))
  resultadoB.append(np.polyval(matriz[i], b))
  i += 1
  size -= 1

def residuez(b,a,tol=1e-3,rtype='avg'):
    """Compute partial-fraction expansion of b(z) / a(z).

    If M = len(b) and N = len(a)

            b(z)     b[0] + b[1] z**(-1) + ... + b[M-1] z**(-M+1)
    H(z) = ------ = ----------------------------------------------
            a(z)     a[0] + a[1] z**(-1) + ... + a[N-1] z**(-N+1)

                 r[0]                   r[-1]
         = --------------- + ... + ---------------- + k[0] + k[1]z**(-1) ...
           (1-p[0]z**(-1))         (1-p[-1]z**(-1))

    If there are any repeated roots (closer than tol), then the partial
    fraction expansion has terms like

               r[i]              r[i+1]                    r[i+n-1]
          -------------- + ------------------ + ... + ------------------
          (1-p[i]z**(-1))  (1-p[i]z**(-1))**2         (1-p[i]z**(-1))**n

    See also:  invresz, poly, polyval, unique_roots
    """
    b,a = map(asarray,(b,a))
    gain = a[0]
    brev, arev = b[::-1],a[::-1]
    krev,brev = polydiv(brev,arev)
    if krev == []:
        k = []
    else:
        k = krev[::-1]
    b = brev[::-1]
    p = roots(a)
    r = p*0.0
    pout, mult = unique_roots(p,tol=tol,rtype=rtype)
    p = []
    for n in range(len(pout)):
        p.extend([pout[n]]*mult[n])
    p = asarray(p)
    # Compute the residue from the general formula (for discrete-time)
    #  the polynomial is in z**(-1) and the multiplication is by terms
    #  like this (1-p[i] z**(-1))**mult[i].  After differentiation,
    #  we must divide by (-p[i])**(m-k) as well as (m-k)!
    indx = 0
    for n in range(len(pout)):
        bn = brev.copy()
        pn = []
        for l in range(len(pout)):
            if l != n:
                pn.extend([pout[l]]*mult[l])
        an = atleast_1d(poly(pn))[::-1]
        # bn(z) / an(z) is (1-po[n] z**(-1))**Nn * b(z) / a(z) where Nn is
        # multiplicity of pole at po[n] and b(z) and a(z) are polynomials.
        sig = mult[n]
        for m in range(sig,0,-1):
            if sig > m:
                # compute next derivative of bn(s) / an(s)
                term1 = polymul(polyder(bn,1),an)
                term2 = polymul(bn,polyder(an,1))
                bn = polysub(term1,term2)
                an = polymul(an,an)
            r[indx+m-1] = polyval(bn,1.0/pout[n]) / polyval(an,1.0/pout[n]) \
                          / factorial(sig-m) / (-pout[n])**(sig-m)
        indx += sig
    return r/gain, p, k

 def test_poly_div(self, level=rlevel):
     """Ticket #553"""
     u = np.poly1d([1, 2, 3])
     v = np.poly1d([1, 2, 3, 4, 5])
     q, r = np.polydiv(u, v)
     assert_equal(q*v + r, u)

def residue(b, a, tol=1e-3, rtype='avg'):
    """
    Compute partial-fraction expansion of b(s) / a(s).

    If ``M = len(b)`` and ``N = len(a)``, then the partial-fraction
    expansion H(s) is defined as::

              b(s)     b[0] s**(M-1) + b[1] s**(M-2) + ... + b[M-1]
      H(s) = ------ = ----------------------------------------------
              a(s)     a[0] s**(N-1) + a[1] s**(N-2) + ... + a[N-1]

               r[0]       r[1]             r[-1]
           = -------- + -------- + ... + --------- + k(s)
             (s-p[0])   (s-p[1])         (s-p[-1])

    If there are any repeated roots (closer together than `tol`), then H(s)
    has terms like::

            r[i]      r[i+1]              r[i+n-1]
          -------- + ----------- + ... + -----------
          (s-p[i])  (s-p[i])**2          (s-p[i])**n

    Returns
    -------
    r : ndarray
        Residues.
    p : ndarray
        Poles.
    k : ndarray
        Coefficients of the direct polynomial term.

    See Also
    --------
    invres, numpy.poly, unique_roots

    """

    b, a = map(asarray, (b, a))
    rscale = a[0]
    k, b = polydiv(b, a)
    p = roots(a)
    r = p * 0.0
    pout, mult = unique_roots(p, tol=tol, rtype=rtype)
    p = []
    for n in range(len(pout)):
        p.extend([pout[n]] * mult[n])
    p = asarray(p)
    # Compute the residue from the general formula
    indx = 0
    for n in range(len(pout)):
        bn = b.copy()
        pn = []
        for l in range(len(pout)):
            if l != n:
                pn.extend([pout[l]] * mult[l])
        an = atleast_1d(poly(pn))
        # bn(s) / an(s) is (s-po[n])**Nn * b(s) / a(s) where Nn is
        # multiplicity of pole at po[n]
        sig = mult[n]
        for m in range(sig, 0, -1):
            if sig > m:
                # compute next derivative of bn(s) / an(s)
                term1 = polymul(polyder(bn, 1), an)
                term2 = polymul(bn, polyder(an, 1))
                bn = polysub(term1, term2)
                an = polymul(an, an)
            r[indx + m - 1] = polyval(bn, pout[n]) / polyval(an, pout[n]) \
                          / factorial(sig - m)
        indx += sig
    return r / rscale, p, k

def poly2lsf(a):
    """Prediction polynomial to line spectral frequencies.

    converts the prediction polynomial specified by A,
    into the corresponding line spectral frequencies, LSF. 
    normalizes the prediction polynomial by A(1).

    .. doctest::

        >>> a = [1.0000  ,  0.6149   , 0.9899   , 0.0000 ,   0.0031,   -0.0082
        >>> lsf = poly2lsf(a)
        >>> lsf =  array([  0.7842,    1.5605 ,   1.8776 ,   1.8984,    2.3593])

    .. seealso:: lsf2poly, poly2rc, poly2qc, rc2is
    """

    #Line spectral frequencies are not defined for complex polynomials.

    # Normalize the polynomial

    a = numpy.array(a)
    if a[0] != 1:
        a/=a[0]

    if max(numpy.abs(numpy.roots(a))) >= 1.0:
        print("hello world")

    # Form the sum and differnce filters

    p  = len(a)-1   # The leading one in the polynomial is not used
    a1 = numpy.concatenate((a, numpy.array([0])))        
    a2 = a1[-1::-1]
    P1 = a1 - a2        # Difference filter
    Q1 = a1 + a2        # Sum Filter 
    divisorOdd = [1,0,-1]
    divisorEven1 = [1,-1]
    divisorEven2 = [1,1]

    # If order is even, remove the known root at z = 1 for P1 and z = -1 for Q1
    # If odd, remove both the roots from P1

    if p%2: # Odd order
        P, r = scipy.signal.deconvolve(P1,divisorOdd)
        Q = Q1
    else:          # Even order 
        P, r = numpy.polydiv(P1, divisorEven1)
        Q, r = numpy.polydiv(Q1, divisorEven2)
    
    rP  = numpy.roots(P)
    rQ  = numpy.roots(Q)

    aP  = numpy.angle(rP[1::2])
    aQ  = numpy.angle(rQ[1::2])

    lsf = sorted(numpy.concatenate((-aP,-aQ)))
    print("a",a)
    print("p",p)
    print("a1",a1)
    print("a2",a2)
    print("P1",P1)
    print("Q1",Q1)
    print("r",r)
    print("P",P)
    print("Q",Q)
    print("rP",rP)
    print("rQ",rQ)
    print("aP",aP)
    print("aQ",aQ)

    return lsf