Python Poly.all_coeffs方法代码示例

# 需要导入模块: from sympy import Poly [as 别名]
# 或者: from sympy.Poly import all_coeffs [as 别名]
    def _eval_expand_func(self, **hints):
        from sympy import exp, I, floor, Add, Poly, Dummy, exp_polar, unpolarify
        z, s, a = self.args
        if z == 1:
            return zeta(s, a)
        if s.is_Integer and s <= 0:
            t = Dummy('t')
            p = Poly((t + a)**(-s), t)
            start = 1/(1 - t)
            res = S(0)
            for c in reversed(p.all_coeffs()):
                res += c*start
                start = t*start.diff(t)
            return res.subs(t, z)

        if a.is_Rational:
            # See section 18 of
            #   Kelly B. Roach.  Hypergeometric Function Representations.
            #   In: Proceedings of the 1997 International Symposium on Symbolic and
            #   Algebraic Computation, pages 205-211, New York, 1997. ACM.
            # TODO should something be polarified here?
            add = S(0)
            mul = S(1)
            # First reduce a to the interaval (0, 1]
            if a > 1:
                n = floor(a)
                if n == a:
                    n -= 1
                a -= n
                mul = z**(-n)
                add = Add(*[-z**(k - n)/(a + k)**s for k in xrange(n)])
            elif a <= 0:
                n = floor(-a) + 1
                a += n
                mul = z**n
                add = Add(*[z**(n - 1 - k)/(a - k - 1)**s for k in xrange(n)])

            m, n = S([a.p, a.q])
            zet = exp_polar(2*pi*I/n)
            root = z**(1/n)
            return add + mul*n**(s - 1)*Add(
                *[polylog(s, zet**k*root)._eval_expand_func(**hints)
                  / (unpolarify(zet)**k*root)**m for k in xrange(n)])

        # TODO use minpoly instead of ad-hoc methods when issue 2789 is fixed
        if z.func is exp and (z.args[0]/(pi*I)).is_Rational or z in [-1, I, -I]:
            # TODO reference?
            if z == -1:
                p, q = S([1, 2])
            elif z == I:
                p, q = S([1, 4])
            elif z == -I:
                p, q = S([-1, 4])
            else:
                arg = z.args[0]/(2*pi*I)
                p, q = S([arg.p, arg.q])
            return Add(*[exp(2*pi*I*k*p/q)/q**s*zeta(s, (k + a)/q)
                         for k in xrange(q)])

        return lerchphi(z, s, a)

# 需要导入模块: from sympy import Poly [as 别名]
# 或者: from sympy.Poly import all_coeffs [as 别名]
def gauss_lobatto_points(p):
    """
    Returns the list of Gauss-Lobatto points of the order 'p'.
    """
    x = Symbol("x")
    print "creating"
    e = legendre(p, x).diff(x)
    e = Poly(e, x)
    print "polydone"
    if e == 0:
        return []
    print "roots"
    if e == 1:
        r = []
    else:
        with workdps(40):
            r, err = polyroots(e.all_coeffs(), error=True)
        #if err > 1e-40:
        #    raise Exception("Internal Error: Root is not precise")
    print "done"
    p = []
    p.append("-1.0")
    for x in r:
        if abs(x) < 1e-40: x = 0
        p.append(str(x))
    p.append("1.0")
    return p

# 需要导入模块: from sympy import Poly [as 别名]
# 或者: from sympy.Poly import all_coeffs [as 别名]
 def linear_arg(arg):
     """ Test if arg is of form a*s+b, raise exception if not. """
     if not arg.is_polynomial(s):
         raise exception(fact)
     p = Poly(arg, s)
     if p.degree() != 1:
         raise exception(fact)
     return p.all_coeffs()

# 需要导入模块: from sympy import Poly [as 别名]
# 或者: from sympy.Poly import all_coeffs [as 别名]
def check_convergence(numer, denom, n):
    """
    Returns (h, g, p) where
    -- h is:
        > 0 for convergence of rate 1/factorial(n)**h
        < 0 for divergence of rate factorial(n)**(-h)
        = 0 for geometric or polynomial convergence or divergence

    -- abs(g) is:
        > 1 for geometric convergence of rate 1/h**n
        < 1 for geometric divergence of rate h**n
        = 1 for polynomial convergence or divergence

        (g < 0 indicates an alternating series)

    -- p is:
        > 1 for polynomial convergence of rate 1/n**h
        <= 1 for polynomial divergence of rate n**(-h)

    """
    from sympy import Poly
    npol = Poly(numer, n)
    dpol = Poly(denom, n)
    p = npol.degree()
    q = dpol.degree()
    rate = q - p
    if rate:
        return rate, None, None
    constant = dpol.LC() / npol.LC()
    if abs(constant) != 1:
        return rate, constant, None
    if npol.degree() == dpol.degree() == 0:
        return rate, constant, 0
    pc = npol.all_coeffs()[1]
    qc = dpol.all_coeffs()[1]
    return rate, constant, (qc - pc)/dpol.LC()

# 需要导入模块: from sympy import Poly [as 别名]
# 或者: from sympy.Poly import all_coeffs [as 别名]
def gauss_lobatto_points(p):
    """
    Returns the list of Gauss-Lobatto points of the order 'p'.
    """
    x = Symbol("x")
    e = (1-x**2)*legendre(p, x).diff(x)
    e = Poly(e, x)
    if e == 0:
        return []
    with workdps(40):
        r, err = polyroots(e.all_coeffs(), error=True)
    if err > 1e-40:
        raise Exception("Internal Error: Root is not precise")
    p = []
    for x in r:
        if abs(x) < 1e-40: x = 0
        p.append(str(x))
    return p

