Python numpy.polydiv函数代码示例

 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)