def tustin(tf, sRate):
	s =	getMapping(tf)['s']
	T = 1/sRate
	poles = sympy.roots(sympy.denom(tf), s, multiple=True)
	centroid = np.mean(np.abs(poles))
	invz = sympy.Symbol('invz', real=True)
	# normalized center frequency derived from poles of filter SISO transfer function
	w0 = 2*np.pi*centroid/(sRate*2*np.pi)
	# modified bilinear transform w/ frequency warping
	bt = w0/sympy.tan(w0*T/2) * ((1-invz)/(1+invz))
	dt = sympy.simplify(tf.subs({s: bt}))
	b = sympy.Poly(sympy.numer(dt)).all_coeffs()[::-1]
	a = sympy.Poly(sympy.denom(dt)).all_coeffs()[::-1]
	normalize = lambda x: float(x/a[0])
	return (map(normalize, b), map(normalize, a))

def poles(G):
    '''
    Return the poles of a multivariable transfer function system. Applies
    Theorem 4.4 (p135).

    Parameters
    ----------
    G : numpy matrix (n x n)
        The transfer function G(s) of the system.

    Returns
    -------
    zero : array
        List of zeros.

    Example
    -------
    >>> def G(s):
    ...     return 1 / (s + 2) * numpy.matrix([[s - 1,  4],
    ...                                       [4.5, 2 * (s - 1)]])
    >>> poles(G)
    [-2.00000000000000]

    Note
    ----
    Not applicable for a non-squared plant, yet.
    '''

    s = sympy.Symbol('s')
    G = sympy.Matrix(G(s))  # convert to sympy matrix object
    det = sympy.simplify(G.det())
    pole = sympy.solve(sympy.denom(det))
    return pole

def find_simple_recurrence_vector(l):
    """
    This function is used internally by other functions from the
    sympy.concrete.guess module. While most users may want to rather use the
    function find_simple_recurrence when looking for recurrence relations
    among rational numbers, the current function may still be useful when
    some post-processing has to be done.

    The function returns a vector of length n when a recurrence relation of
    order n is detected in the sequence of rational numbers v.

    If the returned vector has a length 1, then the returned value is always
    the list [0], which means that no relation has been found.

    While the functions is intended to be used with rational numbers, it should
    work for other kinds of real numbers except for some cases involving
    quadratic numbers; for that reason it should be used with some caution when
    the argument is not a list of rational numbers.

    Examples
    ========

    >>> from sympy.concrete.guess import find_simple_recurrence_vector
    >>> from sympy import fibonacci
    >>> find_simple_recurrence_vector([fibonacci(k) for k in range(12)])
    [1, -1, -1]

    See also
    ========

    See the function sympy.concrete.guess.find_simple_recurrence which is more
    user-friendly.

    """
    q1 = [0]
    q2 = [Integer(1)]
    b, z = 0, len(l) >> 1
    while len(q2) <= z:
        while l[b]==0:
            b += 1
            if b == len(l):
                c = 1
                for x in q2:
                    c = lcm(c, denom(x))
                if q2[0]*c < 0: c = -c
                for k in range(len(q2)):
                    q2[k] = int(q2[k]*c)
                return q2
        a = Integer(1)/l[b]
        m = [a]
        for k in range(b+1, len(l)):
            m.append(-sum(l[j+1]*m[b-j-1] for j in range(b, k))*a)
        l, m = m, [0] * max(len(q2), b+len(q1))
        for k in range(len(q2)):
            m[k] = a*q2[k]
        for k in range(b, b+len(q1)):
            m[k] += q1[k-b]
        while m[-1]==0: m.pop() # because trailing zeros can occur
        q1, q2, b = q2, m, 1
    return [0]

