def test_issue841():
    from sympy import expand_mul
    from sympy.abc import k

    assert expand_mul(integrate(exp(-x ** 2) * exp(I * k * x), (x, -oo, oo))) == sqrt(pi) * exp(-k ** 2 / 4)
    a, d = symbols("a d", positive=True)
    assert expand_mul(integrate(exp(-a * x ** 2 + 2 * d * x), (x, -oo, oo))) == sqrt(pi) * exp(d ** 2 / a) / sqrt(a)

def test_integrate_omit_var():
    y = Symbol('y')

    assert integrate(x) == x**2/2

    raises(ValueError, lambda: integrate(2))
    raises(ValueError, lambda: integrate(x*y))

def test_integrate_linearterm_pow():
    # check integrate((a*x+b)^c, x)  --  issue 3499
    y = Symbol('y', positive=True)
    # TODO: Remove conds='none' below, let the assumption take care of it.
    assert integrate(x**y, x, conds='none') == x**(y + 1)/(y + 1)
    assert integrate((exp(y)*x + 1/y)**(1 + sin(y)), x, conds='none') == \
        exp(-y)*(exp(y)*x + 1/y)**(2 + sin(y)) / (2 + sin(y))

	def setThreshold(self):
		import sympy
		# f(x) is probability density of processing times 
		# Using Bounded Pareto,
			# f(x) = {0 if x < 0, main.customEquation if x >= 0}
		# Probability jobs will be in class 1, ie. their processing time will be between Lower and Threshold
		# Prob(Lower <= procTime <= Threshold) = integral{from L to T}f(x)dx >> 0.8 (should be very close to 1)
		# Most jobs will be small, only a few will be big

		# ???
		# lambda * integral{from L to U} x*f(x) dx < numServers

		#HOW TO FIND THRESHOLD
		# integral{from L to T} x*f(x)dx / integral{from L to U} x*f(x)dx = 0.8 (becasue 80% of jobs will be small)
		# solve for T
		x, L, U, T, alpha = sympy.symbols('x L U T alpha', real=True)
		
		numerator = sympy.integrate((x*BoundedParetoDist.Function), (x, L, T))
		denominator = sympy.integrate((x*BoundedParetoDist.Function), (x, L, U))
		expected = numerator/denominator

		#expected = expected.subs(dict(alpha=JobClass.BPArray[0], L=JobClass.BPArray[1], U=JobClass.BPArray[2]))
		expected = expected.subs([(alpha, JobClass.BPArray[0]), (L, JobClass.BPArray[1]), (U, JobClass.BPArray[2])])
		#print "expected subed"
		#print expected
		
		#Equate expected = 0.8, solve for T
		#thresholdRange = sympy.solve(sympy.Eq(expected, 0.8), T)
		#MachineClass.Threshold = thresholdRange[1] #set only positive number

		##FORCE THRESHOLD AS LOGICAL VALUE
		MachineClass.Threshold = 800000
		GUI.writeToConsole(self.master, "Class threshold = %s"%MachineClass.Threshold)

