def _integrify_basis( basis ):
    """
    Take a vector basis with sympy.Rational in the entries,
    and scale them so they are all integers by finding the least
    common multiple of the denominators and multiplying by that.
    """

    def _gcd(a, b):
	"""Return greatest common divisor using Euclid's Algorithm."""
	while b:
	    a, b = b, a % b
	return a

    def _lcm(a, b):
	"""Return lowest common multiple."""
	return a * b // _gcd(a, b)

    def _lcmm(*args):
	"""Return lcm of args."""   
	return reduce(_lcm, args)

    #This is intended for lists of type sympy.Rational
    assert( use_symbolic_math )

    new_basis = []
    for vec in basis:
        #Make a list of the denominators
        denominators = [sympy.fraction(e)[1] for e in vec]
        #Find the least common multiple
        least_common_multiple = _lcmm( *denominators )
        #Multiply all the entries by that, make it a python Fraction object
        new_vec = [ Fraction( int(e*least_common_multiple), 1 ) for e in vec]
        new_basis.append(new_vec)

    return new_basis

    def set_primitive_set(self):
        """
        Set up deap set of primitives, needed for genetic programming
        """
        self.pset = gp.PrimitiveSet(self.get_unique_str(), 1)
        self.pset.renameArguments(ARG0=str(self.var))

        # add these symbols into primitive set
        for S in self.series_to_sum.free_symbols - {self.var}:
            self.pset.addTerminal(S, name=str(S))

        # Add basic operations into the
        self.use_func(
            (operator.mul, 2), (operator.div, 2), (operator.sub, 2), (operator.add, 2)
        )

        # find unique number from series
        unique_num = set()
        for s_term in self.partitioned_series:
            unique_num.update(S for S in sympy.postorder_traversal(s_term) if S.is_Number)

        # convert numbers into fractions and extract nominator and denominator separately
        unique_num = itertools.chain(*(sympy.fraction(S) for S in unique_num))
        self.unique_num = sorted(set(unique_num))

        return self

    def fix_expressions(cc_num, common_denom_expr, lmatrix,
                        species_independent, species_dependent):

        fix_denom = SymcaToolBox.get_fix_denom(
            lmatrix,
            species_independent,
            species_dependent
        )
        fix = False
        cd_num, cd_denom = fraction(common_denom_expr)
        ret2 = cd_num


        if type(cc_num) is list:
            new_cc_num = cc_num[:]
        else:
            new_cc_num = cc_num[:, :]

        for i, each in enumerate(new_cc_num):
             new_cc_num[i] = ((each * cd_denom)).expand()

        for each in new_cc_num:
            for symb in fix_denom.atoms(Symbol):
                if symb in each.atoms(Symbol):
                    fix = True
                    break
            if fix: break

        if fix:
            for i, each in enumerate(new_cc_num):
                new_cc_num[i] = (each / fix_denom).expand()

            ret2 = (cd_num / fix_denom).expand()

        return new_cc_num, ret2

def eval_trigsubstitution(theta, func, rewritten, substep, restriction, integrand, symbol):
    func = func.subs(sympy.sec(theta), 1/sympy.cos(theta))

    trig_function = list(func.find(TrigonometricFunction))
    assert len(trig_function) == 1
    trig_function = trig_function[0]
    relation = sympy.solve(symbol - func, trig_function)
    assert len(relation) == 1
    numer, denom = sympy.fraction(relation[0])

    if isinstance(trig_function, sympy.sin):
        opposite = numer
        hypotenuse = denom
        adjacent = sympy.sqrt(denom**2 - numer**2)
        inverse = sympy.asin(relation[0])
    elif isinstance(trig_function, sympy.cos):
        adjacent = numer
        hypotenuse = denom
        opposite = sympy.sqrt(denom**2 - numer**2)
        inverse = sympy.acos(relation[0])
    elif isinstance(trig_function, sympy.tan):
        opposite = numer
        adjacent = denom
        hypotenuse = sympy.sqrt(denom**2 + numer**2)
        inverse = sympy.atan(relation[0])

    substitution = [
        (sympy.sin(theta), opposite/hypotenuse),
        (sympy.cos(theta), adjacent/hypotenuse),
        (sympy.tan(theta), opposite/adjacent),
        (theta, inverse)
    ]
    return sympy.Piecewise(
        (_manualintegrate(substep).subs(substitution).trigsimp(), restriction)
    )