def e2nd(expression):
	""" basic helper function that accepts a sympy expression, expands it,
	attempts to simplify it, and returns a numerator and denomenator pair for the instantiation of a scipy
	LTI system object. """
	expression = expression.expand()
	expression = expression.cancel()
	n = sympy.Poly(sympy.numer(expression), s).all_coeffs()
	d = sympy.Poly(sympy.denom(expression), s).all_coeffs()
	n = map(float, n)
	d = map(float, d)
	return (n, d)

    def __new__(cls, b, e, evaluate=True):
        from sympy.functions.elementary.exponential import exp_polar
        from sympy.functions import log
        # don't optimize "if e==0; return 1" here; it's better to handle that
        # in the calling routine so this doesn't get called
        b = _sympify(b)
        e = _sympify(e)
        if evaluate:
            if e is S.Zero:
                return S.One
            elif e is S.One:
                return b
            elif S.NaN in (b, e):
                if b is S.One:  # already handled e == 0 above
                    return S.One
                return S.NaN
            else:
                if e.func == log:
                    if len(e.args) == 2:
                        lbase = e.args[1]
                    else:
                        lbase = S.Exp1
                    if lbase == b:
                        return e.args[0]

                if e is Mul and e.args[1].func == log:
                    if len(e.args[1].args) == 2:
                        lbase = e.args[1].args[1]
                    else:
                        lbase = S.Exp1
                    if lbase == b:
                        return e.args[1].args[0] ** e.args[0]

                # recognize base as E
                if not e.is_Atom and b is not S.Exp1 and b.func is not exp_polar:
                    from sympy import numer, denom, log, sign, im, factor_terms
                    c, ex = factor_terms(e, sign=False).as_coeff_Mul()
                    den = denom(ex)
                    if den.func is log and den.args[0] == b:
                        return S.Exp1**(c*numer(ex))
                    elif den.is_Add:
                        s = sign(im(b))
                        if s.is_Number and s and den == \
                                log(-factor_terms(b, sign=False)) + s*S.ImaginaryUnit*S.Pi:
                            return S.Exp1**(c*numer(ex))

                obj = b._eval_power(e)
                if obj is not None:
                    return obj
        obj = Expr.__new__(cls, b, e)
        obj.is_commutative = (b.is_commutative and e.is_commutative)
        return obj

    def Pow(expr, assumptions):
        """
        Imaginary**integer/odd  -> Imaginary
        Imaginary**integer/even -> Real if integer % 2 == 0
        b**Imaginary            -> !Imaginary if exponent is an integer multiple of I*pi/log(b)
        Imaginary**Real         -> ?
        Negative**even root     -> Imaginary
        Negative**odd root      -> Real
        Negative**Real          -> Imaginary
        Real**Integer           -> Real
        Real**Positive          -> Real
        """
        if expr.is_number:
            return AskImaginaryHandler._number(expr, assumptions)

        if expr.base.func == C.exp:
            if ask(Q.imaginary(expr.base.args[0]), assumptions):
                if ask(Q.imaginary(expr.exp), assumptions):
                    return False
                i = expr.base.args[0]/I/pi
                if ask(Q.integer(2*i), assumptions):
                    return ask(Q.imaginary(((-1)**i)**expr.exp), assumptions)

        if ask(Q.imaginary(expr.base), assumptions):
            if ask(Q.integer(expr.exp), assumptions):
                odd = ask(Q.odd(expr.exp), assumptions)
                if odd is not None:
                    return odd
                return

        if ask(Q.imaginary(expr.exp), assumptions):
            imlog = ask(Q.imaginary(C.log(expr.base)), assumptions)
            if imlog is not None:
                return False  # I**i -> real; (2*I)**i -> complex ==> not imaginary

        if ask(Q.real(expr.base), assumptions):
            if ask(Q.real(expr.exp), assumptions):
                if ask(Q.rational(expr.exp) & Q.even(denom(expr.exp)), assumptions):
                    return ask(Q.negative(expr.base), assumptions)
                elif ask(Q.integer(expr.exp), assumptions):
                    return False
                elif ask(Q.positive(expr.base), assumptions):
                    return False
                elif ask(Q.negative(expr.base), assumptions):
                    return True

