Faster rules

Project.toml

@@ -27,8 +27,8 @@ FresnelIntegrals = "0.2.0"
 HypergeometricFunctions = "0.3.28"
 Nemo = "0.51, 0.52"
 PolyLog = "2.6.0"
-SymbolicUtils = "3"
-Symbolics = "6"
+SymbolicUtils = "4"
+Symbolics = "7"
 julia = "1.10"
 
 [extras]

src/SymbolicIntegration.jl

@@ -22,6 +22,7 @@ include("methods/rule_based/general.jl")
 include("methods/rule_based/frontend.jl")
 include("methods/rule_based/rules_utility_functions.jl")
 include("methods/rule_based/rules_loader.jl")
+include("methods/rule_based/rule2.jl")
 
 # Add method dispatch system
 include("methods.jl")

src/methods.jl

@@ -35,7 +35,7 @@ struct RuleBasedMethod <: AbstractIntegrationMethod
     use_gamma::Bool
     verbose::Bool
 
-    function RuleBasedMethod(; use_gamma::Bool=false, verbose::Bool=true)
+    function RuleBasedMethod(; use_gamma::Bool=false, verbose::Bool=false)
         new(use_gamma, verbose)
     end
 end
@@ -87,17 +87,19 @@ No rule found for ∫(abs(x), x)
 ```
 """
 function integrate(f::Symbolics.Num, x::Symbolics.Num; kwargs...)
+    result = integrate_rule_based(f, x; kwargs...)
+    !contains_int(result) && return result
+
+    printstyled(" > RuleBasedMethod(use_gamma=false, verbose=false) failed, returning $result \n";color=:red)
+    printstyled(" > Sorry we cannot integrate this expression :(\n";color=:red)
+
     result = integrate_risch(f.val, x.val; kwargs...)
     !contains_int(result) && return result
 
     # TODO make them verbose?
     printstyled("\n > RischMethod(use_algebraic_closure=false, catch_errors=true) failed, returning $result \n";color=:red)
     printstyled(" > Trying with RuleBasedMethod(use_gamma=false, verbose=true)...\n\n"; color=:red)
-    result = integrate_rule_based(f, x; kwargs...)
-    !contains_int(result) && return result
-
-    printstyled(" > RuleBasedMethod(use_gamma=false, verbose=true) failed, returning $result \n";color=:red)
-    printstyled(" > Sorry we cannot integrate this expression :(\n";color=:red)
+
 end
 
 """
@@ -113,7 +115,7 @@ function integrate(f::Symbolics.Num, method=nothing; kwargs...)
         @warn "Multiple symbolic variables detect. Please pass the integration variable to the `integrate` function as second argument."
         return nothing
     elseif length(vars) == 1
-        integration_variable = vars[1]
+        integration_variable = first(vars)
     else
         @warn "No integration variable provided"
         return nothing

src/methods/risch/frontend.jl

@@ -160,7 +160,7 @@ function convolution(a::Vector{T}, b::Vector{T}, s::Int; output_size::Int=0) whe
     c
 end
 
-function tan2sincos(f::K, arg::SymbolicUtils.Symbolic, vars::Vector, h::Int=0) where 
+function tan2sincos(f::K, arg::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}, vars::Vector, h::Int=0) where 
     {T<:FieldElement, P<:PolyRingElem{T}, K<:FracElem{P}}
     # This function transforms a Nemo/AbstractAlgebra rational function with
     # variable t representing tan(arg) to a SymbolicUtils expression which is
@@ -269,13 +269,13 @@ end
 
 struct UpdatedArgList <: Exception end
 
-function analyze_expr(f, x::SymbolicUtils.Symbolic)
+function analyze_expr(f, x::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal})
     tanArgs = []
     expArgs = []
     restart = true
     while restart        
         funs = Any[x]
-        vars = SymbolicUtils.Symbolic[x]
+        vars = SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}[x]
         args = Any[x]
         try 
             p = analyze_expr(f, funs, vars, args, tanArgs, expArgs)
@@ -290,7 +290,7 @@ function analyze_expr(f, x::SymbolicUtils.Symbolic)
     end
 end
 
-function analyze_expr(f::SymbolicUtils.Symbolic , funs::Vector, vars::Vector{SymbolicUtils.Symbolic}, 
+function analyze_expr(f::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal} , funs::Vector, vars::Vector{SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}}, 
                       args::Vector, tanArgs::Vector, expArgs::Vector)
     # Handle pure symbols
     if SymbolicUtils.issym(f)
@@ -305,12 +305,12 @@ function analyze_expr(f::SymbolicUtils.Symbolic , funs::Vector, vars::Vector{Sym
         op = operation(f)
         if op in (+, *, /)
             # Handle Add, Mul, Div operations
-            as = arguments(f)
+            as = [SymbolicUtils.unwrap_const(x) for x in arguments(f)]
             ps = [analyze_expr(a, funs, vars, args, tanArgs, expArgs) for a in as]
             return op(ps...)
         elseif op == (^)
             # Handle Pow operations  
-            as = arguments(f)
+            as = [SymbolicUtils.unwrap_const(x) for x in arguments(f)]
             p1 = analyze_expr(as[1], funs, vars, args, tanArgs, expArgs)
             p2 = analyze_expr(as[2], funs, vars, args, tanArgs, expArgs)
             if isa(p2, Integer)
@@ -431,21 +431,21 @@ function analyze_expr(f::SymbolicUtils.Symbolic , funs::Vector, vars::Vector{Sym
     end
 end
 
-function analyze_expr(f::Number , funs::Vector, vars::Vector{SymbolicUtils.Symbolic}, args::Vector, tanArgs::Vector, expArgs::Vector)
+function analyze_expr(f::Number , funs::Vector, vars::Vector{SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}}, args::Vector, tanArgs::Vector, expArgs::Vector)
     # TODO distinguish types of number (rational, real,  complex, etc. )
     return f
 end
 
 # Commented out - these types don't exist in SymbolicUtils 3.x, handled by generic method above
 # function analyze_expr(f::Union{SymbolicUtils.Add, SymbolicUtils.Mul, SymbolicUtils.Div}, funs::Vector, 
-#                       vars::Vector{SymbolicUtils.Symbolic}, args::Vector, tanArgs::Vector, expArgs::Vector)
+#                       vars::Vector{SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}}, args::Vector, tanArgs::Vector, expArgs::Vector)
 #     as = arguments(f)
 #     ps = [analyze_expr(a, funs, vars, args, tanArgs, expArgs) for a in as]
 #     operation(f)(ps...) # apply f
 # end
 
 # Commented out - SymbolicUtils.Pow doesn't exist in SymbolicUtils 3.x, handled by generic method above
-# function analyze_expr(f::SymbolicUtils.Pow, funs::Vector, vars::Vector{SymbolicUtils.Symbolic}, args::Vector, tanArgs::Vector, expArgs::Vector)
+# function analyze_expr(f::SymbolicUtils.Pow, funs::Vector, vars::Vector{SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}}, args::Vector, tanArgs::Vector, expArgs::Vector)
 #     as = arguments(f)
 #     p1 = analyze_expr(as[1], funs, vars, args, tanArgs, expArgs)
 #     p2 = analyze_expr(as[2], funs, vars, args, tanArgs, expArgs)
@@ -463,10 +463,10 @@ is_rational_multiple(a::SymbolicUtils.Mul, b::SymbolicUtils.Mul) =
     a.dict == b.dict && (isa(a.coeff, Integer) || isa(a.coeff, Rational)) &&
                         (isa(b.coeff, Integer) || isa(b.coeff, Rational))
 
-is_rational_multiple(a::SymbolicUtils.Mul, b::SymbolicUtils.Symbolic) =
+is_rational_multiple(a::SymbolicUtils.Mul, b::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}) =
     (isa(a.coeff, Integer) || isa(a.coeff, Rational)) && (b in keys(a.dict)) && isone(a.dict[b])
 
-is_rational_multiple(a::SymbolicUtils.Symbolic, b::SymbolicUtils.Mul) =
+is_rational_multiple(a::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}, b::SymbolicUtils.Mul) =
     (isa(b.coeff, Integer) || isa(b.coeff, Rational)) && (a in keys(b.dict)) && isone(b.dict[a])
 
 rational_multiple(a, b) = error("not a rational multiple")
@@ -478,14 +478,14 @@ function rational_multiple(a::SymbolicUtils.Mul, b::SymbolicUtils.Mul)
     a.coeff//b.coeff
 end
 
-function rational_multiple(a::SymbolicUtils.Mul, b::SymbolicUtils.Symbolic) 
+function rational_multiple(a::SymbolicUtils.Mul, b::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}) 
     if !is_rational_multiple(a, b)
         error("not a rational multiple")        
     end
     return a.coeff//1
 end
 
-function rational_multiple(a::SymbolicUtils.Symbolic, b::SymbolicUtils.Mul)
+function rational_multiple(a::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}, b::SymbolicUtils.Mul)
     if !is_rational_multiple(a, b)
         error("not a rational multiple")        
     end
@@ -510,7 +510,7 @@ function tan_nx(n::Int, x)
     sign_n*a/b
 end
 
-function analyze_expr(f::SymbolicUtils.Term , funs::Vector, vars::Vector{SymbolicUtils.Symbolic}, args::Vector, tanArgs::Vector, expArgs::Vector)    
+function analyze_expr(f::SymbolicUtils.Term , funs::Vector, vars::Vector{SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}}, args::Vector, tanArgs::Vector, expArgs::Vector)    
     op = operation(f)
     a = arguments(f)[1]
     if op == exp
@@ -597,7 +597,7 @@ function analyze_expr(f::SymbolicUtils.Term , funs::Vector, vars::Vector{Symboli
 end
 
 # Generic method for SymbolicUtils 3.x - handles all symbolic expressions
-function transform_symtree_to_mpoly(f::SymbolicUtils.Symbolic, vars::Vector, vars_mpoly::Vector)
+function transform_symtree_to_mpoly(f::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}, vars::Vector, vars_mpoly::Vector)
     # Handle pure symbols
     if SymbolicUtils.issym(f)
         i = findfirst(x -> hash(x)==hash(f), vars)
@@ -609,17 +609,17 @@ function transform_symtree_to_mpoly(f::SymbolicUtils.Symbolic, vars::Vector, var
     op = operation(f)
     if op == (+)
         # Handle Add operations
-        return sum([transform_symtree_to_mpoly(a, vars, vars_mpoly) for a in arguments(f)])
+        return sum([transform_symtree_to_mpoly(SymbolicUtils.unwrap_const(a), vars, vars_mpoly) for a in arguments(f)])
     elseif op == (*)
         # Handle Mul operations
-        return prod([transform_symtree_to_mpoly(a, vars, vars_mpoly) for a in arguments(f)])
+        return prod([transform_symtree_to_mpoly(SymbolicUtils.unwrap_const(a), vars, vars_mpoly) for a in arguments(f)])
     elseif op == (/)
         # Handle Div operations
-        as = arguments(f)
+        as = [SymbolicUtils.unwrap_const(x) for x in arguments(f)]
         return transform_symtree_to_mpoly(as[1], vars, vars_mpoly)//transform_symtree_to_mpoly(as[2], vars, vars_mpoly)
     elseif op == (^)
         # Handle Pow operations
-        as = arguments(f)
+        as = [SymbolicUtils.unwrap_const(x) for x in arguments(f)]
         @assert isa(as[2], Integer)
         if as[2]>=0
             return transform_symtree_to_mpoly(as[1], vars, vars_mpoly)^as[2]
@@ -768,7 +768,7 @@ end
 
 struct AlgebraicNumbersInvolved <: Exception end
 
-function integrate_risch(f::SymbolicUtils.Add, x::SymbolicUtils.Symbolic; useQQBar::Bool=false,
+function integrate_risch(f::SymbolicUtils.Add, x::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}; useQQBar::Bool=false,
     catchNotImplementedError::Bool=true, catchAlgorithmFailedError::Bool=true)
     # For efficiency compute integral of sum as sum of integrals
     g = f.coeff*x
@@ -779,7 +779,7 @@ function integrate_risch(f::SymbolicUtils.Add, x::SymbolicUtils.Symbolic; useQQB
     g
 end
 
-function integrate_risch(f::SymbolicUtils.Symbolic, x::SymbolicUtils.Symbolic; useQQBar::Bool=false,
+function integrate_risch(f::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}, x::SymbolicUtils.BasicSymbolic{SymbolicUtils.SymReal}; useQQBar::Bool=false,
     catchNotImplementedError::Bool=true, catchAlgorithmFailedError::Bool=true)
     try
         p, funs, vars, args = analyze_expr(f, x)

src/methods/rule_based/frontend.jl

@@ -9,7 +9,7 @@ if not, (original problem, false)
 function apply_rule(problem)
     result = nothing
     for (i, rule) in enumerate(RULES)
-        result = rule(problem)
+        result = rule2(rule, problem)
         if result !== nothing
             if result===problem
                 VERBOSE && println("Infinite cycle created by rule $(IDENTIFIERS[i]) applied on ", problem)
@@ -48,9 +48,13 @@ x^-3
 ```
 """
 function ins(expr)
+    println("called ins with $expr, ",typeof(expr))
     t = (@rule (~u)^(~m) => ~)(expr)
-    t!==nothing && return SymbolicUtils.Term{Number}(^,[t[:u],-t[:m]])
-    return SymbolicUtils.Term{Number}(^,[expr,-1])
+    println("t is $t")
+    t!==nothing && return SymbolicUtils.Term{SymbolicUtils.SymReal}(^,[t[:u],-t[:m]])
+    tmp = SymbolicUtils.Term{SymReal}(^,[expr,-1])
+    println("the return is $tmp")
+    return SymbolicUtils.Term{SymbolicUtils.SymReal}(^,[expr,-1])
 end
 
 # TODO add threaded for speed?
@@ -107,11 +111,11 @@ function repeated_prewalk(expr)
     return expr
 end
 
-function integrate_rule_based(integrand::Symbolics.Num, int_var::Symbolics.Num; use_gamma::Bool=false, verbose::Bool=true, kwargs...)
+function integrate_rule_based(integrand::Symbolics.Num, int_var::Symbolics.Num; use_gamma::Bool=false, verbose::Bool=false, kwargs...)
     global VERBOSE
     VERBOSE = verbose
     return simplify(repeated_prewalk(∫(integrand,int_var)))
 end
 
-integrate_rule_based(integrand::SymbolicUtils.BasicSymbolic{Real}, int_var::SymbolicUtils.BasicSymbolic{Real}; kwargs...) =
-    integrate_rule_based(Num(integrand), Num(int_var); kwargs...)
+# integrate_rule_based(integrand::SymbolicUtils.BasicSymbolic{Real}, int_var::SymbolicUtils.BasicSymbolic{Real}; kwargs...) =
+#     integrate_rule_based(Num(integrand), Num(int_var); kwargs...)

src/methods/rule_based/general.jl

@@ -64,7 +64,7 @@ coshintegral(x::Any) = println("hyperbolic cosine integral Chi(z) function (http
 @syms ∫(var1,var2) substitute_after_int(var1, var2, var3)
 
 # very big arrays containing rules and their identifiers
-const RULES = SymbolicUtils.Rule[]
+const RULES = Pair{Expr, Expr}[]
 const IDENTIFIERS = String[]
 
 # to use or not the gamma function in integration results
@@ -80,19 +80,19 @@ all_rules_paths = [
 "9 Miscellaneous/0.1 Integrand simplification rules.jl"
 
 "1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.1 (a+b x)^m.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.2 (a+b x)^m (c+d x)^n.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.3 (a+b x)^m (c+d x)^n (e+f x)^p.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.4 (a+b x)^m (c+d x)^n (e+f x)^p (g+h x)^q.jl"
-# 
-# "1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.1 (a+b x^2)^p.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.2 (c x)^m (a+b x^2)^p.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.3 (a+b x^2)^p (c+d x^2)^q.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.4 (e x)^m (a+b x^2)^p (c+d x^2)^q.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.2 (a+b x)^m (c+d x)^n.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.3 (a+b x)^m (c+d x)^n (e+f x)^p.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.4 (a+b x)^m (c+d x)^n (e+f x)^p (g+h x)^q.jl"
+# 
+"1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.1 (a+b x^2)^p.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.2 (c x)^m (a+b x^2)^p.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.3 (a+b x^2)^p (c+d x^2)^q.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.4 (e x)^m (a+b x^2)^p (c+d x^2)^q.jl"
 # # 5, 6, 7, 8, 9
 # 
-# "1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.1 (a+b x^n)^p.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.2 (c x)^m (a+b x^n)^p.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.3 (a+b x^n)^p (c+d x^n)^q.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.1 (a+b x^n)^p.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.2 (c x)^m (a+b x^n)^p.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.3 (a+b x^n)^p (c+d x^n)^q.jl"
 # "1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.4 (e x)^m (a+b x^n)^p (c+d x^n)^q.jl"
 # "1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.5 (a+b x^n)^p (c+d x^n)^q (e+f x^n)^r.jl"
 # "1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.6 (g x)^m (a+b x^n)^p (c+d x^n)^q (e+f x^n)^r.jl"
@@ -120,8 +120,8 @@ all_rules_paths = [
 # # 1.4.1, 1.4.2
 # 
 # "1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.7 P(x) (a+b x)^m (c+d x)^n (e+f x)^p (g+h x)^q.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.6 P(x) (a+b x)^m (c+d x)^n (e+f x)^p.jl"
-# "1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.5 P(x) (a+b x)^m (c+d x)^n.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.6 P(x) (a+b x)^m (c+d x)^n (e+f x)^p.jl"
+"1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.5 P(x) (a+b x)^m (c+d x)^n.jl"
 # 
 # "1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/1.2.2.5 P(x) (a+b x^2+c x^4)^p.jl"
 # "1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/1.2.2.4 (f x)^m (d+e x^2)^q (a+b x^2+c x^4)^p.jl"
@@ -133,7 +133,7 @@ all_rules_paths = [
 # 
 # 
 # 
-# "2 Exponentials/2.1 (c+d x)^m (a+b (F^(g (e+f x)))^n)^p.jl"
+"2 Exponentials/2.1 (c+d x)^m (a+b (F^(g (e+f x)))^n)^p.jl"
 # "2 Exponentials/2.2 (c+d x)^m (F^(g (e+f x)))^n (a+b (F^(g (e+f x)))^n)^p.jl"
 # "2 Exponentials/2.3 Miscellaneous exponentials.jl"
 #