def solve(a, b):
    s = sympy.fraction(sympy.Rational(a,b))
    if sympy.log(s[1],2) % 1 != 0:
        return 'impossible'
    y = sympy.ceiling(sympy.log(s[1]/s[0],2))
    if y > 40:
        return 'impossible'
    else:
        return str(y)

def transfer_function():
	t, s = symbols('t, s')
	f = ( 4*(s+2)*((s+Rational(2))**3) )/((s+6)*( (s+4)**2 ) )
	H = syslin( f )
	# disp( ilaplace(apart(H), s, t) )  # can't direct

	# http://docs.sympy.org/0.6.7/modules/simplify.html
	numer, _denom = fraction(H)

	disp(numer)
	print roots( numer )  # fixme: !!! complex roots!
	print roots( denom(H) )

    def make_internals_dict(cc_sol, cc_names, common_denom_expr, path_to):
        simpl_dic = {}
        for i, each in enumerate(cc_sol):
            expr = each / common_denom_expr
            expr = SymcaToolBox.maxima_factor(expr, path_to)
            num, denom = fraction(expr)
            if not simpl_dic.has_key(denom):
                simpl_dic[denom] = [[], []]
            simpl_dic[denom][0].append(cc_names[i])
            simpl_dic[denom][1].append(num)

        return simpl_dic

    def _print_Pow(self, power):
        from sympy import fraction
        b, e = power.as_base_exp()
        if power.is_commutative:
            if e is S.NegativeOne:
                return prettyForm("1")/self._print(b)
            n, d = fraction(e)
            if n is S.One and d.is_Atom and not e.is_Integer:
                return self._print_nth_root(b, e)
            if e.is_Rational and e < 0:
                return prettyForm("1")/self._print(b)**self._print(-e)

        # None of the above special forms, do a standard power
        return self._print(b)**self._print(e)

    def maxima_factor(expression, path_to):
        """
        This function is equivalent to the sympy.cancel()
        function but uses maxima instead
        """

        maxima_in_file = join(path_to,'in.txt').replace('\\','\\\\')
        maxima_out_file = join(path_to,'out.txt').replace('\\','\\\\')
        if expression.is_Matrix:
            expr_mat = expression[:, :]
            # print expr_mat
            print 'Simplifying matrix with ' + str(len(expr_mat)) + ' elements'
            for i, e in enumerate(expr_mat):

                sys.stdout.write('*')
                sys.stdout.flush()
                if (i + 1) % 50 == 0:
                    sys.stdout.write(' ' + str(i + 1) + '\n')
                    sys.stdout.flush()
                # print e
                expr_mat[i] = SymcaToolBox.maxima_factor(e, path_to)
            sys.stdout.write('\n')
            sys.stdout.flush()
            return expr_mat
        else:
            batch_string = (
                'stardisp:true;stringout("'
                + maxima_out_file + '",factor(' + str(expression) + '));')
            # print batch_string
            with open(maxima_in_file, 'w') as f:
                f.write(batch_string)

            config = ConfigReader.get_config()
            if config['platform'] == 'win32':
                maxima_command = [config['maxima_path'], '--batch=' + maxima_in_file]
            else:
                maxima_command = ['maxima', '--batch=' + maxima_in_file]

            dn = open(devnull, 'w')
            subprocess.call(maxima_command, stdin=dn, stdout=dn, stderr=dn)
            simplified_expression = ''

            with open(maxima_out_file) as f:
                for line in f:
                    if line != '\n':
                        simplified_expression = line[:-2]
            frac = fraction(sympify(simplified_expression))
            # print frac[0].expand()/frac[1].expand()
            return frac[0].expand() / frac[1].expand()

    def fix_expressions(cc_num, common_denom_expr, lmatrix, species_independent, species_dependent):

        fix_denom = SymcaToolBox.get_fix_denom(
            lmatrix,
            species_independent,
            species_dependent
        )
        #print fix_denom


        cd_num, cd_denom = fraction(common_denom_expr)

        new_cc_num = cc_num[:, :]
        #print type(new_cc_num)
        for i, each in enumerate(new_cc_num):
            new_cc_num[i] = ((each * cd_denom) / fix_denom).expand()

        return new_cc_num, (cd_num / fix_denom).expand()

    def _eval_expand_mul(self, deep=True, **hints):
        from sympy import fraction

        expr = self
        n, d = fraction(expr)
        if d is not S.One:
            expr = n / d._eval_expand_mul(deep=deep, **hints)
            if not expr.is_Mul:
                return expr._eval_expand_mul(deep=deep, **hints)

        plain, sums, rewrite = [], [], False
        for factor in expr.args:
            if deep:
                term = factor.expand(deep=deep, **hints)
                if term != factor:
                    factor = term
                    rewrite = True

            if factor.is_Add:
                sums.append(factor)
                rewrite = True
            else:
                if factor.is_commutative:
                    plain.append(factor)
                else:
                    Wrapper = Basic
                    sums.append(Wrapper(factor))

        if not rewrite:
            return expr
        else:
            plain = Mul(*plain)
            if sums:
                terms = Mul._expandsums(sums)
                args = []
                for term in terms:
                    t = Mul(plain, term)
                    if t.is_Mul and any(a.is_Add for a in t.args):
                        t = t._eval_expand_mul(deep=deep)
                    args.append(t)
                return Add(*args)
            else:
                return plain

    def _eval_expand_mul(self, **hints):
        from sympy import fraction, expand_mul

        # Handle things like 1/(x*(x + 1)), which are automatically converted
        # to 1/x*1/(x + 1)
        expr = self
        n, d = fraction(expr)
        if d.is_Mul:
            expr = n/d._eval_expand_mul(**hints)
            if not expr.is_Mul:
                return expand_mul(expr, deep=False)

        plain, sums, rewrite = [], [], False
        for factor in expr.args:
            if factor.is_Add:
                sums.append(factor)
                rewrite = True
            else:
                if factor.is_commutative:
                    plain.append(factor)
                else:
                    sums.append(Basic(factor))  # Wrapper

        if not rewrite:
            return expr
        else:
            plain = Mul(*plain)
            if sums:
                terms = Mul._expandsums(sums)
                args = []
                for term in terms:
                    t = Mul(plain, term)
                    if t.is_Mul and any(a.is_Add for a in t.args):
                        t = t._eval_expand_mul()
                    args.append(t)
                return Add(*args)
            else:
                return plain