def mutate_rat(value, *, keep_sign=True, mutate_up=None):
    '''Mutate rational number'''

    r = random.random()
    numer, denom = sp.numer(value), sp.denom(value)

    # swap numerator and denominator
    if r < 0.20:
        numer, denom = (mutate_int(denom, mutate_up=mutate_up, keep_sign=keep_sign),
                        mutate_int(numer, mutate_up=mutate_up))
    # mutate numerator
    elif r < 0.65:
        numer = mutate_int(numer, mutate_up=mutate_up, keep_sign=keep_sign)
    # mutate denominator
    else:
        denom = mutate_int(denom, mutate_up=mutate_up, keep_sign=keep_sign)

    return sp.Rational(numer, denom or 1)

    def __new__(cls, b, e, evaluate=None):
        if evaluate is None:
            evaluate = global_evaluate[0]
        from sympy.functions.elementary.exponential import exp_polar

        b = _sympify(b)
        e = _sympify(e)
        if evaluate:
            if e is S.Zero:
                return S.One
            elif e is S.One:
                return b
            elif e.is_integer and _coeff_isneg(b):
                if e.is_even:
                    b = -b
                elif e.is_odd:
                    return -Pow(-b, e)
            if b is S.One:
                if e in (S.NaN, S.Infinity, -S.Infinity):
                    return S.NaN
                return S.One
            elif S.NaN in (b, e):  # XXX S.NaN**x -> S.NaN under assumption that x != 0
                return S.NaN
            else:
                # recognize base as E
                if not e.is_Atom and b is not S.Exp1 and b.func is not exp_polar:
                    from sympy import numer, denom, log, sign, im, factor_terms
                    c, ex = factor_terms(e, sign=False).as_coeff_Mul()
                    den = denom(ex)
                    if den.func is log and den.args[0] == b:
                        return S.Exp1**(c*numer(ex))
                    elif den.is_Add:
                        s = sign(im(b))
                        if s.is_Number and s and den == \
                                log(-factor_terms(b, sign=False)) + s*S.ImaginaryUnit*S.Pi:
                            return S.Exp1**(c*numer(ex))

                obj = b._eval_power(e)
                if obj is not None:
                    return obj
        obj = Expr.__new__(cls, b, e)
        obj.is_commutative = (b.is_commutative and e.is_commutative)
        return obj

def symmetryDetection(allVariables, diffEquations, observables, obsFunctions, initFunctions,
						predictions, predFunctions, ansatz = 'uni', pMax = 2, inputs = [], 
						fixed = [], parallel = 1, allTrafos = False, timeTrans = False,
						pretty = True, suffix=''):
	
	n = len(allVariables)
	m = len(diffEquations)
	h = len(observables)

	###########################################################################################
	#############################     prepare equations    ####################################
	###########################################################################################
	sys.stdout.write('Preparing equations...')
	sys.stdout.flush()

	# make infinitesimal ansatz
	infis, diffInfis, rs = makeAnsatz(ansatz, allVariables, m, len(inputs), pMax, fixed)

	# get infinitesimals of time transformation
	if timeTrans:
		rs.append(spy.var('r_T_1'))
		diffInfiT = rs[-1]
		allVariables += [T]
	else:
		diffInfiT = None

	# and convert to polynomial
	infis, diffInfis = transformInfisToPoly(infis, diffInfis, allVariables, rs, parallel, ansatz)
	
	diffInfiT = Apoly(diffInfiT, allVariables, rs)

	### extract numerator and denominator of equations
	#differential equations
	numerators = [0]*m
	denominators = [0]*m
	for k in range(m):
		rational = spy.together(diffEquations[k])
		numerators[k] = Apoly(spy.numer(rational), allVariables, None)
		denominators[k] = Apoly(spy.denom(rational), allVariables, None)

	#observation functions
	obsNumerators = [0]*h
	obsDenominatros = [0]*h
	for k in range(h):
		rational = spy.together(obsFunctions[k])
		obsNumerators[k] = Apoly(spy.numer(rational), allVariables, None)
		obsDenominatros[k] = Apoly(spy.denom(rational), allVariables, None)

	#initial functions
	if len(initFunctions) != 0:
		initNumerators = [0]*m
		initDenominatros = [0]*m
		for k in range(m):
			rational = spy.together(initFunctions[k])
			initNumerators[k] = Apoly(spy.numer(rational), allVariables, None)
			initDenominatros[k] = Apoly(spy.denom(rational), allVariables, None)
	else:
		initNumerators = []
		initDenominatros = []

	### calculate numerator of derivatives of equations
	#differential equatioins
	derivativesNum = [0]*m
	for i in range(m):
		derivativesNum[i] = [0]*n
	for k in range(m):
		for l in range(n):
			derivativesNum[k][l] = Apoly(None, allVariables, None)
			derivativesNum[k][l].add(numerators[k].diff(l).mul(denominators[k]))
			derivativesNum[k][l].sub(numerators[k].mul(denominators[k].diff(l)))

	#observation functions
	obsDerivativesNum = [0]*h
	for i in range(h):
		obsDerivativesNum[i] = [0]*n
	for k in range(h):
		for l in range(n):
			obsDerivativesNum[k][l] = Apoly(None, allVariables, None)
			obsDerivativesNum[k][l].add(obsNumerators[k].diff(l).mul(obsDenominatros[k]))
			obsDerivativesNum[k][l].sub(obsNumerators[k].mul(obsDenominatros[k].diff(l)))

	#initial functions
	if len(initFunctions) != 0:
		initDerivativesNum = [0]*len(initFunctions)
		for i in range(m):
			initDerivativesNum[i] = [0]*n
		for k in range(m):
			for l in range(n):
				initDerivativesNum[k][l] = Apoly(None, allVariables, None)
				initDerivativesNum[k][l].add(initNumerators[k].diff(l).mul(initDenominatros[k]))
				initDerivativesNum[k][l].sub(initNumerators[k].mul(initDenominatros[k].diff(l)))
	else:
		initDerivativesNum = []

	sys.stdout.write('\rPreparing equations...done\n')
	sys.stdout.flush()

	###########################################################################################
	############################     build linear system    ###################################
	###########################################################################################
