def test_ratsimp():
    x = Symbol("x")
    y = Symbol("y")
    e = 1/x+1/y
    assert e != (x+y)/(x*y)
    assert ratsimp(e) == (x+y)/(x*y)

    e = 1/(1+1/x)
    assert ratsimp(e) == x/(x+1)
    assert ratsimp(exp(e)) == exp(x/(x+1))

def test_pmint_logexp():
    f = (1 + x + x*exp(x))*(x + log(x) + exp(x) - 1)/(x + log(x) + exp(x))**2/x
    g = log(x**2 + 2*x*exp(x) + 2*x*log(x) + exp(2*x) + 2*exp(x)*log(x) + log(x)**2)/2 + 1/(x + exp(x) + log(x))

    # TODO: Optimal solution is g = 1/(x + log(x) + exp(x)) + log(x + log(x) + exp(x)),
    # but SymPy requires a lot of guidance to properly simplify heurisch() output.

    assert ratsimp(heurisch(f, x)) == g

def test_pmint_logexp():
    if ON_TRAVIS:
        # See https://github.com/sympy/sympy/pull/12795
        skip("Too slow for travis.")

    f = (1 + x + x*exp(x))*(x + log(x) + exp(x) - 1)/(x + log(x) + exp(x))**2/x
    g = log(x + exp(x) + log(x)) + 1/(x + exp(x) + log(x))

    assert ratsimp(heurisch(f, x)) == g

 def capacitance_matrix_variables(self, symbolic=False):
     """
     Calculates the capacitance matrix for the energy term of the qubit Lagrangian in the variable respresentation.
     """                        
     
     if symbolic:
         C = self.linear_coordinate_transform.T*self.capacitance_matrix(symbolic)*self.linear_coordinate_transform
         C = sympy.Matrix([sympy.nsimplify(sympy.ratsimp(x)) for x in C]).reshape(*(C.shape))
     else:
         C = np.einsum('ji,jk,kl->il', self.linear_coordinate_transform,self.capacitance_matrix(symbolic),self.linear_coordinate_transform)
     return C

def test_trigsimp1():
    x, y = symbols('x,y')

    assert trigsimp(1 - sin(x)**2) == cos(x)**2
    assert trigsimp(1 - cos(x)**2) == sin(x)**2
    assert trigsimp(sin(x)**2 + cos(x)**2) == 1
    assert trigsimp(1 + tan(x)**2) == 1/cos(x)**2
    assert trigsimp(1/cos(x)**2 - 1) == tan(x)**2
    assert trigsimp(1/cos(x)**2 - tan(x)**2) == 1
    assert trigsimp(1 + cot(x)**2) == 1/sin(x)**2
    assert trigsimp(1/sin(x)**2 - 1) == cot(x)**2
    assert trigsimp(1/sin(x)**2 - cot(x)**2) == 1

    assert trigsimp(5*cos(x)**2 + 5*sin(x)**2) == 5
    assert trigsimp(5*cos(x/2)**2 + 2*sin(x/2)**2) in \
                [2 + 3*cos(x/2)**2, 5 - 3*sin(x/2)**2]

    assert trigsimp(sin(x)/cos(x)) == tan(x)
    assert trigsimp(2*tan(x)*cos(x)) == 2*sin(x)
    assert trigsimp(cot(x)**3*sin(x)**3) == cos(x)**3
    assert trigsimp(y*tan(x)**2/sin(x)**2) == y/cos(x)**2
    assert trigsimp(cot(x)/cos(x)) == 1/sin(x)

    assert trigsimp(sin(x + y) + sin(x - y)) == 2*sin(x)*cos(y)
    assert trigsimp(sin(x + y) - sin(x - y)) == 2*sin(y)*cos(x)
    assert trigsimp(cos(x + y) + cos(x - y)) == 2*cos(x)*cos(y)
    assert trigsimp(cos(x + y) - cos(x - y)) == -2*sin(x)*sin(y)
    assert ratsimp(trigsimp(tan(x + y) - tan(x)/(1 - tan(x)*tan(y)))) == \
        -tan(y)/(tan(x)*tan(y) -1)

    assert trigsimp(sinh(x + y) + sinh(x - y)) == 2*sinh(x)*cosh(y)
    assert trigsimp(sinh(x + y) - sinh(x - y)) == 2*sinh(y)*cosh(x)
    assert trigsimp(cosh(x + y) + cosh(x - y)) == 2*cosh(x)*cosh(y)
    assert trigsimp(cosh(x + y) - cosh(x - y)) == 2*sinh(x)*sinh(y)
    assert ratsimp(trigsimp(tanh(x + y) - tanh(x)/(1 + tanh(x)*tanh(y)))) == \
        tanh(y)/(tanh(x)*tanh(y) + 1)

    assert trigsimp(cos(0.12345)**2 + sin(0.12345)**2) == 1
    e = 2*sin(x)**2 + 2*cos(x)**2
    assert trigsimp(log(e), deep=True) == log(2)