def test_issue_8368():
    assert integrate(exp(-s*x)*cosh(x), (x, 0, oo)) == \
        Piecewise(
            (   pi*Piecewise(
                    (   -s/(pi*(-s**2 + 1)),
                        Abs(s**2) < 1),
                    (   1/(pi*s*(1 - 1/s**2)),
                        Abs(s**(-2)) < 1),
                    (   meijerg(
                            ((S(1)/2,), (0, 0)),
                            ((0, S(1)/2), (0,)),
                            polar_lift(s)**2),
                        True)
                ),
                And(
                    Abs(periodic_argument(polar_lift(s)**2, oo)) < pi,
                    cos(Abs(periodic_argument(polar_lift(s)**2, oo))/2)*sqrt(Abs(s**2)) - 1 > 0,
                    Ne(s**2, 1))
            ),
            (
                Integral(exp(-s*x)*cosh(x), (x, 0, oo)),
                True))
    assert integrate(exp(-s*x)*sinh(x), (x, 0, oo)) == \
        Piecewise(
            (   -1/(s + 1)/2 - 1/(-s + 1)/2,
                And(
                    Ne(1/s, 1),
                    Abs(periodic_argument(s, oo)) < pi/2,
                    Abs(periodic_argument(s, oo)) <= pi/2,
                    cos(Abs(periodic_argument(s, oo)))*Abs(s) - 1 > 0)),
            (   Integral(exp(-s*x)*sinh(x), (x, 0, oo)),
                True))

    def compute_psi_stats(self):
        # define some normal distributions
        mus = [sp.var("mu_%i" % i, real=True) for i in range(self.input_dim)]
        Ss = [sp.var("S_%i" % i, positive=True) for i in range(self.input_dim)]
        normals = [
            (2 * sp.pi * Si) ** (-0.5) * sp.exp(-0.5 * (xi - mui) ** 2 / Si) for xi, mui, Si in zip(self._sp_x, mus, Ss)
        ]

        # do some integration!
        # self._sp_psi0 = ??
        self._sp_psi1 = self._sp_k
        for i in range(self.input_dim):
            print "perfoming integrals %i of %i" % (i + 1, 2 * self.input_dim)
            sys.stdout.flush()
            self._sp_psi1 *= normals[i]
            self._sp_psi1 = sp.integrate(self._sp_psi1, (self._sp_x[i], -sp.oo, sp.oo))
            clear_cache()
        self._sp_psi1 = self._sp_psi1.simplify()

        # and here's psi2 (eek!)
        zprime = [sp.Symbol("zp%i" % i) for i in range(self.input_dim)]
        self._sp_psi2 = self._sp_k.copy() * self._sp_k.copy().subs(zip(self._sp_z, zprime))
        for i in range(self.input_dim):
            print "perfoming integrals %i of %i" % (self.input_dim + i + 1, 2 * self.input_dim)
            sys.stdout.flush()
            self._sp_psi2 *= normals[i]
            self._sp_psi2 = sp.integrate(self._sp_psi2, (self._sp_x[i], -sp.oo, sp.oo))
            clear_cache()
        self._sp_psi2 = self._sp_psi2.simplify()

 def E(expr):
     res1 = integrate(expr*exponential(x, rate)*normal(y, mu1, sigma1),
                      (x, 0, oo), (y, -oo, oo), meijerg=True)
     res2 = integrate(expr*exponential(x, rate)*normal(y, mu1, sigma1),
                     (y, -oo, oo), (x, 0, oo), meijerg=True)
     assert expand_mul(res1) == expand_mul(res2)
     return res1

def test_integrate_linearterm_pow():
    # check integrate((a*x+b)^c, x)  --  #400
    y = Symbol("y")
    assert integrate(x ** y, x) == x ** (y + 1) / (y + 1)
    assert integrate((exp(y) * x + 1 / y) ** (1 + sin(y)), x) == exp(-y) * (exp(y) * x + 1 / y) ** (2 + sin(y)) / (
        2 + sin(y)
    )

def test_issue_1791():
    z = Symbol("z", positive=True)
    assert integrate(exp(-log(x) ** 2), x) == pi ** (S(1) / 2) * erf(-S(1) / 2 + log(x)) * exp(S(1) / 4) / 2
    assert integrate(exp(log(x) ** 2), x) == -I * pi ** (S(1) / 2) * erf(I * log(x) + I / 2) * exp(-S(1) / 4) / 2
    assert integrate(exp(-z * log(x) ** 2), x) == pi ** (S(1) / 2) * erf(
        z ** (S(1) / 2) * log(x) - 1 / (2 * z ** (S(1) / 2))
    ) * exp(S(1) / (4 * z)) / (2 * z ** (S(1) / 2))