#.........这里部分代码省略.........

m2_2 = G[:,[0,2]]
m2 = m2_1.row_join(m2_2)
m2_3 = G[:,[0,1]]
m2 = m2.row_join(m2_3)
print(m2_1)
print(m2_2)
print(m2_3)

# combine order one and two minors

print('---All poles---')
m1count = np.shape(m1)[0]
#find out how many roots there are in total between all 1st order minors
n1Roots = 0
for m in range(m1count):
    n1Roots = n1Roots + sp.degree(sp.denom(m1[m]),s)

#find out how many roots there are in total between all 2nd order minors
n2Roots = 0
r,c = np.shape(m2)
for i in range(r):
    for j in range(c):
        n2Roots = n2Roots + sp.degree(sp.denom(m2[i,j]),s)
print(n2Roots)

# calculate and find common roots 
roots_denom = [None]*(n1Roots + n2Roots)
n = 0
for m in range(m1count):
    denom_roots = sp.solve(sp.denom(m1[m]))
    for i  in range(len(denom_roots)):

def test_issue_5933():
    from sympy import Polygon, RegularPolygon, denom
    x = Polygon(*RegularPolygon((0, 0), 1, 5).vertices).centroid.x
    assert abs(denom(x).n()) > 1e-12
    assert abs(denom(radsimp(x))) > 1e-12  # in case simplify didn't handle it

