Menu
▾
▴
Re: [Maxima-discuss] Simplifying gamma functions?
|
From: Oscar B. <osc...@gm...> - 2022-07-21 23:37:25
|
I'm not sure what an r-test file is. Is that a file for a specific kind of unit test runner? The explicit unit tests for gammasimp in SymPy are here: https://github.com/sympy/sympy/blob/master/sympy/simplify/tests/test_gammasimp.py The gammasimp function is also used as part of other functions such as combsimp which is itself used as part of SymPy's general simplify function. Tests for some other simplification functions can be found in the same directory. Oscar On Thu, 21 Jul 2022 at 19:27, Barton Willis <wi...@un...> wrote: > Does SymPy have an rtest-like file for testing the code that simplifies > products of gamma expressions? > > > > --Barton > > > > *From: *Barton Willis via Maxima-discuss > <max...@li...> > *Sent: *Thursday, July 21, 2022 7:24 AM > *To: *Oscar Benjamin <osc...@gm...>; Raymond Toy > <toy...@gm...> > *Cc: *<max...@li...> > <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 > <max...@li...> > *Sent: *Thursday, July 21, 2022 6:59 AM > *To: *Oscar Benjamin <osc...@gm...>; Raymond Toy > <toy...@gm...> > *Cc: *<max...@li...> > <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 > <max...@li...> > *Sent: *Wednesday, July 20, 2022 10:06 PM > *To: *Oscar Benjamin <osc...@gm...>; Raymond Toy > <toy...@gm...> > *Cc: *<max...@li...> > <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 <osc...@gm...> > *Sent: *Wednesday, July 20, 2022 4:29 PM > *To: *Raymond Toy <toy...@gm...> > *Cc: *<max...@li...> > <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$> > > > > > > > > > |