def T_exact_symbolic(verbose=False):
    """Compute the exact solution formula via sympy."""
    # sol1: solution for t < t_star,
    # sol2: solution for t > t_star
    import sympy as sym
    T0 = sym.symbols('T0')
    k = sym.symbols('k', positive=True)
    # Piecewise linear T_sunction
    t, t_star, C0, C1 = sym.symbols('t t_star C0 C1')
    T_s = C0
    I = sym.integrate(sym.exp(k*t)*T_s, (t, 0, t))
    sol1 = T0*sym.exp(-k*t) + k*sym.exp(-k*t)*I
    sol1 = sym.simplify(sym.expand(sol1))
    if verbose:
        # Some debugging print
        print 'solution t < t_star:', sol1
        #print sym.latex(sol1)
    T_s = C1
    I = sym.integrate(sym.exp(k*t)*C0, (t, 0, t_star)) + \
        sym.integrate(sym.exp(k*t)*C1, (t, t_star, t))
    sol2 = T0*sym.exp(-k*t) + k*sym.exp(-k*t)*I
    sol2 = sym.simplify(sym.expand(sol2))
    if verbose:
        print 'solution t > t_star:', sol2
        #print sym.latex(sol2)

    # Convert to numerical functions
    exact0 = sym.lambdify([t, C0, k, T0],
                          sol1, modules='numpy')
    exact1 = sym.lambdify([t, C0, C1, t_star, k, T0],
                          sol2, modules='numpy')
    return exact0, exact1

    def __init__(
            self,
            tau,  # A float
            start=None,
            ):
        print('''
tau={0}'''.format(tau))
        if start==None:
            start = tau -.05 + 2*tau
        x,lam = sympy.symbols('x lam'.split())
        self.tau = tau
        k = 0
        term = 1.0 + x*0   # \frac{1}{k!} * (-\frac{x}{lam})^k
        pk = 0             # sum_{n=0}^k term_k
        cum_prod = 1       # prod_{n=0}^k pk(tau)
        integral = 0       # integral_0^? rho(s) ds
        while tau*k < 1.0:
            pk += term              # Polynomial of x/lam
            pk_tau = pk.subs(x,tau) # Polynomial of tau/lam
            print('''  k={0}
  term={1}
  pk={2}
  '''.format(k, term, pk))
            if tau*(k+1) >= 1.0:
                d = 1.0 - tau*k
                integral += sympy.integrate(pk, (x, 0, d))/cum_prod
                break
            integral += sympy.integrate(pk, (x, 0, tau))/cum_prod
            # integral is a rational function of tau and lam
            cum_prod *= pk_tau      # Polynomial of tau/lam
            k += 1
            term *= -x/(k*lam)      # Monomial of x/lam
        print('''  cum_prod={2}  Solving
  {0}={1}'''.format(lam,integral.simplify(), cum_prod))
        self.eigenvalue = sympy.solve(lam-integral, lam)

    def integrar(self):

        f = self.ui.lineEdit_int_funcion.text()

        f_var = [Symbol(self.ui.lineEdit_int_var_1.text()),
            Symbol(self.ui.lineEdit_int_var_2.text()),
            Symbol(self.ui.lineEdit_int_var_3.text())]

        f_min = [(self.ui.lineEdit_int_min_1.text()),
            (self.ui.lineEdit_int_min_2.text()),
            (self.ui.lineEdit_int_min_3.text())]

        f_max = [(self.ui.lineEdit_int_max_1.text()),
            (self.ui.lineEdit_int_max_2.text()),
            (self.ui.lineEdit_int_max_3.text())]

        for i in range(0, int(self.ui.spinBox_int.text())):
            if self.ui.radioButton_int_def.isChecked() is True:
                f = integrate(f, (f_var[i], f_min[i], f_max[i]))
                self.ui.label_int_solucion.setText(str(f))
                print (f)

            else:
                f = integrate(f, f_var[i])
                self.ui.label_int_solucion.setText(str(f))
                print (f)