def gen_lobatto(max_order):
    assert max_order > 2

    x = sm.symbols('x')

    lobs = [0, 1]
    lobs[0] = (1 - x) / 2
    lobs[1] = (1 + x) / 2

    dlobs = [lob.diff('x') for lob in lobs]

    legs = [sm.legendre(0, 'y')]
    clegs = [sm.ccode(legs[0])]
    dlegs = [sm.legendre(0, 'y').diff('y')]
    cdlegs = [sm.ccode(dlegs[0])]

    clobs = [sm.ccode(lob) for lob in lobs]
    cdlobs = [sm.ccode(dlob) for dlob in dlobs]

    denoms = [] # for lobs.

    for ii in range(2, max_order + 1):
        coef = sm.sympify('sqrt(2 * (2 * %s - 1)) / 2' % ii)
        leg = sm.legendre(ii - 1, 'y')

        pleg = leg.as_poly()
        coefs = pleg.all_coeffs()
        denom = max(sm.denom(val) for val in coefs)

        cleg = sm.ccode(sm.horner(leg*denom)/denom)

        dleg = leg.diff('y')
        cdleg = sm.ccode(sm.horner(dleg*denom)/denom)

        lob = sm.simplify(coef * sm.integrate(leg, ('y', -1, x)))
        lobnc = sm.simplify(sm.integrate(leg, ('y', -1, x)))

        plobnc = lobnc.as_poly()
        coefs = plobnc.all_coeffs()
        denom = sm.denom(coef) * max(sm.denom(val) for val in coefs)

        clob = sm.ccode(sm.horner(lob*denom)/denom)

        dlob = lob.diff('x')
        cdlob = sm.ccode(sm.horner(dlob*denom)/denom)

        legs.append(leg)
        clegs.append(cleg)
        dlegs.append(dleg)
        cdlegs.append(cdleg)
        lobs.append(lob)
        clobs.append(clob)
        dlobs.append(dlob)
        cdlobs.append(cdlob)
        denoms.append(denom)

    coef = sm.sympify('sqrt(2 * (2 * %s - 1)) / 2' % (max_order + 1))
    leg = sm.legendre(max_order, 'y')

    pleg = leg.as_poly()
    coefs = pleg.all_coeffs()
    denom = max(sm.denom(val) for val in coefs)

    cleg = sm.ccode(sm.horner(leg*denom)/denom)

    dleg = leg.diff('y')
    cdleg = sm.ccode(sm.horner(dleg*denom)/denom)

    legs.append(leg)
    clegs.append(cleg)
    dlegs.append(dleg)
    cdlegs.append(cdleg)

    kerns = []
    ckerns = []
    dkerns = []
    cdkerns = []
    for ii, lob in enumerate(lobs[2:]):
        kern = sm.simplify(lob / (lobs[0] * lobs[1]))
        dkern = kern.diff('x')

        denom = denoms[ii] / 4
        ckern = sm.ccode(sm.horner(kern*denom)/denom)
        cdkern = sm.ccode(sm.horner(dkern*denom)/denom)

        kerns.append(kern)
        ckerns.append(ckern)
        dkerns.append(dkern)
        cdkerns.append(cdkern)

    return (legs, clegs, dlegs, cdlegs,
            lobs, clobs, dlobs, cdlobs,
            kerns, ckerns, dkerns, cdkerns,
            denoms)

	def simplify(self):
		self._tuple = tuple(sympy.factor(ex) for ex in self._tuple)
		m = functools.reduce(lambda x, y: x * sympy.denom(y), self._tuple, 1)
		self._tuple = tuple(ex * m for ex in self._tuple)
		d = functools.reduce(sympy.gcd, self._tuple)
		self._tuple = tuple(sympy.factor(ex / d) for ex in self._tuple)

tf = ahkabHelpers.reduceTF(tf, mycircuit)

sRate = 44.1e3

fs = lambda x: sympy.N(abs(tf.subs({s:sympy.I*x})))
ws = np.logspace(1, np.log10(sRate), 5e2)
mags = map(fs,ws)

print tf
(b, a) = ahkabHelpers.tustin(tf, sRate)
print b,a

pylab.semilogx(ws, map(fs, ws), 'v', label="from transfer function")
pylab.semilogx(r['ac']['w'][::10], np.abs(r['ac']['VU1o'][::10]), '-', label='from AC simulation')
pylab.vlines(np.abs(sympy.roots(sympy.denom(tf), s, multiple=True)), 0, 1, 'r')

# build white noise input vector, normalized to \pm 1.0
x = list(2*np.random.random(sRate)-1.0)
# allocate output vector for y[n-2] indexing to work
y = [0.0]*len(x)

# Direct Form I
for n in range(len(x)):
	y[n] = b[0]*x[n] + b[1]*x[n-1] + b[2]*x[n-2] - a[1]*y[n-1] - a[2]*y[n-2]

# frequencies to radians
freqs = (np.fft.fftfreq(len(y), 1/sRate))[0:len(y)/2] * (2*np.pi)
mags = 	(np.abs(np.fft.fft(np.array(y))))[0:len(y)/2]