def get_roots(expr):
    num, den = sympy.fraction(expr)
    return sympy.solve(den, s), sympy.solve(num, s)

def get_roots(expr):
    """Given the transfer function ``expr``, returns ``poles, zeros``.
    """
    num, den = sympy.fraction(expr)
    return sympy.solve(den, s), sympy.solve(num, s)

def test_radsimp():
    r2=sqrt(2)
    r3=sqrt(3)
    r5=sqrt(5)
    r7=sqrt(7)
    assert radsimp(1/r2) == \
        sqrt(2)/2
    assert radsimp(1/(1 + r2)) == \
        -1 + sqrt(2)
    assert radsimp(1/(r2 + r3)) == \
        -sqrt(2) + sqrt(3)
    assert fraction(radsimp(1/(1 + r2 + r3))) == \
        (-sqrt(6) + sqrt(2) + 2, 4)
    assert fraction(radsimp(1/(r2 + r3 + r5))) == \
        (-sqrt(30) + 2*sqrt(3) + 3*sqrt(2), 12)
    assert fraction(radsimp(1/(1 + r2 + r3 + r5))) == \
        (-34*sqrt(10) -
        26*sqrt(15) -
        55*sqrt(3) -
        61*sqrt(2) +
        14*sqrt(30) +
        93 +
        46*sqrt(6) +
        53*sqrt(5), 71)
    assert fraction(radsimp(1/(r2 + r3 + r5 + r7))) == \
        (-50*sqrt(42) - 133*sqrt(5) - 34*sqrt(70) -
        145*sqrt(3) + 22*sqrt(105) + 185*sqrt(2) +
        62*sqrt(30) + 135*sqrt(7), 215)
    z = radsimp(1/(1 + r2/3 + r3/5 + r5 + r7))
    assert len((3616791619821680643598*z).args) == 16
    assert radsimp(1/z) == 1/z
    assert radsimp(1/z, max_terms=20).expand() == 1 + r2/3 + r3/5 + r5 + r7
    assert radsimp(1/(r2*3)) == \
        sqrt(2)/6
    assert radsimp(1/(r2*a + r3 + r5 + r7)) == 1/(r2*a + r3 + r5 + r7)
    assert radsimp(1/(r2*a + r2*b + r3 + r7)) == \
        ((sqrt(42)*(a + b) +
        sqrt(3)*(-a**2 - 2*a*b - b**2 - 2)  +
        sqrt(7)*(-a**2 - 2*a*b - b**2 + 2)  +
        sqrt(2)*(a**3 + 3*a**2*b + 3*a*b**2 - 5*a + b**3 - 5*b))/
        ((a**4 + 4*a**3*b + 6*a**2*b**2 - 10*a**2  +
        4*a*b**3 - 20*a*b + b**4 - 10*b**2 + 4)))/2
    assert radsimp(1/(r2*a + r2*b + r2*c + r2*d)) == \
        (sqrt(2)/(a + b + c + d))/2
    assert radsimp(1/(1 + r2*a + r2*b + r2*c + r2*d)) == \
        ((sqrt(2)*(-a - b - c - d) + 1)/
        (-2*a**2 - 4*a*b - 4*a*c - 4*a*d - 2*b**2 -
        4*b*c - 4*b*d - 2*c**2 - 4*c*d - 2*d**2 + 1))
    assert radsimp((y**2 - x)/(y - sqrt(x))) == \
        sqrt(x) + y
    assert radsimp(-(y**2 - x)/(y - sqrt(x))) == \
        -(sqrt(x) + y)
    assert radsimp(1/(1 - I + a*I)) == \
        (I*(-a + 1) + 1)/(a**2 - 2*a + 2)
    assert radsimp(1/((-x + y)*(x - sqrt(y)))) == (x + sqrt(y))/((-x + y)*(x**2 - y))
    e = (3 + 3*sqrt(2))*x*(3*x - 3*sqrt(y))
    assert radsimp(e) == 9*x*(1 + sqrt(2))*(x - sqrt(y))
    assert radsimp(1/e) == (-1 + sqrt(2))*(x + sqrt(y))/(9*x*(x**2 - y))
    assert radsimp(1 + 1/(1 + sqrt(3))) == Mul(S(1)/2, 1 + sqrt(3), evaluate=False)
    A = symbols("A", commutative=False)
    assert radsimp(x**2 + sqrt(2)*x**2 - sqrt(2)*x*A) == x**2 + sqrt(2)*(x**2 - x*A)
    assert radsimp(1/sqrt(5 + 2 * sqrt(6))) == -sqrt(2) + sqrt(3)
    assert radsimp(1/sqrt(5 + 2 * sqrt(6))**3) == -11*sqrt(2) + 9*sqrt(3)

    # coverage not provided by above tests
    assert collect_const(2*sqrt(3) + 4*a*sqrt(5)) == Mul(2, (2*sqrt(5)*a + sqrt(3)), evaluate=False)
    assert collect_const(2*sqrt(3) + 4*a*sqrt(5), sqrt(3)) == 2*(2*sqrt(5)*a + sqrt(3))
    assert collect_const(sqrt(2)*(1 + sqrt(2)) + sqrt(3) + x*sqrt(2)) == \
        sqrt(2)*(x + 1 + sqrt(2)) + sqrt(3)