# 需要导入模块: from sympy import Poly [as 别名]
# 或者: from sympy.Poly import all_coeffs [as 别名]

#.........这里部分代码省略.........
        b_ *= s_multiplier

    # 2)
    numer, denom = f.as_numer_denom()
    numer = Mul.make_args(numer)
    denom = Mul.make_args(denom)
    args = zip(numer, repeat(True)) + zip(denom, repeat(False))

    facs = []
    dfacs = []
    # *_gammas will contain pairs (a, c) representing Gamma(a*s + c)
    numer_gammas = []
    denom_gammas = []
    # exponentials will contain bases for exponentials of s
    exponentials = []
    def exception(fact):
        return IntegralTransformError("Inverse Mellin", f, "Unrecognised form '%s'." % fact)
    while args:
        fact, is_numer = args.pop()
        if is_numer:
            ugammas, lgammas = numer_gammas, denom_gammas
            ufacs, lfacs = facs, dfacs
        else:
            ugammas, lgammas = denom_gammas, numer_gammas
            ufacs, lfacs = dfacs, facs

        def linear_arg(arg):
            """ Test if arg is of form a*s+b, raise exception if not. """
            if not arg.is_polynomial(s):
                raise exception(fact)
            p = Poly(arg, s)
            if p.degree() != 1:
                raise exception(fact)
            return p.all_coeffs()

        # constants
        if not fact.has(s):
            ufacs += [fact]
        # exponentials
        elif fact.is_Pow or isinstance(fact, exp_):
            if fact.is_Pow:
                base = fact.base
                exp  = fact.exp
            else:
                base = exp_polar(1)
                exp  = fact.args[0]
            if exp.is_Integer:
                cond = is_numer
                if exp < 0:
                    cond = not cond
                args += [(base, cond)]*abs(exp)
                continue
            elif not base.has(s):
                a, b = linear_arg(exp)
                if not is_numer:
                    base = 1/base
                exponentials += [base**a]
                facs += [base**b]
            else:
                raise exception(fact)
        # linear factors
        elif fact.is_polynomial(s):
            p = Poly(fact, s)
            if p.degree() != 1:
                # We completely factor the poly. For this we need the roots.
                # Now roots() only works in some cases (low degree), and RootOf

# 需要导入模块: from sympy import Poly [as 别名]
# 或者: from sympy.Poly import all_coeffs [as 别名]
print 'W0: ', W0
print 'B: ', B, '\n'
# Defining a list [0,0,...1] of size N+1 required to generate the coefficients of Nth order Cheyshev Polynomial
c = zeros(N+1)
c[-1] = 1.0
#
TnS = polynomial.chebyshev.Chebyshev(c)
TnS_coeff = polynomial.chebyshev.cheb2poly(TnS.coef)
# print 'Coefficients: ', TnS_coeff
CnW = 0.0
for i in xrange(int(N+1)):
	CnW = CnW + TnS_coeff[i]*(-1j*s)**i
# HcS_sq = 1.0/(1.0 + (eps**2)*(CnW**2))
denom_all = Poly(1.0 + (eps**2)*(CnW**2), s)
# print denom_all
denom_all = denom_all.all_coeffs()
leading_coeff = denom_all[0]
denom_all = [i/leading_coeff for i in denom_all]
norm_fac = (leading_coeff/abs(leading_coeff))*pow(abs(leading_coeff),0.5)
poles_all = roots(denom_all)
poles_LHP = zeros(N,dtype=complex64)
j = 0
denominator_LHP = 1.0 + 0j
for i in xrange(int(2*N)):
	if ((poles_all[i]).real<0):
		poles_LHP[j] = poles_all[i]
		j = j+1
		denominator_LHP = expand(denominator_LHP*(s-poles_all[i]))
# print poles
denom_poly = Poly(denominator_LHP,s)
denom_poly = denom_poly.all_coeffs()

# 需要导入模块: from sympy import Poly [as 别名]
# 或者: from sympy.Poly import all_coeffs [as 别名]
# plot(abs(freq_res_LPF), (s,0,100))
print 'HcS: ', HcS
HcS_Bandpass = HcS.subs(s, transform)
# freq_res_BPF = HcS_Bandpass.subs(s,1j*2.0*pi*s)
# plot(abs(freq_res_BPF), (s,0,50))
print 'HcS_Bandpass: ', simplify(HcS_Bandpass)
# bilinear = 2.0*f_sampling*(1.0 - x)/(1.0 + x)
bilinear = (1.0 - x)/(1.0 + x)
Hz = HcS_Bandpass.subs(s, bilinear)
Hz = simplify(Hz)
Hz = cancel(Hz)
print 'Hz: ', Hz
# print 'Numerator: ', fraction(Hz)[0]
# print 'Denominator: ', fraction(Hz)[1]
ax = Poly(fraction(Hz)[0], x) # Separating out the numerator
ax = ax.all_coeffs() # Array of size [(order of ax) + 1] with ax[0] being coeff of x^n
ax = list(reversed(ax)) # Reversing it so that ax[0] is coeff of order x^0
ay = Poly(fraction(Hz)[1], x) # Separating out the denominator
ay = ay.all_coeffs()
ay = list(reversed(ay))
norm_fac = ay[0]
# print 'ay[0]: ', ay[0]
ax = [i/norm_fac for i in ax]
ay = [i/norm_fac for i in ay]
print 'ax: ', ax
print 'ay: ', ay
# print 'Length of ax: ', len(ax)
# print 'Length of ay: ', len(ay)
#
# ax[i] and ay[i] are the coefficients of x[n-i] and y[n-i]  respectively in the difference equation.
# The difference equation is implemented by the function 'filter' on any specified x[n]