# normalize magnitude to 1.0

def ac_analysis(param_d, param_l, instance, file_sufix):
    """ Performs ac analysis

        param_d: substitutions for symbols, or named parameters for plot
        
        param_l: expresions to plot

        format is:
        .ac expresion0 [expresion1 expresion2 ...] sweep = parameter_to_sweep [symmbol_or_option0 = value0
        symmbol_or_option1 = value1 ...]
        
        expresions are on positiona parameters list can hold expresions containing parameters or functions of nodal
        voltages [v(node)] and element and port currents [i(element) isub(port)]

        named parameters (param_d) can contain options for analysis or substitutions for symbols, substituions will be
        done as they are without any parsing, be aware of symbol and option names clashes.

        Config options:
        fstart:         first value of frequency for AC analysis [float]
        fstop:          last value of frequency for AC analysis [float]
        fscale:         scale for frequency points [linear | log]
        npoints:        numbers of points for ac analysis [integer]
        yscale:         scale for y-axis to being displayed on [linear or log]
        hold:           hold plot for next analysis and don't save it to file [yes | no]
        show_poles:     show poles of function on plot [yes | no]
        show_zeroes:    show zeros of function on plot [yes | no]
        title:          display title above ac plot [string]
        show_legend:    show legend on plot [yes | no]
        xkcd:           style plot to be xkcd like scetch
    """
    warnings.filterwarnings('ignore')  # Just getting rid of those fake casting from complex warnings
    s, w = sympy.symbols(('s', 'w'))
    config = {'fstart': 1,
              'fstop': 1e6,
              'fscale': 'log',
              'npoints': 100,
              'yscale': 'log',
              'type': 'amp',
              'hold': 'no',
              'show_poles': 'yes',
              'show_zeros': 'yes',
              'title': None,
              'show_legend': 'no',
              'xkcd': 'no'}

    for config_name in config.keys():
        if config_name in param_d:
            config.update({config_name: param_d[config_name]})
            param_d.pop(config_name)

    subst = []
    for symbol, value in param_d.iteritems():
        tokens = scs_parser.parse_param_expresion(value)
        try:
            value = float(sympy.sympify(scs_parser.params2values(tokens, instance.paramsd),sympy.abc._clash))
        except ValueError:
            raise scs_errors.ScsAnalysisError("Passed subsitution for %s is not a number")

        subst.append((symbol, value))

    if config['fscale'] == 'log':
        fs = np.logspace(np.log10(float(config['fstart'])), np.log10(float(config['fstop'])), int(config['npoints']))
    elif config['fscale'] == 'linear':
        fs = np.linspace(float(config['fstart']), float(config['fstop']), int(config['npoints']))
    else:
        raise scs_errors.ScsAnalysisError(("Option %s for fscale invalid!" % config['yscale']))

    if config['yscale'] != 'log' and config['yscale'] != 'linear':
        raise scs_errors.ScsAnalysisError(("Option %s for yscale invalid!" % config['fscale']))

    filename = "%s.results" % file_sufix

    with open(filename, 'a') as fil:
        if config['xkcd'] == 'yes':
            plt.xkcd()
        plt.hold(True)

        for expresion in param_l:
            fil.write("%s: %s \n---------------------\n" % ('AC analysis of', expresion))
            tokens = scs_parser.parse_analysis_expresion(expresion)
            value0 = sympy.factor(sympy.sympify(scs_parser.results2values(tokens, instance),sympy.abc._clash), s).simplify()
            fil.write("%s = %s \n\n" % (expresion, str(value0)))
            denominator = sympy.denom(value0)
            numerator = sympy.numer(value0)
            poles = sympy.solve(denominator, s)
            zeros = sympy.solve(numerator, s)
            poles_r = sympy.roots(denominator, s)
            zeros_r = sympy.roots(numerator, s)
            gdc = str(value0.subs(s, 0).simplify())
            fil.write('G_DC = %s\n\n' % gdc)

            p = 0
            titled = 1
            for pole, degree in poles_r.iteritems():
                if pole == 0:
                    titled *= s ** degree
                else:
                    titled *= (s / sympy.symbols("\\omega_p%d" % p) + 1)
                p += 1
            z = 0