def test_fraction():
    x, y, z = map(Symbol, 'xyz')

    assert fraction(Rational(1, 2)) == (1, 2)

    assert fraction(x) == (x, 1)
    assert fraction(1/x) == (1, x)
    assert fraction(x/y) == (x, y)
    assert fraction(x/2) == (x, 2)

    assert fraction(x*y/z) == (x*y, z)
    assert fraction(x/(y*z)) == (x, y*z)

    assert fraction(1/y**2) == (1, y**2)
    assert fraction(x/y**2) == (x, y**2)

    assert fraction((x**2+1)/y) == (x**2+1, y)
    assert fraction(x*(y+1)/y**7) == (x*(y+1), y**7)

    assert fraction(exp(-x), exact=True) == (exp(-x), 1)

L=getLagrangian(M,m2,gamma,alpha1,alpha2,psi,*A)

'''L=L.subs({A[0]:0, A[3]:0}).doit()
Axeqm=fieldEqn(L,A[2],M.x)
Axeqm=Axeqm.subs(A[2],d).series(d,n=2).subs(sp.Order(d**2),0).subs(d,A[2])
Axeqm=(Axeqm/(M.x[0]**2*(M.x[0]**3-1))).ratsimp().collect(A[2])
Axeqm=sum(sp.simplify(e,ratio=1).collect(A[1]) for e in Axeqm.args)
#Axeqm=Axeqm.subs({A[2]:d,A[2].diff(M.x[0]):d2,A[2].diff(M.x[0]):d3}).series(d)
sp.pprint(Axeqm)
totex.save(Axeqm.subs({psi:sp.Symbol('psi'),A[1]:sp.Symbol('phi'),A[2]:sp.Symbol('A_x')}),'Axeqm')'''

L=L.subs({gamma:0,alpha1:0,alpha2:0,A[0]:0,A[3]:0}).doit()
#psieqm=simpSum( sp.fraction((fieldEqn(L,psi,M.x)/2).together())[0].collect(psi) )
#psieqm=(fieldEqn(L,psi,M.x)/(-M.r**2*2*M.f)).expand().collect(psi)