def test_ratsimp():
    f, g = 1/x + 1/y, (x + y)/(x*y)

    assert f != g and ratsimp(f) == g

    f, g = 1/(1 + 1/x), 1 - 1/(x + 1)

    assert f != g and ratsimp(f) == g

    f, g = x/(x + y) + y/(x + y), 1

    assert f != g and ratsimp(f) == g

    f, g = -x - y - y**2/(x + y) + x**2/(x + y), -2*y

    assert f != g and ratsimp(f) == g

    f = (a*c*x*y + a*c*z - b*d*x*y - b*d*z - b*t*x*y - b*t*x - b*t*z + e*x)/(x*y + z)
    G = [a*c - b*d - b*t + (-b*t*x + e*x)/(x*y + z),
         a*c - b*d - b*t - ( b*t*x - e*x)/(x*y + z)]

    assert f != g and ratsimp(f) in G

    A = sqrt(pi)

    B = log(erf(x) - 1)
    C = log(erf(x) + 1)

    D = 8 - 8*erf(x)

    f = A*B/D - A*C/D + A*C*erf(x)/D - A*B*erf(x)/D + 2*A/D

    assert ratsimp(f) == A*B/8 - A*C/8 - A/(4*erf(x) - 4)

def test_pmint_rat():
    # TODO: heurisch() is off by a constant: -3/4. Possibly different permutation
    # would give the optimal result?

    def drop_const(expr, x):
        if expr.is_Add:
            return Add(*[ arg for arg in expr.args if arg.has(x) ])
        else:
            return expr

    f = (x**7 - 24*x**4 - 4*x**2 + 8*x - 8)/(x**8 + 6*x**6 + 12*x**4 + 8*x**2)
    g = (4 + 8*x**2 + 6*x + 3*x**3)/(x**5 + 4*x**3 + 4*x) + log(x)

    assert drop_const(ratsimp(heurisch(f, x)), x) == g

def test_action_verbs():
    assert nsimplify((1/(exp(3*pi*x/5)+1))) == (1/(exp(3*pi*x/5)+1)).nsimplify()
    assert ratsimp(1/x + 1/y) == (1/x + 1/y).ratsimp()
    assert trigsimp(log(x), deep=True) == (log(x)).trigsimp(deep = True)
    assert radsimp(1/(2+sqrt(2))) == (1/(2+sqrt(2))).radsimp()
    assert powsimp(x**y*x**z*y**z, combine='all') == (x**y*x**z*y**z).powsimp(combine='all')
    assert simplify(x**y*x**z*y**z) == (x**y*x**z*y**z).simplify()
    assert together(1/x + 1/y) == (1/x + 1/y).together()
    assert separate((x*(y*z)**3)**2) == ((x*(y*z)**3)**2).separate()
    assert collect(a*x**2 + b*x**2 + a*x - b*x + c, x) == (a*x**2 + b*x**2 + a*x - b*x + c).collect(x)
    assert apart(y/(y+2)/(y+1), y) == (y/(y+2)/(y+1)).apart(y)
    assert combsimp(y/(x+2)/(x+1)) == (y/(x+2)/(x+1)).combsimp()
    assert factor(x**2+5*x+6) == (x**2+5*x+6).factor()
    assert refine(sqrt(x**2)) == sqrt(x**2).refine()
    assert cancel((x**2+5*x+6)/(x+2)) == ((x**2+5*x+6)/(x+2)).cancel()

 def print_Div(self, rule):
     with self.new_step():
         f, g = rule.numerator, rule.denominator
         fp, gp = f.diff(rule.symbol), g.diff(rule.symbol)
         x = rule.symbol
         ff = sympy.Function("f")(x)
         gg = sympy.Function("g")(x)
         qrule_left = sympy.Derivative(ff / gg, rule.symbol)
         qrule_right = sympy.ratsimp(sympy.diff(sympy.Function("f")(x) /
                                                sympy.Function("g")(x)))
         qrule = sympy.Eq(qrule_left, qrule_right)
         self.append("Apply the quotient rule, which is:")
         self.append(self.format_math_display(qrule))
         self.append("{} and {}.".format(self.format_math(sympy.Eq(ff, f)),
                                         self.format_math(sympy.Eq(gg, g))))
         self.append("To find {}:".format(self.format_math(ff.diff(rule.symbol))))
         with self.new_level():
             self.print_rule(rule.numerstep)
         self.append("To find {}:".format(self.format_math(gg.diff(rule.symbol))))
         with self.new_level():
             self.print_rule(rule.denomstep)
         self.append("Now plug in to the quotient rule:")
         self.append(self.format_math(diff(rule)))