#.........这里部分代码省略.........

def ALGO4(As,Bs,Cs,Ds,do_test):
    
    
    
#-------------------STEP 1------------------------------
    if not is_row_proper(As):
        Us,Ast=row_proper(As)
        
    else:
        Us,Ast=eye(As.rows),As
    
    Bst=Us*Bs
    Bst=expand(Bst)
    r=Ast.cols
    #-------------------STEP 2------------------------------
    K=simplify(Ast.inv()*Bst)  #very important 
    Ys=zeros(K.shape)
    for i,j in  product(range(K.rows),range(K.cols)):
        Ys[i,j],q=div(numer(K[i,j]),denom(K[i,j]))       
    
    B_hat=Bst-Ast*Ys
    
    #-------------------END STEP 2------------------------------
    #-------------------STEP 3------------------------------
    Psi=diag(*[[s**( mc.row_degrees(Ast,s)[j]  -i -1) for i in range( mc.row_degrees(Ast,s)[j])] for j in range(r)]).T
    S=diag(*[s**(rho) for rho in mc.row_degrees(Ast,s)])
    Ahr=mc.highest_row_degree_matrix(Ast,s)
    Help=Ast-S*Ahr
    
    SOL={}
    numvar=Psi.rows*Psi.cols
    alr=symbols('a0:%d'%numvar)
    Alr=Matrix(Psi.cols,Psi.rows,alr)
    RHS=Psi*Alr
    for i,j in  product(range(Help.rows),range(Help.cols)):                 #diagonal explain later
        SOL.update(solve_undetermined_coeffs(Eq(Help[i,j],RHS[i,j]),alr,s))
    
    Alr=Alr.subs(SOL)    #substitute(SOL)
    
    
    Aoc=Matrix(BlockDiagMatrix(*[Matrix(rho, rho, lambda i,j: KroneckerDelta(i+1,j))for rho in mc.row_degrees(Ast,s)]))
    Boc=eye(sum(mc.row_degrees(Ast,s)))
    Coc=Matrix(BlockDiagMatrix(*[SparseMatrix(1,rho,{(x,0):1 if x==0 else 0 for x in range(rho)}) for rho in mc.row_degrees(Ast,s)]))

    A0=Aoc-Alr*Ahr.inv()*Matrix(Coc)
    C0=Ahr.inv()*Coc



    SOL={}
    numvar=Psi.cols*Bst.cols
    b0=symbols('b0:%d'%numvar)
    B0=Matrix(Psi.cols,Bst.cols,b0)
    RHS=Psi*B0

    for i,j in  product(range(B_hat.rows),range(B_hat.cols)):                 #diagonal explain later
        SOL.update(solve_undetermined_coeffs(Eq(B_hat[i,j],RHS[i,j]),b0,s))
    B0=B0.subs(SOL)    #substitute(SOL)

    LHS_matrix=simplify(Cs*C0)                                        #left hand side of the equation (1)

    sI_A=s*eye(A0.cols)- A0
    max_degree=mc.find_degree(LHS_matrix,s)                   #get the degree of the matrix at the LHS
                                                  #which is also the maximum degree for the coefficients of Λ(s)
    #---------------------------Creating Matrices Λ(s) and C -------------------------------------

    Lamda=[]
    numvar=((max_degree))*A0.cols
    a=symbols('a0:%d'%numvar)
    for i in range(A0.cols):                                        # paratirisi den douleuei to prin giat;i otra oxi diagonios
        p=sum(a[n +i*(max_degree)]*s**n for n in range(max_degree))  # we want variables one degree lower because we are multiplying by first order monomials
        Lamda.append(p)
        
    Lamda=Matrix(Cs.rows,A0.cols,Lamda)                                #convert the list to Matrix
    
    c=symbols('c0:%d'%(Lamda.rows*Lamda.cols))
    C=Matrix(Lamda.rows,Lamda.cols,c)
    #-----------------------------------------
    
    RHS_matrix=Lamda*sI_A +C                            #right hand side of the equation (1)
     
    '''
    -----------Converting equation (1) to a system of linear -----------
    -----------equations, comparing the coefficients of the  -----------
    -----------polynomials in both sides of the equation (1) -----------
    '''
     EQ=[Eq(LHS_matrix[i,j],expand(RHS_matrix[i,j])) for i,j in product(range(LHS_matrix.rows),range(LHS_matrix.cols)) ]