def benefit_from_demand(x, p, demand):
    """Converts the demand curve to the benefit.  It assumes that the
    demand is a x=d(p) function, where the x= is implicit.

    >>> sp.var('x p')
    (x, p)
    >>> sp.simplify(benefit_from_demand(x, p, 10/p -1) -
    ...             10*sp.log(x+1))
    0
    >>> sp.simplify(benefit_from_demand(x, p, sp.Eq(x, 10/p -1)) -
    ...             10*sp.log(x+1))
    0
    >>> sp.simplify(benefit_from_demand(x, p, 100-p) -
    ...             (-x**2/2 + 100*x))
    0
    >>> benefit_from_demand(x, p, sp.Piecewise((0, p < 0),
    ...                                        (-p + 100, p <= 100),
    ...                                        (0, True)))
    -x**2/2 + 100*x
    """
    if isinstance(demand, sp.relational.Relational):
        return sp.integrate(sp.solve(demand, p)[0], (x, 0, x))
    substracting = sp.solve(demand-x, p)
    if substracting:
        toint = substracting[0]
    else:
        substracting = sp.solve(demand, p)
        if substracting:
            toint = substracting[0] - x
        else:
            return None

    return sp.integrate(toint, (x, 0, x))

def integrate_f_over_polygon_code(f):
    """Generate code that will integrate a function over a polygon

    Parameters
    ----------
    f : sympy expression
        Expression to be integrated over the polygon.

    Returns
    -------
    out : str
        Multiline string of function code

    """

    x, y, z=sympy.symbols('x,y,z')
    x0, x1, y0, y1, z0, z1=sympy.symbols('x0,x1,y0,y1,z0,z1')
    t = sympy.symbols('t')
    s2 = [(x, x0+t*(x1-x0)), (y, y0+t*(y1-y0)), (z, z0+t*(z1-z0))]

    ff = sympy.integrate(f,x)
    ff = ff.subs(s2)
    ff = sympy.integrate(ff, (t,0,1))
    ff *= (y1 - y0) #if integrating in y direction 1st then *-(x1-x0)

    ff = replace_x0_and_x1_with_vect(ff)

    template ="""def ifxy(pts):
    "Integrate f = {} over polygon"

    x, y, z = xyz_from_pts(pts, True)

    return np.sum({})"""

    return template.format(str(f), ff)

def case0(f, N=3):
    B = 1 - x ** 3
    dBdx = sm.diff(B, x)

    # Compute basis functions and their derivatives
    phi = {0: [x ** (i + 1) * (1 - x) for i in range(N + 1)]}
    phi[1] = [sm.diff(phi_i, x) for phi_i in phi[0]]

    def integrand_lhs(phi, i, j):
        return phi[1][i] * phi[1][j]

    def integrand_rhs(phi, i):
        return f * phi[0][i] - dBdx * phi[1][i]

    Omega = [0, 1]

    u_bar = solve(integrand_lhs, integrand_rhs, phi, Omega, verbose=True, numint=False)
    u = B + u_bar
    print "solution u:", sm.simplify(sm.expand(u))

    # Calculate analytical solution

    # Solve -u''=f by integrating f twice
    f1 = sm.integrate(f, x)
    f2 = sm.integrate(f1, x)
    # Add integration constants
    C1, C2 = sm.symbols("C1 C2")
    u_e = -f2 + C1 * x + C2
    # Find C1 and C2 from the boundary conditions u(0)=0, u(1)=1
    s = sm.solve([u_e.subs(x, 0) - 1, u_e.subs(x, 1) - 0], [C1, C2])
    # Form the exact solution
    u_e = -f2 + s[C1] * x + s[C2]
    print "analytical solution:", u_e
    # print 'error:', u - u_e  # many terms - which cancel
    print "error:", sm.expand(u - u_e)

def test_issue_841():
    a, b, c, d = symbols('a:d', positive=True, bounded=True)
    assert integrate(
        exp(-x**2 + I*c*x), x) == sqrt(pi)*erf(x - I*c/2)*exp(-c**S(2)/4)/2
    assert integrate(exp(a*x**2 + b*x + c), x) == \
          I*sqrt(pi)*erf(-I*x*sqrt(
              a) - I*b/(2*sqrt(a)))*exp(c)*exp(-b**2/(4*a))/(2*sqrt(a))