def test_ratsimp_X1():
    e = -x-y-(x+y)**(-1)*y**2+(x+y)**(-1)*x**2
    assert e != -2*y
    assert ratsimp(e) == -2*y

def test_ratsimp():
    assert ratsimp(A*B - B*A) == A*B - B*A

             formula = recursion(m,n,base)
             formula = (formula.subs([(alpha,1/a),(beta,1/b)]))
             functions.append(("%s_%d_%d"%(name,m,n), formula))
             # Analytical expressions are singular for alpha=beta.
             # When alpha \approx beta, we therefore compute the expression
             # by writing beta=alpha+epsilon and writing a Taylor expansion
             # in epsilon.
             if m==0 and n==0:
                 # The Taylor expansion in epsilon is computed only for the
                 # base integral...
                 form = base.subs(beta,alpha+epsilon)
                 series_exp = 1
                 for itaylor in xrange(1,4*mmax+2):
                     series_exp += (-epsilon*R)**itaylor/factorial(itaylor)
                 form = form.subs(exp(-R*(alpha+epsilon)),exp(-R*alpha)*series_exp)
                 form = ratsimp(form)
                 taylor_expansion = simplify(form.subs(epsilon,0))
                 for itaylor in xrange(4*mmax+1):
                     form = ratsimp(diff(form,epsilon))
                     taylor_expansion += simplify(form.subs(epsilon, 0)*epsilon**(itaylor+1)/factorial(itaylor+1))
                 #print taylor_expansion
             # ...other Taylor expansions are derived using a modified recursion scheme.
             formula_taylor = simplify(recursiontaylor(m,n,taylor_expansion))
             functions.append(("%s_taylor_%d_%d"%(name,m,n), formula_taylor.subs([(alpha,1/a),(epsilon,b)])))
 # Write routines in Fortran source code
 print "Writing source code..."
 codegen(functions, "F95", "../src/libslater",to_files=True,
      argument_sequence=(a,b,R), project='libslater')
 # Write python wrapper
 with open('../src/libslater.py','w') as f:
     f.write('''#!/usr/bin/env python\n\n''')

def test_pmint_erf():
    f = exp(-x**2)*erf(x)/(erf(x)**3 - erf(x)**2 - erf(x) + 1)
    g = sqrt(pi)*log(erf(x) - 1)/8 - sqrt(pi)*log(erf(x) + 1)/8 - sqrt(pi)/(4*erf(x) - 4)

    assert ratsimp(heurisch(f, x)) == g

def integral_basis(f,x,y):
    """
    Compute the integral basis {b1, ..., bg} of the algebraic function
    field C[x,y] / (f).
    """
    # If the curve is not monic then map y |-> y/lc(x) where lc(x)
    # is the leading coefficient of f
    T  = sympy.Dummy('T')
    d  = sympy.degree(f,y)
    lc = sympy.LC(f,y)
    if x in lc:
        f = sympy.ratsimp( f.subs(y,y/lc)*lc**(d-1) )
    else:
        f = f/lc
        lc = 1

    #
    # Compute df
    #
    p = sympy.Poly(f,[x,y])
    n = p.degree(y)
    res = sympy.resultant(p,p.diff(y),y)
    factors = sympy.factor_list(res)[1]
    df = [k for k,deg in factors if (deg > 1) and (sympy.LC(k) == 1)]

    #
    # Compute series truncations at appropriate x points
    #
    alpha = []
    r = []
    for l in range(len(df)):
        k = df[l]
        alphak = sympy.roots(k).keys()[0]
        rk = compute_series_truncations(f,x,y,alphak,T)
        alpha.append(alphak)
        r.append(rk)

        
    #
    # Main Loop
    #
    a = [sympy.Dummy('a%d'%k) for k in xrange(n)]
    b = [1]
    for d in range(1,n):
        bd = y*b[-1]

        for l in range(len(df)):
            k = df[l]
            alphak = alpha[l]
            rk = r[l]

            found_something = True            
            while found_something:
                # construct system of equations consisting of the coefficients
                # of negative powers of (x-alphak) in the substitutions
                # A(r_{k,1}),...,A(r_{k,n})
                A  = (sum(ak*bk for ak,bk in zip(a,b)) + bd) / (x - alphak)

                coeffs = []
                for rki in rk:
                    # substitute and extract coefficients
                    A_rki  = A.subs(y,rki)
                    coeffs.extend(_negative_power_coeffs(A_rki, x, alphak))
               
                # solve the coefficient equations for a0,...,a_{d-1}
                coeffs = [coeff.as_numer_denom()[0] for coeff in coeffs]
                sols = sympy.solve_poly_system(coeffs, a[:d])
                if sols is None or sols == []:
                    found_something = False
                else:
                    sol  = sols[0]
                    bdm1 = sum( sol[i]*bk for i,bk in zip(range(d),b) )
                    bd   = (bdm1 + bd) / k
                        
        # bd found. Append to list of basis elements
        b.append( bd )

    # finally, convert back to singularized curve if necessary
    for i in xrange(1,len(b)):
        b[i] = b[i].subs(y,y*lc).ratsimp()

    return b

