summaryrefslogtreecommitdiff
path: root/macros/luatex/latex/luacas/tex/calculus/integralexpression.lua
diff options
context:
space:
mode:
Diffstat (limited to 'macros/luatex/latex/luacas/tex/calculus/integralexpression.lua')
-rw-r--r--macros/luatex/latex/luacas/tex/calculus/integralexpression.lua934
1 files changed, 934 insertions, 0 deletions
diff --git a/macros/luatex/latex/luacas/tex/calculus/integralexpression.lua b/macros/luatex/latex/luacas/tex/calculus/integralexpression.lua
new file mode 100644
index 0000000000..5267127264
--- /dev/null
+++ b/macros/luatex/latex/luacas/tex/calculus/integralexpression.lua
@@ -0,0 +1,934 @@
+--- @class IntegralExpression
+--- An expression for the integral of an expression.
+--- @field symbol SymbolExpression
+--- @field expression Expression
+--- @field upper Expression
+--- @field lower Expression
+--- @field attempts table<number, Expression>
+--- @field results table<number, Expression>
+--- @field enhancedsubstitution Integer
+--- @field recursive boolean
+
+IntegralExpression = {}
+__IntegralExpression = {}
+
+
+--------------------------
+-- Static functionality --
+--------------------------
+
+--- Recursive part of the indefinite integral operator. Returns nil if the expression could not be integrated.
+--- We switch to prodcedural programming here because it is more natural.
+--- @param integral IntegralExpression
+--- @return Expression|nil
+function IntegralExpression.integrate(integral)
+ integral.expression = integral.expression:autosimplify()
+
+ if not integral.recursive and #integral.attempts > 0 then
+ return Copy(integral:lock(Expression.NIL, true))
+ end
+
+ -- print(integral.expression)
+
+ local F = IntegralExpression.table(integral)
+ if F then return F end
+
+ -- If we see the same integrand again, and hasn't been solved already, then the integral can't be solved
+ local resultindex = Contains(integral.attempts, integral.expression)
+ if resultindex then
+ return integral.results[resultindex]
+ end
+ local newindex = #integral.attempts+1
+ integral.attempts[newindex] = integral.expression
+
+ -- print("Evalutaing: " .. tostring(integral.expression))
+
+ F = IntegralExpression.linearproperties(integral)
+ if F then
+ -- print("Linear Properties")
+ integral.results[newindex] = F
+ return F
+ end
+
+ -- local exp = integral.expression
+ -- local sym = integral.symbol
+ -- local es = integral.enhancedsubstitution
+ -- integral.enhancedsubstitution = Integer.zero()
+ F = IntegralExpression.substitutionmethod(integral)
+ if F then
+ -- print("u-Substitution")
+ integral.results[newindex] = F
+ return F
+ end
+ -- integral.expression = exp
+ -- integral.symbol = sym
+ -- integral.enhancedsubstitution = es
+
+
+ F = IntegralExpression.rationalfunction(integral)
+ if F then
+ -- print("Rational Function")
+ integral.results[newindex] = F
+ return F
+ end
+
+ F = IntegralExpression.partsmethod(integral)
+ if F then
+ -- print("Parts")
+ integral.results[newindex] = F
+ return F
+ end
+
+ F = IntegralExpression.eulersformula(integral)
+ if F then
+ -- print("Euler's formula")
+ integral.results[newindex] = F
+ return F
+ end
+
+ local expanded = integral.expression:expand()
+ if integral.expression ~= expanded then
+ integral.expression = expanded
+ F = IntegralExpression.integrate(integral)
+ if F then
+ -- print("Expanded")
+ integral.results[newindex] = F
+ return F
+ end
+ end
+
+ expanded = (Integer.one()/((Integer.one()/integral.expression):autosimplify():expand())):autosimplify()
+ if integral.expression ~= expanded then
+ integral.expression = expanded
+ F = IntegralExpression.integrate(integral)
+ if F then
+ -- print("Inverse Expanded")
+ integral.results[newindex] = F
+ return F
+ end
+ end
+
+ F = IntegralExpression.enhancedsubstitutionmethod(integral)
+ if F then
+ -- print("Enhanced u-Substitution")
+ integral.results[newindex] = F
+ return F
+ end
+
+ return nil
+end
+
+--- A table of basic integrals, returns nil if the integrand isn't in the table.
+--- @param integral IntegralExpression
+--- @return Expression|nil
+function IntegralExpression.table(integral)
+ local integrand = integral.expression
+ local symbol = integral.symbol
+
+ -- Constant integrand rule - int(c, x) = c*x
+ if integrand:freeof(symbol) then
+ return integrand*symbol
+ end
+
+ if integrand:type() == SymbolExpression then
+
+ -- int(x, x) = x^2/2
+ if integrand == symbol then
+ return integrand ^ Integer(2) / Integer(2)
+ end
+
+ -- Constant integrand rule again
+ return integrand*symbol
+ end
+
+ if integrand:type() == BinaryOperation then
+
+ if integrand.operation == BinaryOperation.POW then
+ -- int(1/x, x) = ln(x)
+ if integrand.expressions[1] == symbol and integrand.expressions[2] == Integer(-1) then
+ return LN(symbol)
+ end
+
+ -- Cavalieri's formula - int(x^n, x) = x^(n+1)/(n+1)
+ if integrand.expressions[1] == symbol and integrand.expressions[2]:freeof(symbol) then
+ return symbol ^ (integrand.expressions[2] + Integer.one()) / (integrand.expressions[2] + Integer.one())
+ end
+
+ -- int(n^x, x) = n^x/ln(n)
+ if integrand.expressions[1]:freeof(symbol) and integrand.expressions[2] == symbol then
+ return integrand / LN(integrand.expressions[1])
+ end
+
+ -- int(csc(x)^2, x) = -cot(x)
+ if integrand.expressions[1] == CSC(symbol) and integrand.expressions[2] == Integer(2) then
+ return -COT(symbol)
+ end
+
+ -- int(sec(x)^2, x) = tan(x)
+ if integrand.expressions[1] == SEC(symbol) and integrand.expressions[2] == Integer(2) then
+ return TAN(symbol)
+ end
+ end
+
+ if integrand.operation == BinaryOperation.MUL and #integrand.expressions == 2 then
+ -- int(tan(x)sec(x), x) = sec(x)
+ if integrand.expressions[1] == TAN(symbol) and integrand.expressions[2] == SEC(symbol) then
+ return SEC(symbol)
+ end
+
+ -- int(csc(x)cot(x), x) = -csc(x)
+ if integrand.expressions[1] == CSC(symbol) and integrand.expressions[2] == COT(symbol) then
+ return -CSC(symbol)
+ end
+ end
+
+ return nil
+ end
+
+ if integrand:type() == Logarithm then
+ -- int(log_n(x), x) = (x*ln(x)-x)/ln(n)
+ if integrand.base:freeof(symbol) and integrand.expression == symbol then
+ return (symbol * LN(symbol) - symbol) / LN(integrand.base)
+ end
+
+ return nil
+ end
+
+ if integrand:type() == TrigExpression then
+ if integrand == SIN(symbol) then
+ return -COS(symbol)
+ end
+
+ if integrand == COS(symbol) then
+ return SIN(symbol)
+ end
+
+ if integrand == TAN(symbol) then
+ return -LN(COS(symbol))
+ end
+
+ if integrand == CSC(symbol) then
+ return -LN(CSC(symbol)+COT(symbol))
+ end
+
+ if integrand == SEC(symbol) then
+ return LN(SEC(symbol) + TAN(symbol))
+ end
+
+ if integrand == COT(symbol) then
+ return LN(SIN(symbol))
+ end
+
+ if integrand == ARCSIN(symbol) then
+ return symbol*ARCSIN(symbol) + (Integer.one()-symbol^(Integer(2)))^(Integer.one()/Integer(2))
+ end
+
+ if integrand == ARCCOS(symbol) then
+ return symbol*ARCCOS(symbol) - (Integer.one()-symbol^(Integer(2)))^(Integer.one()/Integer(2))
+ end
+
+ if integrand == ARCTAN(symbol) then
+ return symbol*ARCTAN(symbol) - (Integer.one()/Integer(2))*LN(Integer.one()+symbol^Integer(2))
+ end
+
+ if integrand == ARCCSC(symbol) then
+ return symbol*ARCCSC(symbol) + LN(symbol*(Integer.one()+(Integer.one()-symbol^(Integer(-2)))^(Integer.one()/Integer(2))))
+ end
+
+ if integrand == ARCSEC(symbol) then
+ return symbol*ARCSEC(symbol) - LN(symbol*(Integer.one()+(Integer.one()-symbol^(Integer(-2)))^(Integer.one()/Integer(2))))
+ end
+
+ if integrand == ARCCOT(symbol) then
+ return symbol*ARCCOT(symbol) + (Integer.one()/Integer(2))*LN(Integer.one()+symbol^Integer(2))
+ end
+ end
+
+ return nil
+end
+
+--- Uses linearity to break up the integral and integrate each piece.
+--- @param integral IntegralExpression
+--- @return Expression|nil
+function IntegralExpression.linearproperties(integral)
+ local expression = integral.expression
+ local symbol = integral.symbol
+ local es = integral.enhancedsubstitution
+
+ if expression:type() == BinaryOperation then
+ if expression.operation == BinaryOperation.MUL then
+ local freepart = Integer.one()
+ local variablepart = Integer.one()
+ for _, term in ipairs(expression.expressions) do
+ if term:freeof(symbol) then
+ freepart = freepart*term
+ else
+ variablepart = variablepart*term
+ end
+ end
+ if freepart == Integer.one() then
+ return nil
+ end
+ integral.expression = variablepart
+ integral.symbol = symbol
+ integral.enhancedsubstitution = es
+ local F = IntegralExpression.integrate(integral)
+ if F then
+ return freepart*F
+ end
+ return nil
+ end
+
+ if expression.operation == BinaryOperation.ADD then
+ local sum = Integer.zero()
+ for _, term in ipairs(expression.expressions) do
+ integral.expression = term
+ integral.symbol = symbol
+ integral.enhancedsubstitution = es
+ local F = IntegralExpression.integrate(integral)
+ if F then
+ sum = sum + F
+ else
+ return nil
+ end
+
+ end
+ return sum
+ end
+ end
+
+ return nil
+end
+
+--- Attempts u-substitutions to evaluate the integral.
+--- @param integral IntegralExpression
+--- @return Expression|nil
+function IntegralExpression.substitutionmethod(integral)
+ local expression = integral.expression
+ local symbol = integral.symbol
+ local es = integral.enhancedsubstitution
+
+ local P = IntegralExpression.trialsubstitutions(expression)
+ local F = nil
+ local i = 1
+
+ while not F and i <= #P do
+ local g = P[i]
+ if g ~= symbol and not g:freeof(symbol) then
+ local subsymbol = SymbolExpression("u")
+ if symbol == SymbolExpression("u") then
+ subsymbol = SymbolExpression("v")
+ end
+ local u = (expression / (DerivativeExpression(g, symbol))):autosimplify()
+ u = u:substitute({[g]=subsymbol}):autosimplify()
+
+ --factor u and cancel like non-constant terms
+ u = u:factor():autosimplify()
+
+ if u:freeof(symbol) then
+ integral.expression = u
+ integral.symbol = subsymbol
+ integral.enhancedsubstitution = es
+ F = IntegralExpression.integrate(integral)
+ if F then
+ if integral.recursive then
+ F = F:substitute({[subsymbol]=g})
+ end
+ return F
+ end
+ end
+ end
+ i = i + 1
+ end
+
+ return F
+end
+
+--- Attempts u-substitutions to evaluate the integral, including solving for the original variable and substituting the result into the expression.
+--- @param integral IntegralExpression
+--- @return Expression|nil
+function IntegralExpression.enhancedsubstitutionmethod(integral)
+ local expression = integral.expression
+ local symbol = integral.symbol
+ local es = integral.enhancedsubstitution
+
+ local P = IntegralExpression.trialsubstitutions(expression)
+ local F = nil
+ local i = 1
+
+ while not F and i <= #P do
+ local g = P[i]
+ if g ~= symbol and not g:freeof(symbol) then
+ local subsymbol = SymbolExpression("u")
+ if symbol == SymbolExpression("u") then
+ subsymbol = SymbolExpression("v")
+ end
+ local u = (expression / (DerivativeExpression(g, symbol))):autosimplify()
+ u = u:substitute({[g]=subsymbol}):autosimplify()
+
+ --factor u and cancel like non-constant terms
+ u = u:factor():autosimplify()
+
+ if integral.enhancedsubstitution > Integer.zero() then
+ local f = Equation(subsymbol, g):solvefor(symbol)
+ if f.lhs == symbol then
+ u = u:substitute({[symbol]=f.rhs}):autosimplify()
+ integral.expression = u
+ integral.symbol = subsymbol
+ integral.enhancedsubstitution = integral.enhancedsubstitution - Integer.one()
+ F = IntegralExpression.integrate(integral)
+ if F then
+ if integral.recursive then
+ F = F:substitute({[subsymbol]=g})
+ end
+ return F
+ end
+ integral.enhancedsubstitution = integral.enhancedsubstitution + Integer.one()
+ end
+ end
+ end
+ i = i + 1
+ end
+
+ return F
+end
+
+--- Generates a list of possible u-substitutions to attempt
+--- @param expression Expression
+--- @return table<number, Expression>
+function IntegralExpression.trialsubstitutions(expression)
+ local substitutions = {}
+
+ -- Recursive part - evaluates each term in a product.
+ if expression:type() == BinaryOperation and expression.operation == BinaryOperation.MUL then
+ substitutions[#substitutions+1] = expression
+ for _, term in ipairs(expression.expressions) do
+ substitutions = JoinArrays(substitutions, IntegralExpression.trialsubstitutions(term))
+ end
+ end
+
+ --Recursive part - evaluates each term in a sum.
+ if expression:type() == BinaryOperation and expression.operation == BinaryOperation.ADD then
+ substitutions[#substitutions+1] = expression
+ for _,term in ipairs(expression.expressions) do
+ substitutions = JoinArrays(substitutions, IntegralExpression.trialsubstitutions(term))
+ end
+ end
+
+ -- Function forms and arguments of function forms (includes a recursive part)
+ if expression:type() == TrigExpression or expression:type() == Logarithm then
+ substitutions[#substitutions+1] = expression
+ if not expression.expression:isatomic() then
+ substitutions[#substitutions+1] = expression.expression
+ end
+ substitutions = JoinArrays(substitutions, IntegralExpression.trialsubstitutions(expression.expression))
+ end
+
+ -- Bases and exponents of powers
+ if expression:type() == BinaryOperation and expression.operation == BinaryOperation.POW then
+ substitutions[#substitutions+1] = expression
+ -- Atomic expressions are technically valid substitutions, but they won't be useful
+ if not expression.expressions[1]:isatomic() then
+ --substitutions[#substitutions+1] = expression.expressions[1]
+ substitutions = JoinArrays(substitutions, IntegralExpression.trialsubstitutions(expression.expressions[1]))
+ end
+ if not expression.expressions[2]:isatomic() then
+ --substitutions[#substitutions+1] = expression.expressions[2]
+ substitutions = JoinArrays(substitutions, IntegralExpression.trialsubstitutions(expression.expressions[2]))
+ end
+ end
+
+ return substitutions
+end
+
+
+--- Uses Lazard, Rioboo, Rothstein, and Trager's method to integrate rational functions.
+--- This is mostly to try to avoid factoring and finding the roots of the full denominator whenever possible.
+--- @param integral IntegralExpression
+--- @return Expression|nil
+function IntegralExpression.rationalfunction(integral)
+ local expression = integral.expression
+ local symbol = integral.symbol
+ local es = integral.enhancedsubstitution
+
+ -- Type checking and conversion to polynomial type.
+ local f, g, fstat, gstat
+ if expression:type() == BinaryOperation and expression.operation == BinaryOperation.POW and expression.expressions[2] == Integer(-1) then
+ g, gstat = expression.expressions[1]:topolynomial()
+ if not gstat then
+ return nil
+ end
+ f = PolynomialRing({Integer.one()}, g.symbol)
+ else
+ if expression:type() ~= BinaryOperation or expression.operation ~= BinaryOperation.MUL or expression.expressions[3] then
+ return nil
+ end
+ if expression.expressions[2]:type() == BinaryOperation and expression.expressions[2].operation == BinaryOperation.POW and expression.expressions[2].expressions[2] == Integer(-1) then
+ if expression.expressions[1].topolynomial ~=nil and expression.expressions[2].expressions[1].topolynomial ~=nil then
+ f, fstat = expression.expressions[1]:topolynomial()
+ g, gstat = expression.expressions[2].expressions[1]:topolynomial()
+ end
+ elseif expression.expressions[1]:type() == BinaryOperation and expression.expressions[1].operation == BinaryOperation.POW and expression.expressions[1].expressions[2] == Integer(-1) then
+ if expression.expressions[2].topolynomial ~= nil and expression.expressions[1].expressions[1].topolynomial ~= nil then
+ f, fstat = expression.expressions[2]:topolynomial()
+ g, gstat = expression.expressions[1].expressions[1]:topolynomial()
+ end
+ else
+ return nil
+ end
+
+ if not fstat or not gstat or f.symbol ~= symbol.symbol or g.symbol ~= symbol.symbol then
+ return nil
+ end
+ end
+
+ -- Explicit handling of degree 1 or less over a binomial.
+ do
+ local disc = g.coefficients[1]*g.coefficients[1]-Integer(4)*g.coefficients[2]*g.coefficients[0]
+ if f.degree <= Integer.one() and g.degree == Integer(2) and disc < Integer.zero() then
+ return (f.coefficients[1] * LN(g.coefficients[0] + g.coefficients[1] * symbol + g.coefficients[2] * symbol ^ Integer(2))/(Integer(2) * g.coefficients[2]) + (Integer(2)*f.coefficients[0]*g.coefficients[2] - f.coefficients[1]*g.coefficients[1]) / (g.coefficients[2] * sqrt(Integer(4)*g.coefficients[0]*g.coefficients[2] - g.coefficients[1] ^ Integer(2))) * ARCTAN((Integer(2)*g.coefficients[2]*symbol+g.coefficients[1]) / sqrt(Integer(4)*g.coefficients[0]*g.coefficients[2]-g.coefficients[1] ^ Integer(2)))):autosimplify()
+ end
+ end
+
+ -- If the polynomials are not relatively prime, divides out the common factors.
+ local gcd = PolynomialRing.gcd(f, g)
+ if gcd ~= Integer.one() then
+ f, g = f // gcd, g // gcd
+ end
+
+ -- Seperates out the polynomial part and rational part and integrates the polynomial part.
+ local q, h = f:divremainder(g)
+ integral.expression = q
+ integral.symbol = symbol
+ integral.enhancedsubstitution = es
+ U = IntegralExpression.integrate(integral)
+
+ if h == Integer.zero() then
+ return U
+ end
+
+ -- Performs partial fraction decomposition into square-free denominators on the rational part.
+ local gg = g:squarefreefactorization()
+ local pfd = PolynomialRing.partialfractions(h, g, gg)
+
+ -- Hermite reduction.
+ local V = Integer.zero()
+ for _, term in ipairs(pfd.expressions) do
+ local i = #term.expressions
+ if i > 1 then
+ for j = 1, i-1 do
+ local n = term.expressions[j].expressions[1]
+ local d = term.expressions[j].expressions[2].expressions[1]
+ local p = term.expressions[j].expressions[2].expressions[2]
+
+ local _, s, t = PolynomialRing.extendedgcd(d, d:derivative())
+ s = s * n
+ t = t * n
+ V = V - t / ((p-Integer.one()) * BinaryOperation.POWEXP({d, p-Integer.one()}))
+ term.expressions[j+1].expressions[1] = term.expressions[j+1].expressions[1] + s + t:derivative() / (p-Integer.one())
+ end
+ end
+ end
+
+ --Lazard-Rioboo-Trager method.
+ local W = Integer.zero()
+ for _, term in ipairs(pfd.expressions) do
+ local a = term.expressions[#term.expressions].expressions[1]
+ local b = term.expressions[1].expressions[2].expressions[1]
+ local y = a - b:derivative() * PolynomialRing({Integer.zero(), Integer.one()}, "_")
+ local r = PolynomialRing.resultant(b, y)
+
+
+ local rr = r:squarefreefactorization()
+ local remainders = PolynomialRing.monicgcdremainders(b, y)
+ for pos, factor in ipairs(rr.expressions) do
+ if pos > 1 then
+ local re = factor.expressions[1]
+ local e = factor.expressions[2]
+ local roots = re:roots()
+ for _, root in ipairs(roots) do
+ local w
+ for _, remainder in ipairs(remainders) do
+ if remainder.degree == e then
+ w = remainder
+ break
+ end
+ end
+ W = W + root*LN(w:substitute({[SymbolExpression("_")] = root}))
+ end
+ end
+ end
+ end
+
+ return U + V + W
+end
+
+
+--- Attempts integration by parts for expressions with a polynomial factor in them. Other product expressions use Euler's formula.
+--- @param integral IntegralExpression
+--- @return Expression|nil
+function IntegralExpression.partsmethod(integral)
+ local expression = integral.expression
+ local symbol = integral.symbol
+ local es = integral.enhancedsubstitution
+
+ if expression:type() ~= BinaryOperation or expression.operation ~= BinaryOperation.MUL then
+ return
+ end
+
+ local u
+ local vp = Integer.one()
+ --looking for ILATE
+ for _, exp in ipairs(expression:subexpressions()) do
+ if exp:type() == TrigExpression and (exp.name == "arctan" or exp.name == "arccos" or exp.name == "arcsin" or exp.name == "arccot" or exp.name == "arcsec" or exp.name == "arccsc") then
+ u = exp
+ else
+ vp = vp * exp
+ end
+ end
+
+ if not u or u:freeof(symbol) then
+ goto skipI
+ else
+ vp = vp:autosimplify()
+ end
+
+ --if vp:type() == Logarithm or vp.topolynomial or (vp:type() == TrigExpression and (vp.name == "cos" or vp.name == "sin")) or (vp.operation == BinaryOperation.POW and vp.expressions[1]:freeof(symbol)) then
+ if select(2,vp:topolynomial()) then
+ integral.expression = vp
+ integral.symbol = symbol
+ integral.enhancedsubstitution = es
+ local v = IntegralExpression.integrate(integral)
+ if not v then
+ goto skipI
+ end
+
+ local up = DerivativeExpression(u, symbol):autosimplify()
+
+ integral.expression = v*up
+ integral.symbol = symbol
+ integral.enhancedsubstitution = es
+ local vup = IntegralExpression.integrate(integral)
+ if not vup then
+ goto skipI
+ end
+
+ local result = u*v - vup
+
+ return result:autosimplify()
+ end
+ ::skipI::
+
+ local u
+ local vp = Integer.one()
+ --looking for LATE
+ for _, exp in ipairs(expression:subexpressions()) do
+ if exp:type() == Logarithm then
+ u = exp
+ else
+ vp = vp * exp
+ end
+ end
+
+ if not u or u:freeof(symbol) then
+ goto skipL
+ else
+ vp = vp:autosimplify()
+ end
+
+ --if vp.topolynomial or (vp:type() == TrigExpression and (vp.name == "cos" or vp.name == "sin")) or (vp.operation == BinaryOperation.POW and vp.expressions[1]:freeof(symbol)) then
+ if select(2,vp:topolynomial()) then
+ integral.expression = vp
+ integral.symbol = symbol
+ integral.enhancedsubstitution = es
+ local v = IntegralExpression.integrate(integral)
+ if not v then
+ goto skipL
+ end
+
+ local up = DerivativeExpression(u, symbol):autosimplify()
+
+ integral.expression = v*up
+ integral.symbol = symbol
+ integral.enhancedsubstitution = es
+ local vup = IntegralExpression.integrate(integral)
+ if not vup then
+ goto skipL
+ end
+
+ local result = u*v - vup
+
+ return result:autosimplify()
+ end
+ ::skipL::
+
+ local u
+ local vp = Integer.one()
+ --looking for ATE
+ for _, exp in ipairs(expression:subexpressions()) do
+ local _, bool = exp:topolynomial()
+ if bool then
+ u = exp
+ else
+ vp = vp * exp
+ end
+ end
+
+ if not u or u:freeof(symbol) then
+ return
+ else
+ vp = vp:autosimplify()
+ end
+
+ if (vp:type() == TrigExpression and (vp.name == "cos" or vp.name == "sin")) or (vp.operation == BinaryOperation.POW and vp.expressions[1]:freeof(symbol)) then
+ local results = {}
+ while u ~= Integer.zero() do
+ integral.expression = vp
+ integral.symbol = symbol
+ integral.enhancedsubstitution = es
+ local v = IntegralExpression.integrate(integral):unlock():autosimplify()
+ if not v then
+ return
+ end
+ local up = DerivativeExpression(u, symbol):autosimplify()
+
+ if not integral.recursive then
+ return (u*v - IntegralExpression(v*up, symbol):nonrecursive():lock(Expression.NIL, true)):autosimplify()
+ end
+
+ results[#results+1] = u*v
+ u = up
+ vp = v
+ end
+
+ local result = results[#results]
+ for i=#results-1,1,-1 do
+ result = results[i] - result
+ end
+
+ return result:autosimplify()
+ end
+end
+
+--- Attempts integration using Euler's formula and kind. Alternative for integration by parts for many expressions.
+--- @param integral IntegralExpression
+--- @return Expression|nil
+function IntegralExpression.eulersformula(integral)
+ local expression = integral.expression
+ local symbol = integral.symbol
+ local es = integral.enhancedsubstitution
+
+ local new = expression:substitute({[COS(symbol)] = (E^(I*symbol) + E^(-I*symbol))/Integer(2),
+ [SIN(symbol)] = (E^(I*symbol) - E^(-I*symbol))/(Integer(2)*I)})
+
+ if new == expression then
+ return
+ end
+
+ integral.expression = new:autosimplify():expand()
+ integral.symbol = symbol
+ integral.enhancedsubstitution = es
+ local complexresult = IntegralExpression.integrate(integral)
+ if not complexresult then
+ return
+ end
+
+ -- TODO: Proper complex number conversion methods
+ local function converttorectangular(exp)
+ exp = exp:expand()
+ local results = {}
+ for index, sub in ipairs(exp:subexpressions()) do
+ results[index] = converttorectangular(sub)
+ end
+ local converted = exp:setsubexpressions(results)
+
+ if converted.operation == BinaryOperation.POW and converted.expressions[1] == E and converted.expressions[2].operation == BinaryOperation.MUL then
+ local ipart
+ local rest = Integer.one()
+ for _, factor in ipairs(converted.expressions[2]:subexpressions()) do
+ if factor == I then
+ ipart = true
+ else
+ rest = rest * factor
+ end
+ end
+ if ipart then
+ return (COS(rest) + I*SIN(rest)):autosimplify()
+ end
+ end
+
+ return converted
+ end
+
+ return converttorectangular(complexresult:autosimplify()):expand():autosimplify()
+
+end
+
+----------------------------
+-- Instance functionality --
+----------------------------
+
+--- Creates a new integral operation with the given symbol and expression.
+--- @param expression Expression
+--- @param symbol SymbolExpression
+--- @param lower Expression
+--- @param upper Expression
+function IntegralExpression:new(expression, symbol, lower, upper)
+ local o = {}
+ local __o = Copy(__ExpressionOperations)
+
+ if not symbol or not expression then
+ error("Send wrong number of parameters: integrals must have a variable to integrate with respect to and an expression to integrate.")
+ end
+
+ if lower and not upper then
+ error("Send wrong number of parameters: definite integrals must have an upper and a lower bound.")
+ end
+
+ o.symbol = symbol
+ o.expression = Copy(expression)
+ o.upper = Copy(upper)
+ o.lower = Copy(lower)
+ o.recursive = true
+
+ o.attempts = {}
+ o.results = {}
+ o.enhancedsubstitution = IntegralExpression.ENHANCEDSUBSTITUTIONRECURSIONLIMIT
+
+ __o.__index = IntegralExpression
+ __o.__tostring = function(a)
+ if a:isdefinite() then
+ return 'int(' .. tostring(a.expression) .. ", " .. tostring(a.symbol) .. ", ".. tostring(a.lower) .. ', ' .. tostring(a.upper) .. ')'
+ end
+ return 'int(' .. tostring(a.expression) .. ", " .. tostring(a.symbol) .. ')'
+ end
+ __o.__eq = function(a, b)
+ -- This shouldn't be needed, since __eq should only fire if both metamethods have the same function, but for some reason Lua always rungs this anyway
+ if not b:type() == IntegralExpression then
+ return false
+ end
+ return a.symbol == b.symbol and a.expression == b.expression and a.upper == b.upper and a.lower == b.lower
+ end
+ o = setmetatable(o, __o)
+
+ return o
+end
+
+--- Returns true if the integral is definite, i.e., has an upper and lower bound.
+--- @return boolean
+function IntegralExpression:isdefinite()
+ return self.upper ~= nil
+end
+
+--- Sets the integral to not autosimplify other integral expressions that are produced by the integration process.
+--- THIS METHOD MUTATES THE OBJECT IT IS CALLED ON.
+function IntegralExpression:nonrecursive()
+ self.recursive = false
+ return self
+end
+
+--- @return Expression
+function IntegralExpression:autosimplify()
+ local arg = IntegralExpression(self.expression, self.symbol)
+ local integrated = IntegralExpression.integrate(arg)
+
+ -- Our expression could not be integrated.
+ if not integrated then
+ return self
+ end
+
+ if not self.recursive then
+ return integrated:autosimplify():unlock(true)
+ end
+
+ if self:isdefinite() then
+ return (integrated:substitute({[self.symbol]=self.upper}) - integrated:substitute({[self.symbol]=self.lower})):autosimplify()
+ end
+
+ return integrated:autosimplify()
+end
+
+
+--- @return table<number, Expression>
+function IntegralExpression:subexpressions()
+ if self:isdefinite() then
+ return {self.expression, self.symbol, self.lower, self.upper}
+ end
+
+ return {self.expression, self.symbol}
+end
+
+--- @param subexpressions table<number, Expression>
+--- @return IntegralExpression
+function IntegralExpression:setsubexpressions(subexpressions)
+ local out = IntegralExpression(subexpressions[1], subexpressions[2], subexpressions[3], subexpressions[4])
+
+ return out;
+end
+
+-- function IntegralExpression:freeof(symbol)
+-- if self:isdefinite() then
+-- return self.expression:freeof(symbol) and self.upper:freeof(symbol) and self.lower:freeof(symbol)
+-- end
+-- return self.expression:freeof(symbol)
+-- end
+
+-- -- Substitutes each expression for a new one.
+-- function IntegralExpression:substitute(map)
+-- for expression, replacement in pairs(map) do
+-- if self == expression then
+-- return replacement
+-- end
+-- end
+-- -- Typically, we only perform substitution on autosimplified expressions, so this won't get called. May give strange results, i.e.,
+-- -- substituting and then evaluating the integral may not return the same thing as evaluating the integral and then substituting.
+-- if self:isdefinite() then
+-- return IntegralExpression(self.symbol, self.expression:substitute(map), self.upper:substitute(map), self.lower:substitute(map))
+-- end
+-- return IntegralExpression(self.symbol, self.expression:substitute(map))
+-- end
+
+--- @param other Expression
+--- @return boolean
+function IntegralExpression:order(other)
+ if other:type() ~= IntegralExpression then
+ return false
+ end
+
+ if self.symbol ~= other.symbol then
+ return self.symbol:order(other.symbol)
+ end
+
+ return self.expression:order(other.expression)
+end
+
+--- @return string
+function IntegralExpression:tolatex()
+ if self:isdefinite() then
+ return '\\int_{' .. self.lower:tolatex() .. '}^{' .. self.upper:tolatex() .. '}{' .. self.expression:tolatex() .. '\\mathop{d' .. self.symbol:tolatex() .. '}}'
+ end
+ return '\\int{' .. self.expression:tolatex() .. '\\mathop{d' .. self.symbol:tolatex() .. '}}'
+end
+
+
+-----------------
+-- Inheritance --
+-----------------
+__IntegralExpression.__index = CompoundExpression
+__IntegralExpression.__call = IntegralExpression.new
+IntegralExpression = setmetatable(IntegralExpression, __IntegralExpression)
+
+----------------------
+-- Static constants --
+----------------------
+INT = function(symbol, expression, lower, upper)
+ return IntegralExpression(symbol, expression, lower, upper)
+end
+
+----------------------
+-- Static constants --
+----------------------
+
+-- Limit for the maximum number of full u-subs to attempts for any integral.
+-- This should be low, since integrals are highly unlikely to need more than 1 or 2 u-subs, and gives exponentially worse performance the higher the number is.
+IntegralExpression.ENHANCEDSUBSTITUTIONRECURSIONLIMIT = Integer(2) \ No newline at end of file