def guess(l, all=False, evaluate=True, niter=2, variables=None):
    """
    This function is adapted from the Rate.m package for Mathematica
    written by Christian Krattenthaler.
    It tries to guess a formula from a given sequence of rational numbers.

    In order to speed up the process, the 'all' variable is set to False by
    default, stopping the computation as some results are returned during an
    iteration; the variable can be set to True if more iterations are needed
    (other formulas may be found; however they may be equivalent to the first
    ones).

    Another option is the 'evaluate' variable (default is True); setting it
    to False will leave the involved products unevaluated.

    By default, the number of iterations is set to 2 but a greater value (up
    to len(l)-1) can be specified with the optional 'niter' variable.
    More and more convoluted results are found when the order of the
    iteration gets higher:

      * first iteration returns polynomial or rational functions;
      * second iteration returns products of rising factorials and their
        inverses;
      * third iteration returns products of products of rising factorials
        and their inverses;
      * etc.

    The returned formulas contain symbols i0, i1, i2, ... where the main
    variables is i0 (and auxiliary variables are i1, i2, ...). A list of
    other symbols can be provided in the 'variables' option; the length of
    the least should be the value of 'niter' (more is acceptable but only
    the first symbols will be used); in this case, the main variable will be
    the first symbol in the list.

    >>> from sympy.concrete.guess import guess
    >>> guess([1,2,6,24,120], evaluate=False)
    [Product(i1 + 1, (i1, 1, i0 - 1))]

    >>> from sympy import symbols
    >>> r = guess([1,2,7,42,429,7436,218348,10850216], niter=4)
    >>> i0 = symbols("i0")
    >>> [r[0].subs(i0,n).doit() for n in range(1,10)]
    [1, 2, 7, 42, 429, 7436, 218348, 10850216, 911835460]
    """
    if any(a==0 for a in l[:-1]):
        return []
    N = len(l)
    niter = min(N-1, niter)
    myprod = product if evaluate else Product
    g = []
    res = []
    if variables == None:
        symb = symbols('i:'+str(niter))
    else:
        symb = variables
    for k, s in enumerate(symb):
        g.append(l)
        n, r = len(l), []
        for i in range(n-2-1, -1, -1):
            ri = rinterp(enumerate(g[k][:-1], start=1), i, X=s)
            if ((denom(ri).subs({s:n}) != 0)
                    and (ri.subs({s:n}) - g[k][-1] == 0)
                    and ri not in r):
              r.append(ri)
        if r:
            for i in range(k-1, -1, -1):
                r = list(map(lambda v: g[i][0]
                      * myprod(v, (symb[i+1], 1, symb[i]-1)), r))
            if not all: return r
            res += r
        l = [Rational(l[i+1], l[i]) for i in range(N-k-1)]
    return res

 def export(self, x):
     return float(sp.numer(x)) / sp.denom(x)

        if G[r, c] != 0:
            m1.append(G[r, c])
            print G[r, c]

print '---Minors of order 2---'
m2_1 = G[:,[1,2]]
m2_2 = G[:,[0,2]]
m2_3 = G[:,[0,1]]
print m2_1
print m2_2
print m2_3

print '---All poles---'
m1count = np.shape(m1)[0]
for m in range(0, m1count):
    print sp.solve(sp.denom(m1[m]))

def M21(s):
    return G[:,[1,2]]
def M12(s):
    return G[:,[0,1]]
def M13(s):
    return G[:,[0,1]]
    
print poles(M21)
print poles(M12)
print poles(M13)

print 'Therefore the poles are -1, 1 and -2'

##Usefull Sage example

def apart_fact(Fs):
    S=0
    Ls = apart(Fs).as_ordered_terms()
    for l in Ls:
        S += numer(l)/factor(denom(l))
    return expand_mul(S)