def test_pmint_logexp():
    f = (1 + x + x*exp(x))*(x + log(x) + exp(x) - 1)/(x + log(x) + exp(x))**2/x
    g = log(x + exp(x) + log(x)) + 1/(x + exp(x) + log(x))

    assert ratsimp(heurisch(f, x)) == g

n=500
x0,x1,y0,y1=sp.symbols('x0 x1 y0 y1',real=True)
sx=[x0,x1]
sy=[y0,y1]
#Kx=1
#Py=1
Kx=x1**2
Py=x1**2
#Phi1NL=-0.25/sp.pi*sp.log((x0-y0)**2+(x1-y1)**2)
Phi1NL = sp.log(((x0-y0)**2+(x1+y1)**2)/((x0-y0)**2+(x1-y1)**2)) / (sp.pi*2*x1*y1)
omgNL=[ Py*sp.diff(Phi1NL,sy[j]) for j in range(2) ]
#Phi2NL=-0.5/sp.pi*sp.atan((x1-y1)/(x0-y0))
Phi2NL = sp.log(((x0-y0)**2+(x1+y1)**2)/((x0-y0)**2+(x1-y1)**2)) * (x0-y0)/(2*sp.pi*x1)
+ (sp.atan((x1-y1)/(x0-y0)) - sp.atan((x1+y1)/(x0-y0))) / (2*sp.pi)

V2NL=[sp.ratsimp(Kx*sp.diff(Phi2NL,sx[j])) for j in range(2)]
VSNL=[sp.ratsimp(Kx*sp.diff(Phi1NL,sx[j])) for j in range(2)]

def Normal(x):
    n=x[0].shape[0]-1
    return np.array([-x[1,1:]+x[1,:n],x[0,1:]-x[0,:n]])

def normal(x):
    nr=Normal(x)
    return nr/np.sqrt(nr[0]**2+nr[1]**2)

def withoutZero(sFunct):
    def wrapped(sExpress,M,minDist=myeps):
        gd=np.sum([(M[i,1]-M[i,0])**2 for i in range(2)],axis=0)>minDist**2
        AA=M[0,0].copy()
        AA[np.logical_not(gd)]+=1

 def emit(name, func):
     for x in sorted(exparg):
         test(name, x, sp.ratsimp(sp.radsimp(func(x).rewrite(sp.exp))),
              no_trigh=True)

#!/usr/bin/python
import numpy as np
import sympy as sp

myeps=0.000001
x0,x1,y0,y1=sp.symbols('x0 x1 y0 y1',real=True)
sx=[x0,x1]
sy=[y0,y1]
_Kx_=1
_Py_=1
_Phi1_=-0.25/sp.pi*sp.log((x0-y0)**2+(x1-y1)**2)
_sOmega_=[_Py_*sp.diff(_Phi1_,sy[j]) for j in range(2)]
_Phi2_=-0.5/sp.pi*sp.atan((x1-y1)/(x0-y0))
_sV2_=[sp.ratsimp(_Kx_*sp.diff(_Phi2_,sx[j])) for j in range(2)]
_sVS_=[_Kx_*sp.diff(_Phi1_,sx[j]) for j in range(2)]

def Normal(x):
    n=x[0].shape[0]-1
    return np.array([-x[1,1:]+x[1,:n],x[0,1:]-x[0,:n]])

def normal(x):
    nr=Normal(x)
    return nr/np.sqrt(nr[0]**2+nr[1]**2)

def withoutZero(sFunct):
    def wrapped(sExpress,M,minDist=myeps):
        gd=np.sum([(M[i,1]-M[i,0])**2 for i in range(2)],axis=0)>minDist**2
        AA=M[0,0].copy()
        AA[np.logical_not(gd)]+=1
        result=np.zeros_like(M[0,0])
        result[gd]=sFunct(sExpress,np.array([[AA,M[0,1]],[M[1,0],M[1,1]]]),minDist)[gd]

def test_ratsimp_X2():
    e = x/(x+y)+y/(x+y)
    assert e != 1
    assert ratsimp(e) == 1

def test_ratsimp2():
    x = Symbol("x")
    e = 1/(1+1/x)
    assert (x+1)*ratsimp(e)/x == 1