share/gamma_simp/rtest_gamma_simp.mac

   1
   2 /* 22 July 2022
   3
   4 Tests adapted from SymPy; see https://github.com/sympy/sympy/blob/master/sympy/simplify/tests/test_gammasimp.py
   5
   6 All rights reserved.
   7
   8 Redistribution and use in source and binary forms, with or without
   9 modification, are permitted provided that the following conditions are met:
  10
  11   a. Redistributions of source code must retain the above copyright notice,
  12      this list of conditions and the following disclaimer.
  13   b. Redistributions in binary form must reproduce the above copyright
  14      notice, this list of conditions and the following disclaimer in the
  15      documentation and/or other materials provided with the distribution.
  16   c. Neither the name of SymPy nor the names of its contributors
  17      may be used to endorse or promote products derived from this software
  18      without specific prior written permission.
  19
  20
  21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  22 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24 ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
  25 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  27 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  28 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30 OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
  31 DAMAGE.
  32
  33 --------------------------------------------------------------------------------
  34
  35 Patches that were taken from the Diofant project (https://github.com/diofant/diofant)
  36 are licensed as:
  37
  38 Copyright (c) 2006-2018 SymPy Development Team,
  39               2013-2022 Sergey B Kirpichev
  40
  41 All rights reserved.
  42
  43 Redistribution and use in source and binary forms, with or without
  44 modification, are permitted provided that the following conditions are met:
  45
  46   a. Redistributions of source code must retain the above copyright notice,
  47      this list of conditions and the following disclaimer.
  48   b. Redistributions in binary form must reproduce the above copyright
  49      notice, this list of conditions and the following disclaimer in the
  50      documentation and/or other materials provided with the distribution.
  51   c. Neither the name of Diofant or SymPy nor the names of its contributors
  52      may be used to endorse or promote products derived from this software
  53      without specific prior written permission.
  54
  55
  56 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  57 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  58 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  59 ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
  60 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  61 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  62 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  63 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  64 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  65 OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
  66 DAMAGE.
  67
  68 */
  69 /* Let's attempt to limit the number of failures that are due to the value
  70    of various option variables. */
  71 (load("gamma_simp"),%piargs : true, ratfac : true, domain : real,m1pbranch : false,
  72 radsubstflag : false, reset(float_approx_equal_tolerance), 0);
  73 0$
  74
  75 gamma_simp(gamma(x));
  76 gamma(x)$
  77
  78 gamma_simp(gamma(x + 1)/x);
  79 gamma(x)$
  80
  81 gamma_simp(gamma(x)/(x - 1));
  82 gamma(x - 1)$
  83
  84 gamma_simp(x*gamma(x));
  85 gamma(x + 1)$
  86
  87 gamma_simp((x + 1)*gamma(x + 1));
  88 gamma(x + 2)$
  89
  90 gamma_simp(gamma(x + y)*(x + y));
  91 gamma(x + y + 1)$
  92
  93 gamma_simp(x/gamma(x + 1));
  94 1/gamma(x)$
  95
  96 gamma_simp((x + 1)^2/gamma(x + 2));
  97 (x + 1)/gamma(x + 1)$
  98
  99 /* Another answer is x*(x+2)*gamma(x). Let's be happy with (x+2)*gamma(x+1) */
 100 gamma_simp(x*gamma(x) + gamma(x + 3)/(x + 2));
 101 (x+2)*gamma(x+1)$
 102
 103 gamma_simp(gamma(2*x)*x);
 104 gamma(2*x + 1)/2$
 105
 106 gamma_simp(gamma(2*x)/(x - 1/2));
 107 2*gamma(2*x - 1)$
 108
 109 gamma_simp(gamma(x)*gamma(1 - x));
 110 %pi/sin(%pi*x)$
 111
 112 gamma_simp(gamma(x)*gamma(-x));
 113 -(%pi/(x*sin(%pi*x)))$
 114
 115 gamma_simp(1/gamma(x + 3)/gamma(1 - x));
 116 sin(%pi*x)/(%pi*x^3+3*%pi*x^2+2*%pi*x)$
 117
 118 gamma_simp((n + 2)!);
 119 (n+2)!$
 120
 121 gamma_simp(makegamma((n+2)!));
 122 gamma(n+3)$
 123
 124 gamma_simp(binomial(n, k));
 125 (binomial(n,k))$
 126
 127 gamma_simp(makegamma(binomial(n, k)));
 128 gamma(n+1)/(gamma(k+1)*gamma(n-k+1))$
 129
 130 gamma_simp(gamma(x)*gamma(x + 1/2)*gamma(y)/gamma(x + y));
 131 2^(-2*x + 1)*sqrt(%pi)*gamma(2*x)*gamma(y)/gamma(x + y)$
 132
 133 gamma_simp(1/gamma(x)/gamma(x - 1/3)/gamma(x + 1/3));
 134 3^(3*x - 3/2)/(2*%pi*gamma(3*x - 1))$
 135
 136 gamma_simp(gamma(1/2 + x/2)*gamma(1 + x/2)/gamma(1 + x)/sqrt(%pi)*2**x);
 137 1$
 138
 139 gamma_simp(gamma(-1/4)*gamma(-3/4));
 140 16*sqrt(2)*%pi/3$
 141
 142 gamma_simp(gamma(2*x)/gamma(x));
 143 (2^(2*x-1)*gamma((2*x+1)/2))/sqrt(%pi)$
 144
 145 (e : (-gamma(k)*gamma(k + 2) + gamma(k + 1)^2)/gamma(k)^2,0);
 146 0$
 147
 148 gamma_simp(e);
 149 -k$
 150
 151 gamma_simp(1/e);
 152 -1/k$
 153
 154 (e : (gamma(x) + gamma(x + 1))/gamma(x),0);
 155 0$
 156
 157 gamma_simp(e);
 158 x + 1$
 159
 160 gamma_simp(1/e);
 161 1/(x + 1)$
 162
 163 (e : (gamma(x) + gamma(x + 2))*(gamma(x - 1) + gamma(x))/gamma(x),0);
 164 0$
 165
 166 /*Apparently, SymPy simplifies this to (x^2 + x + 1)*gamma(x + 1)/(x - 1).
 167 But Maxima gives gamma(x-1)*x*(x^2+x+1). I think Maxima's answer is better.
 168 Guess: the difference is due to a fluke in the order of the matching.*/
 169 gamma_simp(e);
 170 gamma(x-1)*x*(x^2+x+1)$
 171
 172 (e : (-gamma(k)*gamma(k + 2) + gamma(k + 1)^2)/gamma(k)^2,0);
 173 0$
 174
 175 gamma_simp(e^2);
 176 k^2$
 177
 178 gamma_simp(e^2/gamma(k + 1));
 179 k/gamma(k)$
 180
 181 gamma_simp(gamma(2*k)/gamma(k)*gamma(k + 1/2+1/3)*gamma(k + 1/2+2/3));
 182 sqrt(%pi)*3^(-3*k-1)*2^(2*k)*gamma((6*k+3)/2)$
 183
 184 gamma_simp(makegamma((x + 1)*factorial(x)/gamma(y)));
 185 gamma(x + 2)/gamma(y)$
 186
 187 gamma_simp(pochhammer(x + n, k)*binomial(n, k));
 188 pochhammer(x + n, k)*binomial(n, k)$
 189
 190 /* SymPy apparently gives a split rule splitting on the sign of n+x. */
 191 gamma_simp(makegamma(pochhammer(x + n, k)*binomial(n, k)));
 192 gamma(n + 1)*gamma(k + n + x)/(gamma(k + 1)*gamma(n + x)*gamma(-k + n + 1))$
 193
 194 /* I'm lost here:
 195     A, B = symbols('A B', commutative=False)
 196     assert gamma_simp(e*B*A) == gamma_simp(e)*B*A
 197 */
 198
 199 /* The setting of algebraic alters the solution a bit */
 200 block([algebraic : false], gamma_simp(gamma(2*k)/gamma(k)*gamma(-k - 1/2)));
 201 -((%pi*2^(2*k))/((2*sqrt(%pi)*k+sqrt(%pi))*sin((2*%pi*k+%pi)/2)))$
 202
 203 gamma_simp(gamma(k)*gamma(k + 1/3)*gamma(k + 2/3)/gamma(k*3/2));
 204 sqrt(%pi)*3^(1/2-3*k)*2^(3*k)*gamma((3*k+1)/2)$
 205
 206 gamma_simp(gamma(1/4)/gamma(5/4));
 207 4$
 208
 209 gamma_simp(binomial(n + 2, k + 1/2));
 210 binomial(n + 2, k + 1/2)$
 211
 212 gamma_simp(makegamma(binomial(n + 2, k + 1/2)));
 213 gamma(n+3)/(gamma((2*k+3)/2)*gamma((2*n-2*k+5)/2))$
 214
 215 gamma_simp(binomial(n + 2, k + 2.0));
 216 binomial(n + 2, k + 2.0)$
 217
 218 /* The function ratsubst does not respect the value of keepfloat--this
 219    is its documented behavior. The gamma_simp code is a heavy user of
 220    ratsubst--fixing this bug likely isn't worthwhile. */
 221 gamma_simp(makegamma(binomial(n + 2, k + 2.0)));
 222 gamma(n + 3)/(gamma(k + 3.0)*gamma(-k + n + 1))$
 223
 224 gamma_simp(binomial(0, x));
 225 binomial(0, x)$
 226
 227 gamma_simp(makegamma(binomial(0, x)));
 228 sin(%pi*x)/(%pi*x)$
 229
 230 (e : gamma(n + 1/3)*gamma(n + 2/3),0);
 231 0$
 232
 233 gamma_simp(e);
 234 gamma((3*n+1)/3)*gamma((3*n+2)/3)$
 235
 236 gamma_simp(gamma(4*n + 1/2)/gamma(2*n - 3/4));
 237 ((8*n-3)*2^(4*n-5/2)*gamma((8*n+3)/4))/sqrt(%pi)$
 238
 239 (e : gamma(exp(i)),0);
 240 0$
 241
 242 gamma_simp(e);
 243 gamma(exp(i))$
 244
 245 (e : gamma(m + 3),0);
 246 0$
 247
 248 gamma_simp(e);
 249  gamma(m + 3)$
 250
 251 (e : gamma(m + 1)/(gamma(i + 1)*gamma(-i + m + 1)),0);
 252 0$
 253
 254 gamma_simp(e);
 255 gamma(m + 1)/(gamma(i + 1)*gamma(-i + m + 1))$
 256
 257 (declare(p,integer), assume(p>0),0);
 258 0$
 259
 260 gamma_simp(gamma(-p + 4));
 261 gamma(-p + 4)$
 262
 263 /* End of tests from SymPy*/
 264
 265 gamma_simp(2/3);
 266 2/3$
 267
 268 gamma_simp(-2/3);
 269 -2/3$
 270
 271 gamma_simp(1.3);
 272 1.3$
 273
 274 gamma_simp(1.3b9);
 275 1.3b9$
 276
 277 gamma_simp([%pi,%e,%phi,minf,-inf,inf]);
 278 [%pi,%e,%phi,minf,-inf,inf]$
 279
 280 gamma_simp(x+1/2);
 281 x+1/2$
 282
 283 gamma_simp(sin(a+y)-z);
 284 sin(a+y)-z$
 285
 286 (X : gamma_simp(x*gamma(x)) - x*gamma(x),0);
 287 0$
 288
 289 makelist(subst(x=k/5,X),k,1,16);
 290 [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]$
 291
 292 gamma_simp(%pi*x*gamma(x));
 293 %pi*gamma(x+1)$
 294
 295 gamma_simp((a+1/2)*gamma(a+1/2));
 296 gamma((2*a+3)/2)$
 297
 298 gamma_simp((2*a+1)*gamma(a+1/2))/2;
 299 gamma((2*a+3)/2)$
 300
 301 gamma_simp(7*x*gamma(x));
 302 7*gamma(x+1)$
 303
 304 gamma_simp(32*x^2*gamma(x));
 305 32*x*gamma(x+1)$
 306
 307 gamma_simp(32*x^2*gamma(x) - 32*x*gamma(x+1));
 308 0$
 309
 310 gamma_simp(x^2*gamma(x)^2/8 - gamma(x+1)^2/8);
 311 0$
 312
 313 gamma_simp(gamma(2*x)/gamma(x)- (2^(2*x -1) * gamma(x+1/2) / sqrt(%pi)));
 314 0$
 315
 316 gamma_simp(sqrt(%pi)*2^(2*x+2)*gamma(x)*gamma(x+1)*gamma((2*x+3)/2)-2*%pi*gamma(x)*gamma(2*x+2));
 317 0$
 318
 319 gamma_simp(sqrt(%pi)*2^(2*x)*gamma(x)*gamma(x+1)*gamma((2*x+1)/2)-2*%pi*gamma(2*x)*gamma(x+1));
 320 0$
 321
 322 /*See: #2118 scanmap(minfactorial,a!) infinite loop. This bug has been fixed.  */
 323 scanmap(factorial_simp,a!) ;
 324 a!$
 325
 326 /*See: #920 minfactorial doesn't look inside "!". When factorial_expand is true,
 327 the factorial code simplifies the input to n!. OK, but minfactorial cannot
 328 do this simplify this to n! when factorial_expand is false.  */
 329 block([factorial_expand : false],factorial_simp((n!/(n-1)!)!));
 330 n!$
 331
 332 /* See bug #1. */
 333 gamma_simp(x/gamma(x+1));
 334 1/gamma(x)$
 335
 336 gamma_simp(%pi*x*(x+1/2)*(x+3/2)*gamma((2*x+1)/2));
 337 %pi*x*gamma((2*x+5)/2)$
 338
 339 gamma_simp(%pi*x*(x+1/2)*(x+3/2)*gamma(x+1/2));
 340 %pi*x*gamma((2*x+5)/2)$
 341
 342 gamma_simp(%pi*x*(x+1/28)*(x+3/2)*gamma(x+1/2));
 343 (%pi*x*(2*x+3)*(28*x+1)*gamma((2*x+1)/2))/56$
 344
 345 gamma_simp(x/gamma(x+1) - 1/gamma(x));
 346 0$
 347
 348 gamma_simp(gamma(2*x+1)/2 - x*gamma(2*x));
 349 0$
 350
 351 gamma_simp((2^(2*x-1)*gamma(x+1/2))/sqrt(%pi) - gamma(2*x)/gamma(x));
 352 0$
 353
 354 gamma_simp((gamma(x)*gamma(2*(x+1)))/(gamma(2*x)*gamma(x+1)));
 355 2*(2*x+1)$
 356
 357 gamma_simp(gamma(x^2)/(x^2-1));
 358 gamma(x^2-1)$
 359
 360 gamma_simp(gamma(2*x)/(2*x-1));
 361 gamma(2*x-1)$
 362
 363 gamma_simp(gamma(2*x)/(x-1/2));
 364 2*gamma(2*x-1)$
 365
 366 gamma_simp(gamma(1/2 + %i)*gamma(1/2 - %i));
 367 %pi/cosh(%pi)$
 368
 369 gamma_simp(gamma(1/2 + %i*y)*gamma(1/2 - %i*y));
 370 %pi/cosh(%pi*y)$
 371
 372 gamma_simp(gamma(%i*y)*gamma(- %i*y));
 373 %pi/(y*sinh(%pi*y))$
 374
 375 (f(z,b,n) := product(gamma(z + (b+k)/n),k,0,n-1),0);
 376 0$
 377
 378 gamma_simp(f(x,-1,2));
 379 sqrt(%pi)*2^(2-2*x)*gamma(2*x-1)$
 380
 381 gamma_simp(f(x,1/8,2));
 382 sqrt(%pi)*2^(7/8-2*x)*gamma((16*x+1)/8)$
 383
 384 gamma_simp(f(x,1/8,12));
 385 %pi^(11/2)*2^(25/4-24*x)*3^(3/8-12*x)*gamma((96*x+1)/8)$
 386
 387 float_approx_equal(subst(x=2.3, gamma_simp(f(x,1/8,2))) , subst(x=2.3, f(x,1/8,2)));
 388 true$
 389
 390 float_approx_equal(subst(x=2.3, gamma_simp(f(x,1/8,4))) , subst(x=2.3, f(x,1/8,4)));
 391 true$
 392
 393 float_approx_equal(subst(x=1.3, gamma_simp(f(x,2,6))) , subst(x=1.3, f(x,2,6)));
 394 true$
 395
 396 /* Bug #2:  misses simplification for gamma(1/2-x) * gamma(1/2+x) #2 */
 397 gamma_simp(gamma(1/2-x) * gamma(1/2+x));
 398 -(%pi/sin((2*%pi*x-%pi)/2))$
 399
 400 gamma_simp(1/(x*gamma(x)));
 401 1/gamma(x+1)$
 402
 403 /* Bug #5:  gamma_simp(gamma(%i * a)*gamma(-%i*a)) #5 */
 404 gamma_simp(gamma(%i * a)*gamma(-%i*a));
 405 %pi/(a*sinh(%pi*a))$
 406
 407 /* See https://mathematica.stackexchange.com/questions/124546/to-what-extent-can-one-simplify-these-products-of-gamma-functions
 408 The denominator of the result vanishes when n=-1 and when n=-3/2. But the input
 409 isn't defined for either n=-1 or n=-3/2, so I say this is OK.*/
 410
 411 factorial_simp((1/10)*(23 + 10*n)*(1 + n)!*(2 + n)!*(1 + 2*n)!*(3 + 2*n)!*((1/10)*(7 + 10*n))!*
 412   ((1/10)*(9 + 10*n))!*((1/10)*(11 + 10*n))!*((1/10)*(13 + 10*n))!*((1/10)*(17 + 10*n))!*
 413   ((1/10)*(19 + 10*n))!*((1/10)*(21 + 10*n))!);
 414   %pi*(n+1)^3*(n+2)*(10*n+23)*5^(-5*n-9)*2^(4*n+5)*n!^4*((2*n+3)/2-1)!*
 415   ((10*n+17)/2-1)!*((10*n+27)/10-1)!*((10*n+29)/10-1)!*((10*n+31)/10-1)!$
 416
 417 /* This test takes a long time. But it did reveal a bug in gamma_simp (that has
 418    been fixed). So let's keep it.*/
 419   (xxx : gamma(1/2 + k)*gamma(7/10 + k)*gamma(9/10 + k)*gamma(1 + k)^6*
 420   gamma(11/10 + k)*gamma(13/10 + k)*gamma(3/2 + k)*gamma(17/10 + k)*
 421   gamma(19/10 + k)*gamma(21/10 + k)*gamma(23/10 + k)*gamma(5/2 + k)*
 422   gamma(13/5 + k)*gamma(27/10 + k)*gamma(29/10 + k)*gamma(3 + k)*
 423   gamma(31/10 + k)*gamma(33/10 + k)*gamma(7/2 + k)*gamma(18/5 + k)*
 424   gamma(37/10 + k)*gamma(39/10 + k)*gamma(41/10 + k)*gamma(43/10 + k)*
 425   gamma(9/2 + k)*gamma(23/5 + k)*gamma(47/10 + k)*gamma(49/10 + k)*
 426   gamma(51/10 + k)*gamma(53/10 + k)*gamma(28/5 + k), yyy : gamma_simp(xxx),0);
 427   0$
 428
 429 factor(yyy);
 430 (%pi^10*(k+1)*(k+2)*(2*k+1)^4*(2*k+3)^3*(2*k+5)^2*(2*k+7)*
 431 (5*k+13)^3*(5*k+18)^2*(5*k+23)*(10*k+7)^4*(10*k+9)^4*(10*k+11)^4*
 432 (10*k+13)^4*(10*k+17)^3*(10*k+19)^3*(10*k+21)^3*(10*k+23)^3*(10*k+27)^2*
 433 (10*k+29)^2*(10*k+31)^2*(10*k+33)^2*(10*k+37)*(10*k+39)*(10*k+41)*(10*k+43)*
 434 5^(-25*k-56)*gamma(k+1)^7*gamma((5*(2*k+1))/2)^5*gamma((5*k+13)/5)^4)/
 435 (1099511627776)$
 436
 437 trigrat(gamma_simp(subst(k=1,xxx/yyy)));
 438 1$
 439
 440 gamma_simp(subst(k=3/8,xxx/yyy));
 441 1$
 442
 443 block([float_approx_equal_tolerance : float(1/2^48)],
 444    float_approx_equal(subst(k=3/8,xxx), subst(k=3/8,yyy)));
 445 true$
 446
 447
 448 /*See: https://math.stackexchange.com/questions/237775/how-is-this-expression-with-the-gamma-function-simplify */
 449 gamma_simp(gamma(a+b)*gamma(a+1)/(gamma(a)*gamma(a+b+1)));
 450 a/(b+a)$
 451
 452 /*See: https://mathematica.stackexchange.com/questions/164504/simplify-fraction-of-gamma-functions  */
 453 gamma_simp(sqrt(gamma(2*n)/gamma(2*n+3)));
 454 1/(2*sqrt(n*(2*n^2+3*n+1)))$
 455
 456 /*See: https://resources.wolframcloud.com/FunctionRepository/resources/GammaSimplify  */
 457 gamma_simp(gamma(z-1)+gamma(z));
 458 gamma(z-1)*z$
 459
 460 /* From rtest_gamma */
 461 minfactorial(factorial(factorial(n)/factorial(n-1)));
 462 factorial(factorial(n)/factorial(n-1));
 463
 464 /* From rtest10 */
 465 minfactorial(n!/(n+1)!);
 466 1/(n+1)$
 467
 468 /* pochhammer tests */
 469
 470 (cntx : newcontext(), activate(cntx), declare(m,integer),declare(k,integer), declare(n,integer),declare(N,integer), 0);
 471 0$
 472
 473 gamma_simp(makegamma(pochhammer(a,n) = (a+n-1)* pochhammer(a-1,n)/(a-1)));
 474 gamma(n+a)/gamma(a)=gamma(n+a)/gamma(a)$
 475
 476 gamma_simp(makegamma(pochhammer(a,n) = (a/(a+n))* pochhammer(a+1,n)));
 477 gamma(n+a)/gamma(a)=gamma(n+a)/gamma(a)$
 478
 479 gamma_simp(makegamma(pochhammer(a,n) = (a + n-1)* pochhammer(a,n-1)));
 480 gamma(n+a)/gamma(a)=gamma(n+a)/gamma(a)$
 481
 482 gamma_simp(makegamma(pochhammer(a,n) = pochhammer(a,n+1)/(a+n)));
 483 gamma(n+a)/gamma(a)=gamma(n+a)/gamma(a)$
 484
 485 /* DLMF: http://dlmf.nist.gov/5.2.E7 . A better answer is
 486    pochhammer(a,2*n)=pochhammer(a,2*n)$*/
 487 gamma_simp(makegamma(pochhammer(a,2*n) = 2^(2*n) * pochhammer(a/2,n) * pochhammer((a+1)/2,n)));
 488 gamma(2*n+a)/gamma(a)=gamma(2*n+a)/gamma(a)$
 489
 490 /* DLMF: http://dlmf.nist.gov/5.2.E8 . A better answer is
 491  pochhammer(a,2*n+1)=pochhammer(a,2*n+1)$ */
 492 gamma_simp(makegamma(pochhammer(a,2*n+1) = 2^(2*n+1) * pochhammer(a/2,n+1) * pochhammer((a+1)/2,n)));
 493 gamma(2*n+a+1)/gamma(a)=gamma(2*n+a+1)/gamma(a)$
 494
 495 /* See https://math.stackexchange.com/questions/656505/an-identity-involving-the-pochhammer-symbol
 496 Also see bug: gamma_simp not idempotent #6 */
 497
 498 (xxx : (6*x)!/(3*x)! = ( 1728^x * pochhammer(1/6,x)*pochhammer(1/2,x) * pochhammer(5/6,x)),0);
 499 0$
 500
 501 gamma_simp(makegamma(xxx));
 502 (2^(6*x)*gamma((6*x+1)/2))/sqrt(%pi)=(1728^x*gamma((6*x+1)/2))/(sqrt(%pi)*3^(3*x))$
 503
 504 is(gamma_simp(makegamma(xxx)) = gamma_simp(gamma_simp(makegamma(xxx))));
 505 true$
 506
 507 gamma_simp(makegamma((6*n)!/(3*n)! =( 1728^n * pochhammer(1/6,n)*pochhammer(1/2,n)
 508 * pochhammer(5/6,n))));
 509 (2^(6*n)*gamma((6*n+1)/2))/sqrt(%pi)=(1728^n*gamma((6*n+1)/2))/(sqrt(%pi)*3^(3*n))$
 510
 511 /* end pochhammer tests */
 512
 513 /* binomial coefficient identities  */
 514 gamma_simp(makegamma(binomial(n,k) - (n/k)* binomial(n-1,k-1)));
 515 0$
 516
 517 gamma_simp(makegamma(k*binomial(n,k)  = n * binomial(n-1,k-1)));
 518 gamma(n+1)/(gamma(k)*gamma(n-k+1))=gamma(n+1)/(gamma(k)*gamma(n-k+1))$
 519
 520 gamma_simp(gamma(1/2 + %i*3/8)*gamma(1/2-%i*3/8));
 521 %pi/cosh((3*%pi)/8)$
 522
 523 /* Tests from Ray Rodgers */
 524 is(gamma_simp(makegamma(pochhammer(a+n,n) = pochhammer(a,2*n)/pochhammer(a,n))));
 525 true$
 526
 527 is(gamma_simp(makegamma(pochhammer(k*n+a,n) = pochhammer(a,(k+1)*n)/pochhammer(a,k*n))));
 528 true$
 529
 530 /* gamma_simp doesn't yield a canonical form for this case, but it does
 531 simplify the difference to zero.*/
 532
 533 featurep(n,'integer);
 534 true$
 535
 536 trigexpand(expand(gamma_simp(makegamma(pochhammer(a-n,n) - pochhammer(1-a,n)*(-1)^n))));
 537 0$
 538
 539 /* Same story for this case. But this case needs a final expand--OK*/
 540 trigexpand(expand(gamma_simp(makegamma( pochhammer(a-k*n,n) - (-1)^n*pochhammer((1-a),k*n)/pochhammer(1-a,(k-1)*n)))));
 541 0$
 542
 543 trigexpand(expand(gamma_simp(makegamma(1/pochhammer(a-n,n)  - (-1)^n/pochhammer(1-a,n)))));
 544 0$
 545
 546 /* and more tests */
 547
 548 factorial_simp((1100*(23/10)!*(27/10)!*(29/10)!*(31/10)!)/(8091*(17/10)!*(19/10)!*(21/10)!*(33/10)!));
 549 1$
 550
 551 (xxx : gamma(4*n+1/2)/gamma(2*n-3/4),0);
 552 0$
 553
 554 float_approx_equal(subst(n=1.07,xxx), subst(n=1.07,gamma_simp(xxx)));
 555 true$
 556
 557 gamma_simp(xxx - gamma_simp(xxx));
 558 0$
 559
 560 /* Clean up */
 561 (remvalue(X,e,xxx,yyy),killcontext(cntx), remvalue(cntx), 0);
 562 0$
 563
 564 /* This test is fragile. But it does check if the gamma_simp code introduces
 565 any extraneous global variables */
 566 values;
 567 [nextlayerfactor, facsum_combine]$
 568
 569 /* Did gamma_simp introduce any new contexts?*/
 570 contexts;
 571 [initial,global]$
 572
 573 (reset(radsubstflag,ratfac,domain,m1pbranch,%piargs),0);
 574 0$
 575