Skip to content

Match WMA for Gamma[1+x] and Product[...] #1395

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rocky merged 7 commits into from
Feb 21, 2025
Merged

Match WMA for Gamma[1+x] and Product[...] #1395

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
88f6d61
Match WMA for Gamma[1+x] and Product[...]
rocky Feb 19, 2025
36f8dba
Add a Gamma test
rocky Feb 19, 2025
daaec23
Go over Product user documentation
rocky Feb 19, 2025
e0057f2
Appease some flaky homegrown style code
rocky Feb 19, 2025
25bfeee
Add Wiki links and update docs for Sum and Product
rocky Feb 19, 2025
eff651f
Merge branch 'master' into WMA-product+gamma-compatibility
rocky Feb 21, 2025
e432f5c
Add add FIXMEs to cover of recent discussion
rocky Feb 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 22 additions & 8 deletions mathics/builtin/arithmetic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -713,21 +713,35 @@ class Product(IterationFunction, SympyFunction):
outermost-to-innermost order.
</dl>

>> Product[k, {k, 1, 10}]
= 3628800
>> 10!
= 3628800
'Product[k, {k, i, n}]' is defined in terms of <url>
:Factorial:
/doc/reference-of-built-in-symbols/special-functions/gamma-and-related-functions/factorial/</url>:

>> Product[k, {k, i, n}]
= n! / (-1 + i)!

When $i$ is 1, we get the factorial function:
>> Product[k, {k, 1, n}]
= n!

Or more succinctly:
>> Product[k, {k, n}]
= n!

Symbolic products involving the factorial are evaluated:
>> Product[k, {k, 3, n}]
= n! / 2

Examples of numeric evaluation using more complex functions:
>> Product[x^k, {k, 2, 20, 2}]
= x ^ 110

>> Product[2 ^ i, {i, 1, n}]
= 2 ^ (n / 2 + n ^ 2 / 2)

>> Product[f[i], {i, 1, 7}]
= f[1] f[2] f[3] f[4] f[5] f[6] f[7]

Symbolic products involving the factorial are evaluated:
>> Product[k, {k, 3, n}]
= n! / 2

Evaluate the $n$-th primorial:
>> primorial[0] = 1;
>> primorial[n_Integer] := Product[Prime[k], {k, 1, n}];
Expand Down
3 changes: 1 addition & 2 deletions mathics/builtin/specialfns/gamma.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -332,7 +332,6 @@ class Gamma(MPMathMultiFunction):

rules = {
"Gamma[z_, x0_, x1_]": "Gamma[z, x0] - Gamma[z, x1]",

@mmatera mmatera Feb 19, 2025

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In WMA, this rule only applies if the three arguments are numeric, and at least one is a Real number

@rocky rocky Feb 19, 2025

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I imagine then that this kind of check would be better done inside an evaluation method.

If you are up for adding, feel free to do. But if not, let's leave for another PR.

@mmatera mmatera Feb 19, 2025

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, let's postpone for another PR.

"Gamma[1 + z_]": "z!",

@mmatera mmatera Feb 19, 2025 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The rule could be
"Gamma[z_Integer]": "(z-1)!"

@rocky rocky Feb 19, 2025 •
edited
Loading

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this cause any observable difference in output? If not, I'd recommend removing the rule unless there is a good reason to add it.

@mmatera mmatera Feb 19, 2025

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We have different implementations for Gamma and Factorial, relying on different mpmath functions. Right now, mpmath function gives the same answer, so both outputs are the same.

@rocky rocky Feb 19, 2025

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great - then let's not complicate here when there is no user-observable difference.

"Gamma[Undefined]": "Undefined",
"Gamma[x_, Undefined]": "Undefined",
"Gamma[Undefined, y_]": "Undefined",
Expand Down Expand Up @@ -448,7 +447,7 @@ class Pochhammer(SympyFunction):

rules = {
"Pochhammer[0, 0]": "1",
"Pochhammer[a_, n_]": "Gamma[a + n] / Gamma[a]",
"Pochhammer[a_, n_]": "Factorial[a + n - 1] / Factorial[a - 1]",

@mmatera mmatera Feb 19, 2025 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is what Pochhammer should do:

In[1]:= Pochhammer[a,b]                                                         

Out[1]= Pochhammer[a, b]

In[2]:= Pochhammer[a,12.4]                                                      

Out[2]= Pochhammer[a, 12.4]

In[3]:= Pochhammer[a,12.]                                                       

Out[3]= a (1 + a) (2 + a) (3 + a) (4 + a) (5 + a) (6 + a) (7 + a) (8 + a) 
 
>    (9 + a) (10 + a) (11 + a)

In[4]:= Pochhammer[12,b]                                                        

Out[4]= Pochhammer[12, b]

In[5]:= Pochhammer[12.1,b]                                                      

Out[5]= Pochhammer[12.1, b]

Notice that n must be an Integer to apply this rule.

Suggested change
"Pochhammer[a_, n_]": "Factorial[a + n - 1] / Factorial[a - 1]",
"Pochhammer[a_, n_Integer]": "Factorial[a + n - 1] / Factorial[a - 1]",

or maybe n_?NumerQ

@rocky rocky Feb 19, 2025

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to the examples cited, NumberQ only works when the number can be exactly converted to an Integer.

@mmatera mmatera Feb 19, 2025

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed. In any case, this is the point when we need to decide if the complexity of a rule that matches the exact behavior worth. The exact behavior would be obtained with something like

        "Pochhammer[a_, /;(IntegerQ[n]|| NumberQ[n] && Round[n]==n)]": "Factorial[a + n - 1] / Factorial[a - 1]",

@rocky rocky Feb 19, 2025 •
edited
Loading

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IntegerQ[n] can be dropped following the principle: simple as possible.... But even if we change to Pochhammer[a_, n_(NumberQ[n] && Round[n]==n)], Product[k, {k, i, n}] will not transform as WMA does, in particular: n! / (-1+i)!. So WMA may be using some other mechanism to implement Product.

But let's back up a bit. We got into this as a result of a real problem reported by @aravindh-krishnamoorthy, in his day-to-day use of Mathics3. His PR broke something else. If this PR addresses both his problem and does not break what is already there, then let's go with this.

Note that this PR adds a test to ensure we do not rewrite Gamma as Factorial in the kind of way that @aravindh-krishnamoorthy was not expecting.

In the future, if someone has a problem as a result of Pochhammer getting expanded to factorial, then we can address that issue at that time while maintaining satisfaction of all the other checks in place.

In other words: complicate when there is a tangible need to complicate. Again this is also encompassed by the principle: Simple as possible, but no simpler.

"Derivative[1,0][Pochhammer]": "(Pochhammer[#1, #2]*(-PolyGamma[0, #1] + PolyGamma[0, #1 + #2]))&",
"Derivative[0,1][Pochhammer]": "(Pochhammer[#1, #2]*PolyGamma[0, #1 + #2])&",
}
Expand Down
1 change: 1 addition & 0 deletions test/builtin/specialfns/test_gamma.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
"Overflow",
),
("Gamma[1., 2.]", None, "Gamma[1., 2.]", "needs mpmath for lowergamma"),
("Gamma[1 + x]", None, "Gamma[1 + x]", "Gamma should not expand to factorial"),
],
)
def test_private_doctests_gamma(str_expr, msgs, str_expected, fail_msg):
Expand Down