psieqm, phieqm, Axeqm=[sp.fraction(fieldEqn(L,f,M.x).together())[0].expand().collect(f) for f in [psi,A[1],A[2]]]

t=sp.Dummy('t')
hpsieqm=psieqm.subs(M.x[0],M.zh+t).simplify()

hpsiexp=[hpsieqm.diff(t,deg).limit(t,0).subs(A[1].subs(M.x[0],M.zh),0) for deg in range(2)]
sp.pprint( sp.solve(hpsiexp[1],psi.subs(M.x[0],M.zh))) #.diff(M.x[0]).subs(M.x[0],M.zh)) )
#phieqm=simpSum( (phieqm/(phieqm.subs(A[1].diff(M.x[0],2),d).diff(d))).together().ratsimp().collect(A[1]) )
#sol=sp.solve(phieqm.subs(A[1].diff(M.x[0],2),d),d)
#sp.pprint(sol)
#assert len(sol)==1
#fr=sp.fraction(sol[0])
#sp.pprint(fr[0].collect(A[1])/fr[1].collect(A[1]))
if __name__=='__main__':
    sp.pprint(psieqm)
    sp.pprint(phieqm)

def sym_rpoly(na=None, nb=None, nk=None, nterms=None):
    """Compute the R polynomial as defined above.

    Internally the construction of R is done symbolically. If the numerical
    variables (na, nb, nk) were supplied, these values are substituted at the
    end and a sympy.Poly object with numerical coefficients is returned.  If
    (na, nb, nk) are not given, then nterms *must* be given, and a symbolic
    answer is returned.

    Parameters
    ----------
    na : ndarray, optional
      Numerical array of 'a' coefficients.

    nb : ndarray, optional
      Numerical array of 'b' coefficients.

    k : float, optional
      Numerical value of k.

    nterms : int, optional.
      Number of terms. This is only used if na, nb and nk are *not* given, in
      which case a symbolic answer is returned with nterms total.

    Returns
    -------
    poly : sympy.Poly instance
      A univariate polynomial in x.

    Examples
    --------
    With only nterms, a symbolic polynomial is returned:
    >>> sym_rpoly(nterms=1)
    Poly(x**2 - k*x - b_0*k/a_0, x)

    But if numerical values are supplied, the output polynomial has numerical
    values:
    >>> sym_rpoly([2], [5], 3)
    Poly(x**2 - 3*x - 15/2, x)
    >>> sym_rpoly([1, 2], [4, 6], 1)
    Poly(x**3 + 3*x**2 - 7/2*x - 6, x)
    """
    if na is None and nterms is None:
        raise ValueError('You must provide either arrays or nterms')

    if nterms is None:
        na = np.asarray(na)
        nb = np.asarray(nb)
        # We have input arrays, return numerical answer
        if na.ndim > 1:
            raise ValueError('Input arrays must be one-dimensional')
        nterms = na.shape
        mode = 'num'
    else:
        # We have only a length, return symbolic answer.
        mode = 'sym'
        
    # Create symbolic variables
    k, x = sym.symbols('k x')
    a = symarray(nterms, 'a')
    b = symarray(nterms, 'b')

    # Construct polynomial symbolically
    t = [ ai/(ai*x+bi) for (ai, bi) in zip(a,b)]
    P, Q = sym.fraction(sym.together(sum(t)))
    Rs = sym.Poly(x**2*P -k*Q, x)

    if mode == 'num':
        # Substitute coefficients for numerical values
        Rs = Rs.subs([(k, nk)] + zip(a, na) + zip(b, nb))

    # Normalize
    Rs /= Rs.lead_coeff
    return Rs

## DC Gain
dc_gain = simplify(res_simp).limit('s',0)
print ""
print "DC Gain:"
pprint(dc_gain)
print ""

## AC Transfer function
tf = collect(res_simp,s)
print "AC Transfer function:"
pprint(tf)
print ""

## Denominator of transfer function
tf_denom = fraction(tf)[1]
tf_denom = tf_denom.expand()
print "Denominator of transfer function times R[1]:"
pprint(tf_denom)
print ""

## Poles
tf_poles = solve(tf_denom,s) 
tf_poles = simplify(tf_poles)
print "Poles of transfer function:"
print tf_poles
p1 = tf_poles # Need to exctract one pole
p2 = tf_poles # Need to extract one pole
pprint(p1)
pprint(p2)
print ""

def eqSimp(e):
    return sp.fraction(e.together())[0]