def least_squares_orth(f, phi, Omega):
    """
    Same as least_squares, but for orthogonal
    basis such that one avoids calling up standard
    Gaussian elimination.
    """
    N = len(phi) - 1
    A = [0]*(N+1)       # plain list to hold symbolic expressions
    b = [0]*(N+1)
    x = sm.Symbol('x')
    print '...evaluating matrix...',
    for i in range(N+1):
        print '(%d,%d)' % (i, i)
        A[i] = sm.integrate(phi[i]**2, (x, Omega[0], Omega[1]))

        # Fallback on numerical integration if f*phi is too difficult
        # to integrate
        integrand = phi[i]*f
        I = sm.integrate(integrand,  (x, Omega[0], Omega[1]))
        if isinstance(I, sm.Integral):
            print 'numerical integration of', integrand
            integrand = sm.lambdify([x], integrand)
            I = sm.mpmath.quad(integrand, [Omega[0], Omega[1]])
        b[i] = I
    print 'A:\n', A, '\nb:\n', b
    c = [b[i]/A[i] for i in range(len(b))]
    print 'coeff:', c
    u = 0
    for i in range(len(phi)):
        u += c[i]*phi[i]
    # Alternative:
    # u = sum(c[i,0]*phi[i] for i in range(len(phi)))
    print 'approximation:', u
    return u

 def __init__(self, cu,cv):
   x, y, z = sympy.symbols('x y z')
   self.cu=cu
   self.cv=cv
   self.u=sympy.Lambda((x,y,z),cu)
   self.v=sympy.Lambda((x,y,z),cv)
   self.fu=eval("lambda x,y,z: "+cu)
   self.fv=eval("lambda x,y,z: "+cv)
   if self.u.is_polynomial(x,y,z):
     D=sympy.polys.Poly(cu,x,y,z).as_dict()
     self.du=max(numpy.sum(list(D.keys()),axis=1))
   else:
     self.du=-1
   if self.v.is_polynomial(x,y,z):
     D=sympy.polys.Poly(cv,x,y,z).as_dict()
     self.dv=max(numpy.sum(list(D.keys()),axis=1))
   else:
     self.dv=-1
   self.I=sympy.integrate(
            sympy.integrate(
              sympy.integrate(
                sympy.diff(self.u(x,y,z),x)*sympy.diff(self.v(x,y,z),x)+
                sympy.diff(self.u(x,y,z),y)*sympy.diff(self.v(x,y,z),y)+
                sympy.diff(self.u(x,y,z),z)*sympy.diff(self.v(x,y,z),z),
                (x,0,1)),
              (y,0,1)),
            (z,0,1))

def test_issue524():
    assert integrate(cos((n+1) * x), x)   == sin(x*(n+1)) / (n+1)
    assert integrate(cos((n-1) * x), x)   == sin(x*(n-1)) / (n-1)

    assert integrate(cos((n+1) * x) + cos((n-1) * x), x) == \
                                             sin(x*(n+1)) / (n+1)  + \
                                             sin(x*(n-1)) / (n-1)

def FEM_function_norms(N):

	h = 1.0/N
	print "N = ", N, " |  h = ", h

	# set up mesh points
	x = linspace(0,1,N+1)

	# set up function u(x)
	u = zeros(N+1)
	u[int(N/2)] = 1.0

	print "**** Analytical norms ****"
	import sympy as sym
	X = sym.Symbol('X')
	u1 = sym.Rational(1,h)*X - sym.Rational(1,h)*(0.5-h) 
	u2 = sym.Rational(-1,h)*X + sym.Rational(1,h)*(0.5+h) 

	L2 = sym.integrate(u1*u1,(X,0.5-h,0.5)) + sym.integrate(u2*u2,(X,0.5,0.5+h))
	print "L2: ", sqrt(float(L2))

	H1 = L2 + sym.integrate(sym.diff(u1,X)*sym.diff(u1,X),(X,0.5-h,0.5)) + sym.integrate(sym.diff(u1,X)*sym.diff(u1,X),(X,0.5,0.5+h))
	print "H1: ", sqrt(float(H1))

	print "**** Numerical norms ****"

	num_L2 = sum((u**2)/len(u))
	print  "L2: ", sqrt(num_L2) 

	dudx = (u[1:] - u[:-1])/(x[1]-x[0]) 
	num_H1 = sum(u[1:]**2 + dudx**2 )/len(u)
	print "H1: ", sqrt(num_H1) 

	return L2, H1