Menu
▾
▴
Re: [Maxima-discuss] Simplifying gamma functions?
|
From: Barton W. <wi...@un...> - 2022-07-21 18:27:35
|
Does SymPy have an rtest-like file for testing the code that simplifies products of gamma expressions?
--Barton
From: Barton Willis via Maxima-discuss<mailto:max...@li...>
Sent: Thursday, July 21, 2022 7:24 AM
To: Oscar Benjamin<mailto:osc...@gm...>; Raymond Toy<mailto:toy...@gm...>
Cc: <max...@li...><mailto:max...@li...>
Subject: Re: [Maxima-discuss] Simplifying gamma functions?
Non-NU Email
I’m humiliated:
gamma_product_simp(e):=block([a : map('first, gatherargs(e,'gamma))],
for ak in a do (
if not featurep(ak,'integer) then e : ratsubst(%pi/sin(%pi*ak),gamma(ak)*gamma(1-ak),e)),
a : map('first, gatherargs(e,'gamma)),
for ak in a do (
e : ratsubst(sqrt(%pi) * 2^(1-2*ak) * gamma(2*ak), gamma(ak)*gamma((2*ak+1)/2),e)),
e)$
From: Barton Willis via Maxima-discuss<mailto:max...@li...>
Sent: Thursday, July 21, 2022 6:59 AM
To: Oscar Benjamin<mailto:osc...@gm...>; Raymond Toy<mailto:toy...@gm...>
Cc: <max...@li...><mailto:max...@li...>
Subject: Re: [Maxima-discuss] Simplifying gamma functions?
Non-NU Email
An alert reader mentioned that my code for the gamma(x) * gamma(x+1/2) identity is missing a factor. A less broken version:
load("opsubst")$
gamma_product_simp(e):=block([a : map('first, gatherargs(e,'gamma))],
for ak in a do (
if not featurep(ak,'integer) then e : ratsubst(%pi/sin(%pi*ak),gamma(ak)*gamma(1-ak),e)),
a : map('first, gatherargs(e,'gamma)),
for ak in a do (
e : ratsubst(sqrt(%pi) * 2^(1-2*ak) * gamma(2*ak), gamma(x)*gamma((2*x+1)/2),e)),
e)$
But the function gamma_product_simp isn’t idempotent. Example:
(%i32) gamma_product_simp(gamma(x)*gamma(x+1/2) * gamma(1/3));
(%o32) 2^(1/3)*sqrt(%pi)*gamma(1/3)*gamma(2/3)
(%i33) gamma_product_simp(%);
(%o33) (2^(4/3)*%pi^(3/2))/sqrt(3)
When I get a chance, I’ll fix that. Likely the infinite recursion prone fullratsubst isn’t the cure—rather the cure is to loop over pairs of the arguments of the gamma function and to take notice when a sum of the pairs is 1 or the difference is ½. When such a sum or difference is found, declare a do-over.
I nominate gatherargs, or some variation of it, for membership in the Maxima source. The need to double quote opsubst to load the package has caused frustration to users, including me.
As I recall, minfactorial isn’t idempotent either.
--Barton
From: Barton Willis via Maxima-discuss<mailto:max...@li...>
Sent: Wednesday, July 20, 2022 10:06 PM
To: Oscar Benjamin<mailto:osc...@gm...>; Raymond Toy<mailto:toy...@gm...>
Cc: <max...@li...><mailto:max...@li...>
Subject: Re: [Maxima-discuss] Simplifying gamma functions?
Non-NU Email
First waffle:
(%i122) load("opsubst")$
(%i123) gamma_product_simp(e):=block([a : map('first, gatherargs(e,'gamma))],
for ak in a do (
if not featurep(ak,'integer) then e : ratsubst(%pi/sin(%pi*ak),gamma(ak)*gamma(1-ak),e)),
a : map('first, gatherargs(e,'gamma)),
for ak in a do (
e : ratsubst(gamma(2*ak), gamma(x)*gamma((2*x+1)/2),e)),
e)$
(%i124) gamma_product_simp(gamma(2/3)*gamma(1/3));
(%o124) (2*%pi)/sqrt(3)
(%i125) gamma_product_simp(gamma(x)*gamma(x+1/2));
(%o125) gamma(2*x)
Many things to ponder—the order of the substitutions, ratsubst or fullratsubst, featurep vs askinteger vs … the Legendre triplication formula, and what else?
From: Oscar Benjamin<mailto:osc...@gm...>
Sent: Wednesday, July 20, 2022 4:29 PM
To: Raymond Toy<mailto:toy...@gm...>
Cc: <max...@li...><mailto:max...@li...>
Subject: Re: [Maxima-discuss] Simplifying gamma functions?
Non-NU Email
On Wed, 20 Jul 2022 at 06:02, Raymond Toy <toy...@gm...> wrote:
>
> Is there some builtin function to simplify a product of gamma functions?
>
> I can't simplify gamma(1/4)*gamma(3/4).
>
> For example gamma(1-z)*gamma(z) = %pi/sin(%pi*z). In this case if z = 1/4, we can see that the product is %pi/sin(%pi/4) = sqrt(2)*%pi.
>
> Or maybe note that beta(1/4,3/4) = gamma(1/4)*gamma(3/4)/gamma(1), so gamma(1/4)*gamma(3/4) = beta(1/4,3/4)*gamma(1). And maxima knows that beta(1/4,3/4) = sqrt(2)*%pi.
>
> The beta form can be useful for gamma(x)*gamma(y) where x+y is an integer or half integer.
>
> The reflection formula is useful if y = 1-x.
>
> The duplication formula gamma(z)*gamma(z+1/2) = 2^(1-2*z)*sqrt(%pi)*gamma(2*z) is also useful here when z = 1/4.
Not sure if this is useful but there is a function gammasimp in SymPy
that handles the cases described above:
>>> gammasimp(gamma(z)*gamma(1-z))
pi/sin(pi*z)
>>> gammasimp(gamma(S(1)/4)*gamma(S(3)/4))
sqrt(2)*pi
>>> gammasimp(gamma(z)*gamma(z + S(1)/2))
2*sqrt(2)*2**(-2*z - 1/2)*sqrt(pi)*gamma(2*z)
The code for it is here:
https://urldefense.com/v3/__https://github.com/sympy/sympy/blob/master/sympy/simplify/gammasimp.py__;!!PvXuogZ4sRB2p-tU!FBC3ZGVwsQg5hM3uLxEVQOApRjziHgp4j8F52mOELMyuw_tDI7pTUZq_PMtPfE-Pbo7zbW_yZbt7gQc5NJZYbVUgDA$<https://urldefense.com/v3/__https:/github.com/sympy/sympy/blob/master/sympy/simplify/gammasimp.py__;!!PvXuogZ4sRB2p-tU!FBC3ZGVwsQg5hM3uLxEVQOApRjziHgp4j8F52mOELMyuw_tDI7pTUZq_PMtPfE-Pbo7zbW_yZbt7gQc5NJZYbVUgDA$>
I haven't studied it in detail.
--
Oscar
_______________________________________________
Maxima-discuss mailing list
Max...@li...
https://urldefense.com/v3/__https://lists.sourceforge.net/lists/listinfo/maxima-discuss__;!!PvXuogZ4sRB2p-tU!FBC3ZGVwsQg5hM3uLxEVQOApRjziHgp4j8F52mOELMyuw_tDI7pTUZq_PMtPfE-Pbo7zbW_yZbt7gQc5NJZ6QFyxSQ$<https://urldefense.com/v3/__https:/lists.sourceforge.net/lists/listinfo/maxima-discuss__;!!PvXuogZ4sRB2p-tU!FBC3ZGVwsQg5hM3uLxEVQOApRjziHgp4j8F52mOELMyuw_tDI7pTUZq_PMtPfE-Pbo7zbW_yZbt7gQc5NJZ6QFyxSQ$>
|
