RuleBasedIntegration / Rubi

Rubi for Mathematica
http://rulebasedintegration.org
MIT License
228 stars 23 forks source link

Questions on some rules (which came up in preparation of a Maxima port) #53

Open MSoegtropIMC opened 3 months ago

MSoegtropIMC commented 3 months ago

Dear Rubi Team,

I started a port of Rubi to Maxima, which is based on the idea of the late Albert Rich to automatically convert the rule base to a decision tree - which should then be easily portable to various systems including Maxima. The plan is to write the converter in OCaml, because I think it is easier to do this in a strictly typed language supporting algebraic data types. Also I plan to port it to Coq some day.

As first step I wrote a parser for the subset of Mathematica used in Rubi rules in Menhir (OCaml's yacc). I know that Mathematica parsers exist, but it is easier this way and it also helps me to find out what subset of Mathematica is actually used in Rubi.

I came across a few rules which look questionable to me.

First there are quite a few rules which use \[Star] where I would expect ASCII * aka. Times. E.g.:

https://github.com/RuleBasedIntegration/Rubi/blob/61e9c18ea248061cd83c67882f7c91a73cef912d/Rubi/IntegrationRules/1%20Algebraic%20functions/1.1%20Binomial%20products/1.1.2%20Quadratic/1.1.2.1%20(a%2Bb%20x%5E2)%5Ep.m#L24

Is this a typo (when the file is loaded, the rendering is hard to distinguish from ASCII *) or is this some trick I need to understand? In the latter case, can you please explain it?

Second there are a few rules which use float numbers while almost all rules user integers / rationals: E.g. the x_^2. in:

https://github.com/RuleBasedIntegration/Rubi/blob/61e9c18ea248061cd83c67882f7c91a73cef912d/Rubi/IntegrationRules/1%20Algebraic%20functions/1.1%20Binomial%20products/1.1.2%20Quadratic/1.1.2.4%20(e%20x)%5Em%20(a%2Bb%20x%5E2)%5Ep%20(c%2Bd%20x%5E2)%5Eq.m#L94

Is this intentional?

Best regards,

Michael

stblake commented 3 months ago

Hi Michael,

It's great to hear that someone is taking up the challenge of compiling Rubi's rules into a decision tree!

Regarding [Star] - I believe it was Albert's way of having an inert Times[] to prevent things like partial simplification/expansion, for example

In[878]:= -1*(a + b)

Out[878]= -a - b

So Star[] should probably parse to Times[]. It will be interesting to see what bugs this introduces/eliminates.

It is hard to know if patterns containing Real numbers are intentional or accidental. I would guess accidental as it wouldn't fit with the rest of the design of the Rubi rules. All this guesswork is because Rubi was poorly documented.

You should also be aware that Nasser has a list of Rubi rules which cause infinite recursion. Almost certainly these rules should be removed from Rubi prior to compiling your decision tree.

Cheers,

Sam Blake

On Thu, Jun 20, 2024 at 5:50 AM MSoegtropIMC @.***> wrote:

Dear Rubi Team,

I started a port of Rubi to Maxima, which is based on the idea of the late Albert Rich to automatically convert the rule base to a decision tree - which should then be easily portable to various systems including Maxima. The plan is to write the converter in OCaml, because I think it is easier to do this in a strictly typed language supporting algebraic data types. Also I plan to port it to Coq some day.

As first step I wrote a parser for the subset of Mathematica used in Rubi rules in Menhir (OCaml's yacc). I know that Mathematica parsers exist, but it is easier this way and it also helps me to find out what subset of Mathematica is actually used in Rubi.

I came across a few rules which look questionable to me.

First there are quite a few rules which use [Star] where I would expect ASCII * aka. Times. E.g.:

https://github.com/RuleBasedIntegration/Rubi/blob/61e9c18ea248061cd83c67882f7c91a73cef912d/Rubi/IntegrationRules/1%20Algebraic%20functions/1.1%20Binomial%20products/1.1.2%20Quadratic/1.1.2.1%20(a%2Bb%20x%5E2)%5Ep.m#L24

Is this a typo (when the file is loaded, the rendering is hard to distinguish from ASCII *) or is this some trick I need to understand? In the latter case, can you please explain it?

Second there are a few rules which use float numbers while almost all rules user integers / rationals: E.g. the x_^2. in:

https://github.com/RuleBasedIntegration/Rubi/blob/61e9c18ea248061cd83c67882f7c91a73cef912d/Rubi/IntegrationRules/1%20Algebraic%20functions/1.1%20Binomial%20products/1.1.2%20Quadratic/1.1.2.4%20(e%20x)%5Em%20(a%2Bb%20x%5E2)%5Ep%20(c%2Bd%20x%5E2)%5Eq.m#L94

Is this intentional?

Best regards,

Michael

— Reply to this email directly, view it on GitHub https://github.com/RuleBasedIntegration/Rubi/issues/53, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB5ALCRW6ND4G36T5EMHOALZIHOKBAVCNFSM6AAAAABJST2JW2VHI2DSMVQWIX3LMV43ASLTON2WKOZSGM3DGMBTG43DSOI . You are receiving this because you are subscribed to this thread.Message ID: @.***>

arthurcnorman commented 3 months ago

See https://reference.wolfram.com/language/ref/character/Star.html for what Mathematica says [Star] means. How that participated in sunsequent Mathematica transformations is not my expertise, but in my activity parting the Rubi rulke-set for Reduce I encountered the same syntax. If you are working on a Maxima port and I am in parallel looking at Reduce then it could well make sense for us to exchange experiences and problems encountered - I have a medium length list of thoughts. But the level of detail in what we might be exchanging would seem to me to mean that would be best done by direct email between us not through this forum.

I can parse all the Rubi rules (and indeed test examples) from the mathematica syntax they come in into LIsp prefix form with operators such as "/;" and so on. I have been experimenting with various transformations there but the details of where I am are very fluid and are messy enough that I am not going to post details here to bore everybody else. But if you are active on a project a bit of joining forces might be a good plan? Arthur Norman (find me if you need via the reduce-algebra project on sourceforge).

axkr commented 3 months ago

I started a port of Rubi to Maxima, which is based on the idea of the late Albert Rich to automatically convert the rule base to a decision tree - which should then be easily portable to various systems including Maxima.

I'm also planning to create a generator for converting rules into a Java source code method:

I have some uncommitted changes creating a simplified decision tree (discrimination net). I think the Rubi rules are to complicated for the first step and I will try to support only simple pattern-matching rules for "builtin-function" at the beginning.

Maybe we can work together on the level of "best used data structures" or "known papers" for this topic?

MSoegtropIMC commented 3 months ago

@arthurcnorman @axkr : thank you - indeed it makes sense to work together on this!

@arthurcnorman: are you fine with the basic approach of creating a decision tree / discrimination network offline with a common tool and exporting this into various languages? The advantage of this method I see is that it should be easier to port than reimplementing Mathematica pattern matching and it should also be much faster (if one does the creation of the tree/net offline once).

@axkr: I also had the plan to first look into what MatchPy does (the paper you linked), just as an offline solution. The SymPy Rubi port reportedly also uses MatchPy, but it is said that the online generation of the tree / net is too slow. One could think about generating Maxima / Java code using MatchPy, so that the generation time would not be an issue.

Where do we want to continue to discuss this? I think it would make sense to have a public discussion and IMHO it makes sense to have the discussion somewhere under Rubi, but maybe under Rubi-5 since this was the main intention of Rubi-5 and there might be fewer people who are subscribed to all issues.

Does someone from the Rubi team object?

stblake commented 3 months ago

Here's an interesting use of GPT4o parsing a Rubi rule into Maple syntax. (I do not have Maple so I cannot test the resulting code... )

Prompt:

Can you code this integration rule in Maple syntax? Int[(a + b.x + c.x^2)^p, x_Symbol] := 2((a + bx + cx^2)^(p + 1)/ ((2p + 1)(b + 2cx))) /; FreeQ[{a, b, c, p}, x] && b^2 - 4ac == 0 && p <= -1*

Certainly! Here’s how you can implement the given integration rule in Maple syntax. The rule can be encoded as a custom integration procedure using Maple’s pattern matching and conditional logic.

First, let’s restate the rule: [ \int (a + bx + cx^2)^p \, dx = \frac{2(a + bx + cx^2)^{p+1}}{(2p + 1)(b

Maple Code

with(inttrans):
with(int):
IntegrationRule := proc(f, x)
    local a, b, c, p, integrand, integrandPattern, result;

    # Define the pattern for the integrand
    integrandPattern := (a_ + b_ * x + c_ * x^2)^p_;

    # Check if the integrand matches the pattern and conditions
    if type(f, integrandPattern) then
        # Extract the coefficients and exponent
        a := op(1, f);
        b := op(2, f);
        c := op(3, f);
        p := op(4, f);

        # Check the conditions: FreeQ[{a, b, c, p}, x] && b^2 - 4*a*c = 0
&& p <= -1
        if not depends(a, x) and not depends(b, x) and not depends(c, x)
and not depends(p, x) and
           b^2 - 4*a*c = 0 and p <= -1 then
            # Apply the integration rule
            result := 2 * ((a + b * x + c * x^2)^(p + 1)) / ((2 * p + 1) *
(b + 2 * c * x));
            return result;
        else
            return 'unknown';
        end if;
    else
        return 'unknown';
    end if;
end proc;

# Example usage
f := (3 + 6 * x + 3 * x^2)^(-2);
IntegrationRule(f, x);

Explanation:

  1. Load Necessary Packages: The inttrans and int packages are loaded to use integration functions.
  2. Define the Procedure: IntegrationRule is a procedure that takes an integrand f and a variable x.
  3. Pattern Matching: The integrand pattern (a_ + b_ * x + c_ * x^2)^p_ is defined to match the given form.
  4. Condition Checking: The procedure checks if:
    • a, b, c, and p are free of the variable x.
    • The condition ( b^2 - 4ac = 0 ) holds.
    • The exponent ( p \leq -1 ).
  5. Apply the Integration Rule: If all conditions are met, the integration rule is applied and the result is returned. If any condition is not met, the procedure returns 'unknown'.

Example Usage:

This implementation should work for the specified rule and conditions in Maple. If you need further customization or have additional rules to implement, similar procedures can be created following this structure.

On Fri, Jun 21, 2024 at 2:39 AM MSoegtropIMC @.***> wrote:

@arthurcnorman https://github.com/arthurcnorman @axkr https://github.com/axkr : thank you - indeed it makes sense to work together on this!

@arthurcnorman https://github.com/arthurcnorman: are you fine with the basic approach of creating a decision tree / discrimination network offline with a common tool and exporting this into various languages? The advantage of this method I see is that it should be easier to port than reimplementing Mathematica pattern matching and it should also be much faster (if one does the creation of the tree/net offline once).

@axkr https://github.com/axkr: I also had the plan to first look into what MatchPy does (the paper you linked), just as an offline solution. The SymPy Rubi port reportedly also uses MatchPy, but it is said that the online generation of the tree / net is too slow. One could think about generating Maxima / Java code using MatchPy, so that the generation time would not be an issue.

Where do we want to continue to discuss this? I think it would make sense to have a public discussion and IMHO it makes sense to have the discussion somewhere under Rubi, but maybe under Rubi-5 since this was the main intention of Rubi-5 and there might be fewer people who are subscribed to all issues.

Does someone from the Rubi team object?

— Reply to this email directly, view it on GitHub https://github.com/RuleBasedIntegration/Rubi/issues/53#issuecomment-2181114566, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB5ALCWDQSD5PC7PVBTVAA3ZIMAT7AVCNFSM6AAAAABJST2JW2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCOBRGEYTINJWGY . You are receiving this because you commented.Message ID: @.***>

stblake commented 3 months ago

It looks dubious to me, but I don't have Maple any more. The online documentation says about the op() function ....

If the first argument to the op function is a positive integer i, then the ith operand of e is returned.

In the example, it seems that the ChatGPT program uses if type(f, integrandPattern) then

Extract the coefficients and exponent

    a := op(1, f);
    b := op(2, f);
    c := op(3, f);

    p := op(4, f)

... to match the pattern.

is this the way matching works in Maple???

My own experiments (with Copilot, mostly) suggest that it is happy to provide

buggy programs (in Lisp!) which have a whiff of some nearby idea.

On Thu, Jun 20, 2024 at 4:52 PM Sam Blake @.***> wrote:

Here's an interesting use of GPT4o parsing a Rubi rule into Maple syntax. (I do not have Maple so I cannot test the resulting code... )

Prompt:

Can you code this integration rule in Maple syntax? Int[(a + b.x + c.x^2)^p, x_Symbol] := 2((a + bx + cx^2)^(p + 1)/ ((2p + 1)(b + 2cx))) /; FreeQ[{a, b, c, p}, x] && b^2 - 4ac == 0 && p <= -1*

Certainly! Here’s how you can implement the given integration rule in Maple syntax. The rule can be encoded as a custom integration procedure using Maple’s pattern matching and conditional logic.

First, let’s restate the rule: [ \int (a + bx + cx^2)^p \, dx = \frac{2(a + bx + cx^2)^{p+1}}{(2p + 1)(b

  • 2cx)} ] given that ( b^2 - 4ac = 0 ) and ( p \leq -1 ).

Maple Code

with(inttrans):
with(int):
IntegrationRule := proc(f, x)
    local a, b, c, p, integrand, integrandPattern, result;

    # Define the pattern for the integrand
    integrandPattern := (a_ + b_ * x + c_ * x^2)^p_;

    # Check if the integrand matches the pattern and conditions
    if type(f, integrandPattern) then
        # Extract the coefficients and exponent
        a := op(1, f);
        b := op(2, f);
        c := op(3, f);
        p := op(4, f);

        # Check the conditions: FreeQ[{a, b, c, p}, x] && b^2 - 4*a*c = 0
&& p <= -1
        if not depends(a, x) and not depends(b, x) and not depends(c, x)
and not depends(p, x) and
           b^2 - 4*a*c = 0 and p <= -1 then
            # Apply the integration rule
            result := 2 * ((a + b * x + c * x^2)^(p + 1)) / ((2 * p + 1) *
(b + 2 * c * x));
            return result;
        else
            return 'unknown';
        end if;
    else
        return 'unknown';
    end if;
end proc;

# Example usage
f := (3 + 6 * x + 3 * x^2)^(-2);
IntegrationRule(f, x);

Explanation:

  1. Load Necessary Packages: The inttrans and int packages are loaded to use integration functions.
  2. Define the Procedure: IntegrationRule is a procedure that takes an integrand f and a variable x.
  3. Pattern Matching: The integrand pattern (a_ + b_ * x + c_ * x^2)^p_ is defined to match the given form.
  4. Condition Checking: The procedure checks if:
    • a, b, c, and p are free of the variable x.
    • The condition ( b^2 - 4ac = 0 ) holds.
    • The exponent ( p \leq -1 ).
  5. Apply the Integration Rule: If all conditions are met, the integration rule is applied and the result is returned. If any condition is not met, the procedure returns 'unknown'.

Example Usage:

  • An example usage is provided with the integrand ( (3 + 6x + 3x^2)^{-2} ).
  • The procedure IntegrationRule is called with this integrand and the variable x.

This implementation should work for the specified rule and conditions in Maple. If you need further customization or have additional rules to implement, similar procedures can be created following this structure.

On Fri, Jun 21, 2024 at 2:39 AM MSoegtropIMC @.***> wrote:

@arthurcnorman https://github.com/arthurcnorman @axkr https://github.com/axkr : thank you - indeed it makes sense to work together on this!

@arthurcnorman https://github.com/arthurcnorman: are you fine with the basic approach of creating a decision tree / discrimination network offline with a common tool and exporting this into various languages? The advantage of this method I see is that it should be easier to port than reimplementing Mathematica pattern matching and it should also be much faster (if one does the creation of the tree/net offline once).

@axkr https://github.com/axkr: I also had the plan to first look into what MatchPy does (the paper you linked), just as an offline solution. The SymPy Rubi port reportedly also uses MatchPy, but it is said that the online generation of the tree / net is too slow. One could think about generating Maxima / Java code using MatchPy, so that the generation time would not be an issue.

Where do we want to continue to discuss this? I think it would make sense to have a public discussion and IMHO it makes sense to have the discussion somewhere under Rubi, but maybe under Rubi-5 since this was the main intention of Rubi-5 and there might be fewer people who are subscribed to all issues.

Does someone from the Rubi team object?

— Reply to this email directly, view it on GitHub https://github.com/RuleBasedIntegration/Rubi/issues/53#issuecomment-2181114566, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB5ALCWDQSD5PC7PVBTVAA3ZIMAT7AVCNFSM6AAAAABJST2JW2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCOBRGEYTINJWGY . You are receiving this because you commented.Message ID: @.***>

stblake commented 3 months ago

Has anyone here looked at this?

https://github.com/axkr/symja_android_library/wiki/Porting-Rubi-Integration-rules-to-Symja

(It was just pointed out on the Maxima discussion list. Some of Mathematica is implemented.. This page gets you a web interface

https://matheclipse.org/

which answers some questions right. Int[a*x^2,{x,r,s}] works. ) RJF

axkr commented 3 months ago

Has anyone here looked at this?

yes, despite the name, the https://matheclipse.org/ is the Symja demo page. If you know some Java you can install it locally:

stblake commented 3 months ago

Wow, I love the web interface!

I was surprised to see this time-out:

Integrate[1/(x^5+1), x] Error: 𝖳𝗂𝗆𝖾𝗈𝗎𝗍 𝖾𝗑𝖼𝖾𝖾𝖽𝖾𝖽. 𝖢𝖺𝗅𝖼𝗎𝗅𝖺𝗍𝗂𝗈𝗇 𝖺𝖻𝗈𝗋𝗍𝖾𝖽!

On Mon, Jun 24, 2024 at 7:49 AM Axel Kramer @.***> wrote:

Has anyone here looked at this?

yes, despite the name, the https://matheclipse.org/ is the Symja demo page. If you know some Java you can install it locally:

— Reply to this email directly, view it on GitHub https://github.com/RuleBasedIntegration/Rubi/issues/53#issuecomment-2185333316, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB5ALCQYUXQV3RDR73SCHPTZI47HVAVCNFSM6AAAAABJST2JW2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCOBVGMZTGMZRGY . You are receiving this because you commented.Message ID: @.***>

robert-dodier commented 3 months ago

@arthurcnorman About exchanging ideas with others working on integration rules, I would like to very strongly encourage you to conduct the discussion in a public forum, whatever is convenient for you, even if at this time only one other person is participating. Whatever you discover or encounter is potentially useful to others, even if it is not clear to anyone at present what form that will take.

All the best, Robert Dodier (developer and project administrator for Maxima)

axkr commented 3 months ago

I was surprised to see this time-out: Integrate[1/(x^5+1), x] Error: 𝖳𝗂𝗆𝖾𝗈𝗎𝗍 𝖾𝗑𝖼𝖾𝖾𝖽𝖾𝖽. 𝖢𝖺𝗅𝖼𝗎𝗅𝖺𝗍𝗂𝗈𝗇 𝖺𝖻𝗈𝗋𝗍𝖾𝖽!

Yes the demo on google appengine is very slow and times out after 30 seconds.

Locally on my machine, I get this result: image

axkr commented 3 months ago

@arthurcnorman About exchanging ideas with others working on integration rules, I would like to very strongly encourage you to conduct the discussion in a public forum, whatever is convenient for you, even if at this time only one other person is participating. Whatever you discover or encounter is potentially useful to others, even if it is not clear to anyone at present what form that will take.

the maintainers of the github Rubi organization can enable "discussions" for this purpose in their organization:

or in a specific Rubi repo:

For chatting I personally prefer Discord nowadays:

stblake commented 3 months ago

Suggestion for speeding up (Mathematica's version) of Rubi by maybe 30-40%..

It seemed to me that pattern matching should be way faster when there are better "anchors" in the match. Thus if we substitute a particular name, say VAR for the variable of integration, all the rules are simpler. A rule that now looks like

Int [a.*x^2+ ....., xSymbol] ..... could be replaced by Int[a.* VAR^2+.....]

In some simple tests I get about 33% faster, though I am just testing matching speed. (this trick is used in Moses's SIN integration system in Maxima. SIN relies heavily on pattern matching.) RJF

Message ID: @.***>

stblake commented 3 months ago

Albert probably would probably have hated the inelegant hackery, but perhaps another way to do that would be have only one public rule for

Int [expression_? ExpressionQ, var_SymbolQ] := IntAux [expression, var],

then change the name of the other rules to private IntAux which would then not have to verify that the first argument is not, for example a Boolen or the second argument an Integer.

On Tue, Jun 25, 2024 at 8:05 AM Richard Fateman @.***> wrote:

Suggestion for speeding up (Mathematica's version) of Rubi by maybe 30-40%..

It seemed to me that pattern matching should be way faster when there are better "anchors" in the match. Thus if we substitute a particular name, say VAR for the variable of integration, all the rules are simpler. A rule that now looks like

Int [a.*x^2+ ....., xSymbol] ..... could be replaced by Int[a.* VAR^2+.....]

In some simple tests I get about 33% faster, though I am just testing matching speed. (this trick is used in Moses's SIN integration system in Maxima. SIN relies heavily on pattern matching.) RJF

Message ID: @.***>

stblake commented 3 months ago

On Tue, Jun 25, 2024 at 11:41 AM David Stoutemyer @.***> wrote:

Albert probably would probably have hated the inelegant hackery, but perhaps another way to do that would be have only one public rule for

Int [expression_? ExpressionQ, var_SymbolQ] := IntAux [expression, var],

I'm suggesting a bigger improvement in speed, at least in principle. Yes, those IntAux rules would not have to check (again) that var is a symbol, but those rules would still have to check the names. Thus making up a rule as an example... for x^2 -> x^3/3....

IntAux[ a.*var^2+b,var] := a*var^3/3+b /; FreeQ[{a,b},var]

what if the application is

IntAux[ u^2+5v^2+7w^2+ 2x^2+4y^2+5*z^2, w]

the pattern matcher, unless clever enough to match args out of order and set var==w first, will wander around trying var==u, var==v, var... etc.

So it is better to do this:

Int [expression_? ExpressionQ, var_SymbolQ] := IntAux [expression /. var-> VAR]

and then IntAux would see u^2+5v^2+7VAR^2+ 2x^2+4y^2+5z^2. Much easier to match with IntAux[ a_.VAR^2+b_] := a*VAR^3/3+b /; FreeQ[{a,b},VAR]

one glitch .. if we generate some error message from IntAux and friends, those programs will not know the external name of the variable of integration. They will only know VAR. So we have to save the external name. Easy to do at top level..

Int [expression_? ExpressionQ, var_SymbolQ] := (SavedVAR=var ; IntAux [expression /. var-> VAR])

.............. RJF

then change the name of the other rules to private IntAux which would then

not have to verify that the first argument is not, for example a Boolen or the second argument an Integer.

On Tue, Jun 25, 2024 at 8:05 AM Richard Fateman @.***> wrote:

Suggestion for speeding up (Mathematica's version) of Rubi by maybe 30-40%..

It seemed to me that pattern matching should be way faster when there are better "anchors" in the match. Thus if we substitute a particular name, say VAR for the variable of integration, all the rules are simpler. A rule that now looks like

Int [a.*x^2+ ....., xSymbol] ..... could be replaced by Int[a.* VAR^2+.....]

In some simple tests I get about 33% faster, though I am just testing matching speed. (this trick is used in Moses's SIN integration system in Maxima. SIN relies heavily on pattern matching.) RJF

Message ID: @.***>

arthurcnorman commented 3 months ago

Was it David who commented on what Albert might have hated? Well my understanding is that Albert liked the concept of mapping from the rule-set to a decision tree. So I would feel that his legacy would be honoured by almost any automated conversion so that the original rule-set remained unaltered but the version fed into a matching engine had been transformed in various ways. But systematiaclly and mechanically and in ways that did not alter the semantics. An implmentation for say Maxima etc certainly does that starting with its own parser for the Mathematica syntax and a non-mathematica pattern match engine.

A quite separate issue will be long term maintenance and enhancement of the main rule-set. It was still work in progress whan Albert's work on it ceased so there will be more to do. But piecemeal hacks are not liable to be useful ways to support long life there.

My interest is in potential implementation outside Mathematica. The SYNTAX side is - as we have seen - pretty routine. The Semantics is to my mind a quite different matter. As seen from the 30% cost change from basically trivial adjustments, just how Mathematica pattern matching works and so exactly what it will always do is maybe a bit obcure around the edges. Many Rubi rules use Simplify[] and the "specification" of that is that it tries many (unspecified) transformations on its argument and returns the "simplest" where even that term is not really defined. Then Rubi comes with 4000 lines or so of Mathematica-coded support code. Some of that looks as it does obvious or simple things but replicating other bits outside Mathematica asks for reverse engineering exact mathematica behaviour and re-building it portabily. A serious task. I find mathematica behaviour a bit odd. So eg NumberQ[1/2] is true but LeafCount[1/2] = 3 so the first sort of says it is a single thing while the second says it is composite.

My conclusion is that at present the only stable definition of what Rubi is is the EXACT way it behaves on (the current release of) Mathematica. The rest of us can use the rule set to obtain what may end up as still interesting integration capability, but what we get will not be exactly Rubi. But seeing what the divergences are will be interesting. Then as a very separate activity there can be hand adjustment to the rule set. I do not guarantee that messing with the exact rule for the variable of integration will speed up my version. It could even slow it down! So fixing actual bugs and enhancing capability would come first for me.

End of rant.

Arthur

On Tue, 25 Jun 2024, Sam Blake wrote:

On Tue, Jun 25, 2024 at 11:41 AM David Stoutemyer @.***> wrote:

Albert probably would probably have hated the inelegant hackery, but perhaps another way to do that would be have only one public rule for

Int [expression_? ExpressionQ, var_SymbolQ] := IntAux [expression, var],

I'm suggesting a bigger improvement in speed, at least in principle. Yes, those IntAux rules would not have to check (again) that var is a symbol, but those rules would still have to check the names. Thus making up a rule as an example... for x^2 -> x^3/3....

IntAux[ a.*var^2+b,var] := a*var^3/3+b /; FreeQ[{a,b},var]

what if the application is

IntAux[ u^2+5v^2+7w^2+ 2x^2+4y^2+5*z^2, w]

the pattern matcher, unless clever enough to match args out of order and set var==w first, will wander around trying var==u, var==v, var... etc.

So it is better to do this:

Int [expression_? ExpressionQ, var_SymbolQ] := IntAux [expression /. var-> VAR]

and then IntAux would see u^2+5v^2+7VAR^2+ 2x^2+4y^2+5z^2. Much easier to match with IntAux[ a_.VAR^2+b_] := a*VAR^3/3+b /; FreeQ[{a,b},VAR]

one glitch .. if we generate some error message from IntAux and friends, those programs will not know the external name of the variable of integration. They will only know VAR. So we have to save the external name. Easy to do at top level..

Int [expression_? ExpressionQ, var_SymbolQ] := (SavedVAR=var ; IntAux [expression /. var-> VAR])

.............. RJF

then change the name of the other rules to private IntAux which would then

not have to verify that the first argument is not, for example a Boolen or the second argument an Integer.

On Tue, Jun 25, 2024 at 8:05 AM Richard Fateman @.***> wrote:

Suggestion for speeding up (Mathematica's version) of Rubi by maybe 30-40%..

It seemed to me that pattern matching should be way faster when there are better "anchors" in the match. Thus if we substitute a particular name, say VAR for the variable of integration, all the rules are simpler. A rule that now looks like

Int [a.*x^2+ ....., xSymbol] ..... could be replaced by Int[a.* VAR^2+.....]

In some simple tests I get about 33% faster, though I am just testing matching speed. (this trick is used in Moses's SIN integration system in Maxima. SIN relies heavily on pattern matching.) RJF

Message ID: @.***>

-- Reply to this email directly or view it on GitHub: https://github.com/RuleBasedIntegration/Rubi/issues/53#issuecomment-2189775698 You are receiving this because you were mentioned.

Message ID: @.***> --3735943886-1724912013-1719352376=:27136--

stblake commented 3 months ago

Subtle differences between Mathematica Rubi, Maxima Rubi, Reduce Rubi, etc are to be expected. Similar to the (usually not so) subtle differences in the Risch algorithm implementations in each CAS.

Richard, I see about a 10% improvement on a Rubi rule I picked out at random.

In[3946]:= ClearAll[integrale];

integrale[(a + b.x + c.x^2)^p(d + f.x^2)^q(g. + h.x_), x_Symbol] := h(a + bx + cx^2)^ p(d + fx^2)^(q + 1)/(2f(p + q + 1)) - (1/(2f(p + q + 1))) integrale[(a + bx + cx^2)^(p - 1)(d + fx^2)^q hp(bd) + a(-2gf)(p + q + 1) + (2hp(cd - af) + b(-2gf)(p + q + 1)) x + (hp(-bf) + c(-2gf)(p + q + 1))x^2, x] /; FreeQ[{a, b, c, d, f, g, h, q}, x] && b^2 - 4ac != 0 && p > 0 && p + q + 1 != 0

In[3948]:= integrale[(1 + 3 x + 4 x^2)^(3/2)*(-2 + x^2)^(5/2) (1 - x), x] // RepeatedTiming

Out[3948]= {0.0000586511, -(1/10) (-2 + x^2)^(7/2) (1 + 3 x + 4 x^2)^( 3/2) - 1/ 10 integrale[-10 - 3 x - (71 x^2)/2 + 9 (-2 + x^2)^(5/2) Sqrt[1 + 3 x + 4 x^2], x]}

In[3949]:= ClearAll[integrale];

integrale[(a + b.X + c_.X^2)^p*(d + f.*X^2)^q(g. + h.X)] := h(a + bX + cX^2)^ p(d + fX^2)^(q + 1)/(2f(p + q + 1)) - (1/(2f(p + q + 1))) integrale[(a + bX + cX^2)^(p - 1)(d + fX^2)^q hp(bd) + a(-2gf)(p + q + 1) + (2hp(cd - af) + b(-2gf)(p + q + 1)) X + (hp(-bf) + c(-2gf)(p + q + 1))X^2, X] /; FreeQ[{a, b, c, d, f, g, h, q}, X] && b^2 - 4ac != 0 && p > 0 && p + q + 1 != 0

In[3951]:= integrale[(1 + 3 X + 4 X^2)^(3/2)*(-2 + X^2)^(5/2) (1 - X)] // RepeatedTiming

Out[3951]= {0.0000518786, -(1/10) (-2 + X^2)^(7/2) (1 + 3 X + 4 X^2)^( 3/2) - 1/ 10 integrale[-10 - 3 X - (71 X^2)/2 + 9 (-2 + X^2)^(5/2) Sqrt[1 + 3 X + 4 X^2], X]}

In[3952]:= (0.00005187860107421875- 0.00005865106201171875)/0.00005865106201171875`

Out[3952]= -0.11547

Cheers,

Sam

On Wed, Jun 26, 2024 at 7:53 AM Arthur Norman @.***> wrote:

Was it David who commented on what Albert might have hated? Well my understanding is that Albert liked the concept of mapping from the rule-set to a decision tree. So I would feel that his legacy would be honoured by almost any automated conversion so that the original rule-set remained unaltered but the version fed into a matching engine had been transformed in various ways. But systematiaclly and mechanically and in ways that did not alter the semantics. An implmentation for say Maxima etc certainly does that starting with its own parser for the Mathematica syntax and a non-mathematica pattern match engine.

A quite separate issue will be long term maintenance and enhancement of the main rule-set. It was still work in progress whan Albert's work on it ceased so there will be more to do. But piecemeal hacks are not liable to be useful ways to support long life there.

My interest is in potential implementation outside Mathematica. The SYNTAX side is - as we have seen - pretty routine. The Semantics is to my mind a quite different matter. As seen from the 30% cost change from basically trivial adjustments, just how Mathematica pattern matching works and so exactly what it will always do is maybe a bit obcure around the edges. Many Rubi rules use Simplify[] and the "specification" of that is that it tries many (unspecified) transformations on its argument and returns the "simplest" where even that term is not really defined. Then Rubi comes with 4000 lines or so of Mathematica-coded support code. Some of that looks as it does obvious or simple things but replicating other bits outside Mathematica asks for reverse engineering exact mathematica behaviour and re-building it portabily. A serious task. I find mathematica behaviour a bit odd. So eg NumberQ[1/2] is true but LeafCount[1/2] = 3 so the first sort of says it is a single thing while the second says it is composite.

My conclusion is that at present the only stable definition of what Rubi is is the EXACT way it behaves on (the current release of) Mathematica. The rest of us can use the rule set to obtain what may end up as still interesting integration capability, but what we get will not be exactly Rubi. But seeing what the divergences are will be interesting. Then as a very separate activity there can be hand adjustment to the rule set. I do not guarantee that messing with the exact rule for the variable of integration will speed up my version. It could even slow it down! So fixing actual bugs and enhancing capability would come first for me.

End of rant.

Arthur

On Tue, 25 Jun 2024, Sam Blake wrote:

On Tue, Jun 25, 2024 at 11:41 AM David Stoutemyer @.***> wrote:

Albert probably would probably have hated the inelegant hackery, but perhaps another way to do that would be have only one public rule for

Int [expression_? ExpressionQ, var_SymbolQ] := IntAux [expression, var],

I'm suggesting a bigger improvement in speed, at least in principle. Yes, those IntAux rules would not have to check (again) that var is a symbol, but those rules would still have to check the names. Thus making up a rule as an example... for x^2 -> x^3/3....

IntAux[ a.*var^2+b,var] := a*var^3/3+b /; FreeQ[{a,b},var]

what if the application is

IntAux[ u^2+5v^2+7w^2+ 2x^2+4y^2+5*z^2, w]

the pattern matcher, unless clever enough to match args out of order and set var==w first, will wander around trying var==u, var==v, var... etc.

So it is better to do this:

Int [expression_? ExpressionQ, var_SymbolQ] := IntAux [expression /. var-> VAR]

and then IntAux would see u^2+5v^2+7VAR^2+ 2x^2+4y^2+5z^2. Much easier to match with IntAux[ a_.VAR^2+b_] := a*VAR^3/3+b /; FreeQ[{a,b},VAR]

one glitch .. if we generate some error message from IntAux and friends, those programs will not know the external name of the variable of integration. They will only know VAR. So we have to save the external name. Easy to do at top level..

Int [expression_? ExpressionQ, var_SymbolQ] := (SavedVAR=var ; IntAux [expression /. var-> VAR])

.............. RJF

then change the name of the other rules to private IntAux which would then

not have to verify that the first argument is not, for example a Boolen or the second argument an Integer.

On Tue, Jun 25, 2024 at 8:05 AM Richard Fateman @.***> wrote:

Suggestion for speeding up (Mathematica's version) of Rubi by maybe 30-40%..

It seemed to me that pattern matching should be way faster when there are better "anchors" in the match. Thus if we substitute a particular name, say VAR for the variable of integration, all the rules are simpler. A rule that now looks like

Int [a.*x^2+ ....., xSymbol] ..... could be replaced by Int[a.* VAR^2+.....]

In some simple tests I get about 33% faster, though I am just testing matching speed. (this trick is used in Moses's SIN integration system in Maxima. SIN relies heavily on pattern matching.) RJF

Message ID: @.***>

-- Reply to this email directly or view it on GitHub:

https://github.com/RuleBasedIntegration/Rubi/issues/53#issuecomment-2189775698 You are receiving this because you were mentioned.

Message ID: @.***> --3735943886-1724912013-1719352376=:27136--

— Reply to this email directly, view it on GitHub https://github.com/RuleBasedIntegration/Rubi/issues/53#issuecomment-2190033859, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB5ALCQ63NBSTT5DVUQPNFTZJHRFJAVCNFSM6AAAAABJST2JW2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCOJQGAZTGOBVHE . You are receiving this because you commented.Message ID: @.***>

stblake commented 3 months ago

Pretty much a random stab at a ruleset, I notice that Mathematica can integrate 1/(a+b*x^2)^p wrt x, and its answer, x (a + b x^2)^-p (1 + (b x^2)/a)^p Hypergeometric2F1[1/2, p, 3/ 2, -((b x^2)/a)]

Is fairly small.

Rubi has separate rules for specific values of p.. namely 1/4, 3/4, 1/3, 2/3, 1/6, 5/6. Rubi also works for just plain p and gives the answer shown above.

I got all this from reading in a .nb (notebook) file, which, when read in, has a nice display of what patterns it is looking for.

The related .m file, which is harder to read, says it was automatically generated from the .nb file and should not be edited. But those .m files are the ones I think some of us are using as the ground truth of the rules...

The particular file is

I expect all the variations in the separate rules boil down to simplification of Hypergeometric2F1[1/2,p,3/2,-b/a*x^2]. Using just the one rule is, in effect, mapping back from the table of integrals to the algorithm (actually, rule) that generated the table. Thoughts? Richard

Message ID: @.***>

stblake commented 3 months ago

Albert always wanted every special case to be returned in its best form. Thus, although the Hypergeometric function is "the correct answer" for all values of the parameter, and thus would be the choice of many systems, Albert looks at the special cases separately, for 2 reasons.

  1. Why leave it to the user to take the extra step of asking the system to simplify the Hypergeometric function, when Rubi can do it using its rules? One of the jobs of computer algebra is to bring special cases to the attention of the user, in case the simplification is significant.

  2. Sometimes the most elegant integral expression cannot be obtained by simplifying the general case. The elegant solution may differ by a constant.

     David

On 30/06/2024 7:24 pm, Richard Fateman wrote:

Pretty much a random stab at a ruleset, I notice that Mathematica can integrate 1/(a+b*x^2)^p  wrt x, and its answer, x (a + b x^2)^-p (1 + (b x^2)/a)^p Hypergeometric2F1[1/2, p, 3/   2, -((b x^2)/a)]

Is fairly small.

Rubi has separate rules for specific values of p.. namely 1/4, 3/4, 1/3, 2/3, 1/6, 5/6. Rubi also works for just plain p and gives the answer shown above.

I got all this from reading in a .nb  (notebook) file, which, when read in, has a nice display of what patterns it is looking for.

  The related .m file, which is harder to read, says it was automatically generated from the .nb file and should not be edited.  But those .m files are the ones I think some of us are using as the ground truth of the rules...

The particular file is Rubi

https://github.com/RuleBasedIntegration/Rubi/tree/61e9c18ea248061cd83c67882f7c91a73cef912d

/Rubi

https://github.com/RuleBasedIntegration/Rubi/tree/61e9c18ea248061cd83c67882f7c91a73cef912d/Rubi

/IntegrationRules

https://github.com/RuleBasedIntegration/Rubi/tree/61e9c18ea248061cd83c67882f7c91a73cef912d/Rubi/IntegrationRules

/1 Algebraic functions

https://github.com/RuleBasedIntegration/Rubi/tree/61e9c18ea248061cd83c67882f7c91a73cef912d/Rubi/IntegrationRules/1%20Algebraic%20functions

/1.1 Binomial products

https://github.com/RuleBasedIntegration/Rubi/tree/61e9c18ea248061cd83c67882f7c91a73cef912d/Rubi/IntegrationRules/1%20Algebraic%20functions/1.1%20Binomial%20products

/1.1.2 Quadratic

https://github.com/RuleBasedIntegration/Rubi/tree/61e9c18ea248061cd83c67882f7c91a73cef912d/Rubi/IntegrationRules/1%20Algebraic%20functions/1.1%20Binomial%20products/1.1.2%20Quadratic /

1.1.2.1 (a+b x^2)^p.nb

I expect all the variations in the separate rules boil down to simplification of Hypergeometric2F1[1/2,p,3/2,-b/a*x^2]. Using just the one rule is, in effect, mapping back from the table of integrals to the algorithm (actually, rule) that generated the table. Thoughts? Richard

        Message ID:
        ***@***.***>

--------------AUPe3LIDPSYF60AW02xMgtHr Content-Type: text/html; charset=UTF-8 Content-Transfer-Encoding: 8bit

<!DOCTYPE html>

Albert always wanted every special case to be returned in its best form. Thus, although the Hypergeometric function is "the correct answer" for all values of the parameter, and thus would be the choice of many systems, Albert looks at the special cases separately, for 2 reasons.

1. Why leave it to the user to take the extra step of asking the system to simplify the Hypergeometric function, when Rubi can do it using its rules? One of the jobs of computer algebra is to bring special cases to the attention of the user, in case the simplification is significant.

2. Sometimes the most elegant integral expression *cannot* be obtained by simplifying the general case. The elegant solution may differ by a constant.

 David

On 30/06/2024 7:24 pm, Richard Fateman wrote:
Pretty much a random stab at a ruleset, I notice that Mathematica can integrate 1/(a+b*x^2)^p  wrt x,
and its answer, 
x (a + b x^2)^-p (1 + (b x^2)/a)^p Hypergeometric2F1[1/2, p, 3/
  2, -((b x^2)/a)]

Is fairly small.  

Rubi has separate rules for specific values of p.. namely
1/4, 3/4, 1/3, 2/3, 1/6, 5/6.
Rubi also works for just plain p and gives the answer shown above.

I got all this from reading in a .nb  (notebook) file, which, when read in, has a nice display of what patterns it is looking for.

  The related .m file, which is harder to read, says it was automatically
generated from the .nb file and should not be edited. 
 But those .m files are the ones I think some of us are
using as the ground truth of the rules...

The particular file is 
  • Rubi
  • Rubi
  • IntegrationRules
  • 1 Algebraic functions
  • 1.1 Binomial products
  • 1.1.2 Quadratic
  • 1.1.2.1 (a+b x^2)^p.nb


    I expect all the variations in the separate rules boil down to simplification of Hypergeometric2F1[1/2,p,3/2,-b/a*x^2].
    Using just the one rule is, in effect, mapping back from the table of integrals to the algorithm (actually, rule) that
    generated the table.
    Thoughts?
    Richard






    Message ID: <RuleBasedIntegration/Rubi/issues/53/2181114566@github.com>

    --------------AUPe3LIDPSYF60AW02xMgtHr--

    stblake commented 3 months ago

    Another consolidation. As is well known from the fundamental theorem of algebra, an nth degree polynomial has n roots. Assume then that (polynomial-in-x)^p * (other-polynomials-in-x)^p2 .... is the form of the integrand. Let us find (by factoring) or [given unknown coefficients] hypothesize and name the roots, r1, r2, ... . [ The rules that include conditions like b^2-4ac=0 etc are trying to get around factoring, I think...] Anyway, we get a chance to name the roots and are left with integrating a constant times a product of (x-rk)^pk. In fact, the integrand can be programmed as..

    v[k_]:=Product[(x-r[i])^p[i],{i,1,k}]

    and the task is then to compute Int[v[k],x] for k=1,2,3,4,.... Indeed, for k>=4 neither Mathematica nor Rubi can compute a closed form. For k<=3, they both compute a closed form. In fact, it seems that if we have only the form ans3 for k=3, we can get the answer for k=3 by computing ans3/. p[3]->0.

    So a whole bunch of stuff can probably be consolidated if we say, Oh, it makes integration much simpler if we rewrite your polynomial (s) [to a power] into factored form C(x-r1)^p1.... . The answer is then ... [ in terms of Appel or Hypergeometric or ... simpler functions]

    Now we could offer to simplify stuff in the sense that we know the constant coefficient of the polynomial is the product of the roots, (etc etc).

    This is going further off track from using Rubi as it is, but maybe worth a thought. Note that the matching problem is drastically reduced into producing a canonical factored representation of polynomials and picking off roots and powers.

    RJF

    On Sun, Jun 30, 2024 at 4:24 PM Richard Fateman @.***> wrote:

    Pretty much a random stab at a ruleset, I notice that Mathematica can integrate 1/(a+b*x^2)^p wrt x, and its answer, x (a + b x^2)^-p (1 + (b x^2)/a)^p Hypergeometric2F1[1/2, p, 3/ 2, -((b x^2)/a)]

    Is fairly small.

    Rubi has separate rules for specific values of p.. namely 1/4, 3/4, 1/3, 2/3, 1/6, 5/6. Rubi also works for just plain p and gives the answer shown above.

    I got all this from reading in a .nb (notebook) file, which, when read in, has a nice display of what patterns it is looking for.

    The related .m file, which is harder to read, says it was automatically generated from the .nb file and should not be edited. But those .m files are the ones I think some of us are using as the ground truth of the rules...

    The particular file is

    /1.1.2.1 (a+b x^2)^p.nb

    I expect all the variations in the separate rules boil down to simplification of Hypergeometric2F1[1/2,p,3/2,-b/a*x^2]. Using just the one rule is, in effect, mapping back from the table of integrals to the algorithm (actually, rule) that generated the table. Thoughts? Richard

    Message ID: @.***>