tests/rtesthyp.mac

   1 /*******************************************************************************
   2  *
   3  *   Tests of hypergeometric functions
   4  *
   5  ******************************************************************************/
   6
   7 kill(all);
   8 done;
   9
  10 assume(notequal(z,0));
  11 [notequal(z, 0)];
  12
  13 /*******************************************************************************
  14  *
  15  *  Test of hypergeometric function 0F0
  16  *
  17  ******************************************************************************/
  18
  19 hgfred([],[],z);
  20 %e^z;
  21
  22 /* ---------- Specific values ----------------------------------------------- */
  23
  24 hgfred([],[],0);
  25 1;
  26 hgfred([],[],0.0);
  27 1.0;
  28 hgfred([],[],0.0b0);
  29 1.0b0;
  30
  31 hgfred([],[],1);
  32 %e;
  33 hgfred([],[],2);
  34 %e^2;
  35
  36 /*******************************************************************************
  37  *
  38  *   Hypergeometric 1F0
  39  *
  40  ******************************************************************************/
  41
  42 hgfred([a],[],z);
  43 (1-z)^-a;
  44
  45 /* ---------- Specific values ----------------------------------------------- */
  46
  47 hgfred([a],[],0);
  48 1;
  49 hgfred([a],[],0.0);
  50 1.0;
  51 hgfred([a],[],0.0b0);
  52 1.0b0;
  53
  54 /* --- For x = 1 the result depends on the sign of a --- */
  55
  56 /* a is negative -> 0 */
  57 hgfred([-1],[],1);
  58 0;
  59 hgfred([-1],[],1.0);
  60 0.0;
  61 hgfred([-1],[],1.0b0);
  62 0.0b0;
  63
  64 /* a is positive -> inf or infinity */
  65 limit(hgfred([1],[],x),x,1);
  66 infinity;
  67
  68 limit(hgfred([2],[],x),x,1);
  69 inf;
  70
  71 /* Maxima does not check the sign of a symbolic expression for a.
  72  * The result is always 0.
  73  */
  74
  75 assume(xp>0,xn<0);
  76 [xp>0, xn<0];
  77
  78 /* The following two cases should give a Maxima error. But they do not.
  79  * The reason is, that Maxima always simplifies 0^expr -> 0.
  80  * That is a problem of simpexpt and not of $hgfred.
  81  */
  82 errcatch(hgfred([xp],[],1));
  83 [];
  84 errcatch(hgfred([-xn],[],1));
  85 [];
  86
  87 forget(xp>0, xn<0);
  88 [xp>0, xn<0];
  89
  90 /*******************************************************************************
  91  *
  92  *   Hypergeometric 0F1
  93  *
  94  ******************************************************************************/
  95
  96 /* --- General simplification in terms of bessel_i and bessel_j --- */
  97
  98 hgfred([],[b],x);
  99 bessel_i(b-1,2*sqrt(x))*gamma(b)*x^(1/2-b/2);
 100
 101 hgfred([],[b],-x);
 102 bessel_j(b-1,2*sqrt(x))*gamma(b)*x^(1/2-b/2);
 103
 104 /* ---------- Specific values ----------------------------------------------- */
 105
 106 hgfred([],[b],0);
 107 1;
 108 hgfred([],[b], 0.0);
 109 1.0;
 110 hgfred([],[b], 0.0b0);
 111 1.0b0;
 112
 113 /* --- b = n/2 and z = x^2/4 --- */
 114
 115 /* The results are in terms of bessel_i and bessel_j and a half integral order.
 116  * We expand the Bessel functions and get results in terms of trigonometric
 117  * functions.
 118  */
 119
 120 hgfred([],[1/2],x^2/4),besselexpand:true;
 121 cosh(x);
 122 hgfred([],[1/2],-x^2/4),besselexpand:true;
 123 cos(x);
 124
 125 hgfred([],[-1/2], x^2/4), besselexpand:true;
 126 cosh(x)-x*sinh(x);
 127 hgfred([],[-1/2], -x^2/4), besselexpand:true;
 128 cos(x)+x*sin(x);
 129
 130 hgfred([],[3/2],x^2/4), besselexpand:true;
 131 sinh(x)/x;
 132 hgfred([],[3/2],-x^2/4), besselexpand:true;
 133 sin(x)/x;
 134
 135 hgfred([],[-3/2],x^2/4), besselexpand:true;
 136 (3*x*sinh(x)+(-x^2-3)*cosh(x))/-3;
 137 hgfred([],[-3/2],-x^2/4), besselexpand:true;
 138 (3*x*sin(x)+(3-x^2)*cos(x))/3;
 139
 140 /*******************************************************************************
 141  *
 142  *   Hypergeometric 1F1 - Confluent hypergeometric function
 143  *
 144  ******************************************************************************/
 145
 146 /* ---------- Specific values ----------------------------------------------- */
 147
 148 hgfred([a],[b],0);
 149 1;
 150 hgfred([a],[b],0.0);
 151 1.0;
 152 hgfred([a],[b],0.0b0);
 153 1.0b0;
 154
 155 /* ---------- 1F1(a; a; z) -------------------------------------------------- */
 156
 157 /* A&S 13.6.12 */
 158 hgfred([a],[a],z);
 159 exp(z);
 160
 161 /* ---------- 1F1(a; a+n; z) for n= 1, 2 ------------------------------------ */
 162
 163 /* A&S 13.6.10 */
 164 hgfred([a],[a+1],-z);
 165 a*z^(-a)*gamma_incomplete_lower(a,z);
 166
 167 hgfred([a],[a+1],-z),prefer_gamma_incomplete:true;
 168 (a*gamma_incomplete(a,z)-a*gamma(a))/-z^a;
 169
 170 hgfred([a],[a+2],-z);
 171 z^(-a-1)*((a^2+a)*gamma_incomplete_lower(a,z)*z+(-a^2-a)*gamma_incomplete_lower(a+1,z));
 172
 173 hgfred([a],[a+2],-z),prefer_gamma_incomplete:true;
 174 -z^(-a-1)*(((a^2+a)*gamma_incomplete(a,z)+(-a^2-a)*gamma(a))*z
 175                 +(-a^2-a)*gamma_incomplete(a+1,z)+(a^2+a)*gamma(a+1));
 176
 177 /* ---------- 1F1(a; a-n; z) for n = -1, -2 --------------------------------- */
 178
 179 hgfred([a],[a-1],z);
 180 z*%e^z/(a-1)+%e^z;
 181
 182 hgfred([a],[a-2],z);
 183 z^2*%e^z/((a-2)*(a-1))+2*z*%e^z/(a-2)+%e^z;
 184
 185 /* ---------- Test 1F1 polynomials ------------------------------------------ */
 186
 187 /* First the cases, when we have a symbolic negative integer */
 188
 189 declare(n,integer);
 190 done;
 191 assume(n>0);
 192 [n>0];
 193
 194 hgfred([-n],[1/2],z);
 195 (-1)^n*n!*hermite(2*n,sqrt(z))/(2*n)!;
 196
 197 hgfred([-n],[3/2],z);
 198 (-1)^n*2^((2*n+1)/-2+n-1/2)*n!*hermite(2*n+1,sqrt(z))/((2*n+1)!*sqrt(z));
 199
 200 hgfred([-n],[b],z);
 201 gamma(b)*n!*gen_laguerre(n,b-1,z)/gamma(n+b);
 202
 203 remove(n,integer);
 204 done;
 205 forget(n>0);
 206 [n>0];
 207
 208 /* Check for specific values */
 209
 210 hgfred([-2],[1],z);
 211 /*laguerre(2,z)$*/
 212 z^2/2-2*z+1;
 213
 214 hgfred([-7],[1],z);
 215 /*laguerre(7,z)$*/
 216 -z^7/5040+7*z^6/720-7*z^5/40+35*z^4/24-35*z^3/6+21*z^2/2-7*z+1;
 217
 218 hgfred([-7],[c],z);
 219 /*5040*gen_laguerre(7,c-1,z)*gamma(c)/gamma(c+7)$*/
 220 -z^7/(c*(c+1)*(c+2)*(c+3)*(c+4)*(c+5)*(c+6))
 221  +7*z^6/(c*(c+1)*(c+2)*(c+3)*(c+4)*(c+5))-21*z^5/(c*(c+1)*(c+2)*(c+3)*(c+4))
 222  +35*z^4/(c*(c+1)*(c+2)*(c+3))-35*z^3/(c*(c+1)*(c+2))+21*z^2/(c*(c+1))-7*z/c+1;
 223
 224 /* A&S 13.6.17
 225  * Note that this formula is in terms of He but we want just H.
 226  * A&S 22.5.18 gives the relationship between He and H:
 227  *
 228  * He(n,x) = 2^(-n/2)*H(n,x/sqrt(2))
 229  */
 230 hgfred([-4],[1/2],z^2/2);
 231 /*hermite(8,z/sqrt(2))/1680;*/
 232 z^8/105-4*z^6/15+2*z^4-4*z^2+1;
 233
 234 hgfred([-2],[1/2],z);
 235 /*hermite(4,sqrt(z))/12$*/
 236 4*z^2/3-4*z+1;
 237
 238 hgfred([-2],[3/2],z);
 239 /*hermite(5,sqrt(z))/(120*sqrt(z));*/
 240 4*z^2/15-4*z/3+1;
 241
 242 /* A&S 13.6.18 */
 243 hgfred([-4],[3/2],z^2/2);
 244 /*hermite(9,z/sqrt(2))/(15120*sqrt(2)*z)$*/
 245 z^8/945-4*z^6/105+2*z^4/5-4*z^2/3+1;
 246
 247 /* ---------- 1F1(a; 2*a; x) ------------------------------------------------ */
 248
 249 hgfred([a],[2*a],x);
 250 16^(a/2)*bessel_i((2*a-1)/2,x/2)*gamma((2*a+1)/2)*x^(1/2-a)*%e^(x/2)/2;
 251
 252 hgfred([a],[2*a],-x);
 253 16^(a/2)*bessel_i((2*a-1)/2,x/2)*gamma((2*a+1)/2)*x^(1/2-a)*%e^-(x/2)/2;
 254
 255 /* A&S 13.6.1 */
 256 hgfred([v+1/2],[2*v+1],2*%i*z);
 257 4^(v/2)*bessel_j(v,z)*gamma(v+1)*exp(%i*z)/z^v$
 258
 259 /* A&S 13.6.3 */
 260 hgfred([v+1/2],[2*v+1],2*z);
 261 4^(v/2)*bessel_i(v,z)*gamma(v+1)*exp(z)/z^v$
 262
 263 /* A&S 13.6.4 */
 264 expand(hgfred([n+1],[2*n+2],2*%i*z));
 265 sqrt(2)*4^(n/2)*bessel_j(n+1/2,z)*gamma(n+3/2)*z^(-n-1/2)*%e^(%i*z);
 266
 267 /* A&S 13.6.5 */
 268 expand(hgfred([-n],[-2*n],2*%i*z));
 269 bessel_j(-n-1/2,z)*gamma(1/2-n)*z^(n+1/2)*%e^(%i*z)/(sqrt(2)*4^(n/2));
 270
 271 /* A&S 13.6.6 */
 272 expand(hgfred([n+1],[2*n+2],2*z));
 273 sqrt(2)*4^(n/2)*bessel_i(n+1/2,z)*gamma(n+3/2)*z^(-n-1/2)*%e^z;
 274
 275 /* ---------- 1F1(n+1/2; m; z) for m, n integer ----------------------------- */
 276
 277 hgfred([3/2],[1],x);
 278 ((bessel_i(1,x/2)+bessel_i(0,x/2))*x+bessel_i(0,x/2))*%e^(x/2);
 279
 280 hgfred([3/2],[2],x);
 281 (bessel_i(1,x/2)+bessel_i(0,x/2))*%e^(x/2);
 282
 283 hgfred([3/2],[4],x);
 284 ((4*bessel_i(1,x/2)-4*bessel_i(0,x/2))*x+16*bessel_i(1,x/2))*%e^(x/2)/x^2;
 285
 286 /* ---------- 1F1(1/2-n,m,x) for m, n integer ------------------------------- */
 287
 288 hgfred([-1/2],[1],x);
 289 ((bessel_i(1,x/2)-bessel_i(0,x/2))*x+bessel_i(0,x/2))*%e^(x/2);
 290
 291 hgfred([-1/2],[2],x);
 292 ((2*bessel_i(1,x/2)-2*bessel_i(0,x/2))*x
 293         -bessel_i(1,x/2)+3*bessel_i(0,x/2))
 294         *%e^(x/2)/3;
 295
 296 hgfred([-1/2],[4],x);
 297 ((16*bessel_i(1,x/2)-16*bessel_i(0,x/2))*x^3
 298         +(40*bessel_i(0,x/2)-24*bessel_i(1,x/2))*x^2
 299         +(4*bessel_i(0,x/2)-20*bessel_i(1,x/2))*x-16*bessel_i(1,x/2))
 300         *%e^(x/2)
 301         /(35*x^2);
 302
 303 /* ---------- 1F1(a; 2*a+n) for n an integer -------------------------------- */
 304
 305 hgfred([1/3],[2*1/3+1],x);
 306 (2^(2/3)*gamma(5/6)*bessel_i(5/6,x/2)
 307         -2^(2/3)*bessel_i(-1/6,x/2)*gamma(5/6))
 308         *x^(1/6)*%e^(x/2)/-2;
 309
 310 hgfred([1/3],[2*1/3+2],x);
 311 (2^(2/3)*gamma(5/6)*bessel_i(11/6,x/2)
 312         -5*2^(2/3)*gamma(5/6)*bessel_i(5/6,x/2)
 313         +2^(8/3)*bessel_i(-1/6,x/2)*gamma(5/6))
 314         *x^(1/6)*%e^(x/2)/8;
 315
 316 /* ---------- 1F1(a; 2*a-n) for n an integer -------------------------------- */
 317
 318 hgfred([1/3],[2*1/3-1],x);
 319 (3*bessel_i(-1/6,x/2)+3*bessel_i(-7/6,x/2))
 320         *gamma(5/6)*x^(7/6)*%e^(x/2)/-2^(4/3);
 321
 322 hgfred([1/3],[2*1/3-2],x);
 323 (45*bessel_i(-1/6,x/2)+63*bessel_i(-7/6,x/2)+18*bessel_i(-13/6,x/2))
 324         *gamma(5/6)*x^(13/6)*%e^(x/2)/(7*2^(10/3));
 325
 326 /*******************************************************************************
 327  *
 328  *   Hypergeometric 2F0
 329  *
 330  ******************************************************************************/
 331
 332 /* ---------- Test 2F0 polynomials ------------------------------------------ */
 333
 334 declare(n,integer);
 335 done;
 336 assume(n>0);
 337 [n>0];
 338
 339 hgfred([-n,-n+1/2],[],-1/z^2),radexpand:all;
 340 hermite(2*n,z)/(2^(2*n)*z^(2*n));
 341
 342 hgfred([-n,-n-1/2],[],-1/z^2),radexpand:all;
 343 2^((-2*n-1)/2-(2*n+1)/2)*hermite(2*n+1,z)*z^(-2*n-1);
 344
 345 hgfred([-n,b],[],z);
 346 n!*gen_laguerre(n,-n-b,-1/z)*z^n;
 347
 348 remove(n,integer);
 349 done;
 350 forget(n>0);
 351 [n>0];
 352
 353 /* Test 2F0 polynomials for specific values */
 354
 355 ratsimp(hgfred([-3,-3+1/2],[],z));
 356 /*-hermite(6,%i/sqrt(z))*z^3/64$*/
 357 (15*z^3+90*z^2+60*z+8)/8;
 358
 359 ratsimp(hgfred([-3-1/2,-3],[],z));
 360 /*%i*hermite(7,%i/sqrt(z))*z^(7/2)/128$*/
 361 (105*z^3+210*z^2+84*z+8)/8;
 362
 363 hgfred([-2,-4],[],z);
 364 /*2*gen_laguerre(2,2,-1/z)*z^2$*/
 365 12*(2/(3*z)+1/(12*z^2)+1)*z^2;
 366
 367 /* Test that 2F0([a,b],...) = 2F0([b,a],...) */
 368 hgfred([-4,-2],[],z);
 369 /*2*gen_laguerre(2,2,-1/z)*z^2$*/
 370 12*(2/(3*z)+1/(12*z^2)+1)*z^2;
 371
 372 /*******************************************************************************
 373  *
 374  *   Hypergeometric 2F1
 375  *
 376  ******************************************************************************/
 377
 378 /* ---------- 2F1(a,b; c; 0) ------------------------------------------------ */
 379
 380 hgfred([a,b],[c],0);
 381 1;
 382 hgfred([a,b],[c],0.0);
 383 1.0;
 384 hgfred([a,b],[c],0.0b0);
 385 1.0b0;
 386
 387 /* More examples from bug reports */
 388
 389 /* Bug ID: 1531688 - hgfred([ ], [ ], 0) */
 390 hgfred([3,4],[7],0);
 391 1;
 392 hgfred([1,1],[2],0);
 393 1;
 394
 395 /* Bug ID: 1857562 - hgfred([a,b],[c], 0) -/--> 1 */
 396 hgfred([3,9],[7],0);
 397 1;
 398 hgfred([3,9],[77],0);
 399 1;
 400
 401 /* ---------- 2F1(a,b; a; z) and 2F1(a,b; b; z) ----------------------------- */
 402
 403 hgfred([a,b],[a],z);
 404 (1-z)^-b;
 405
 406 hgfred([a,b],[b],z);
 407 (1-z)^-a;
 408
 409 /* ---------- 2F1(a,b; b-n; z) works for n a positive integer --------------- */
 410
 411 /* This is a test of the routine splitpfq */
 412
 413 hgfred([a,b],[b-1],z);
 414 a*(1-z)^(-a-1)*z/(b-1)+1/(1-z)^a;
 415
 416 hgfred([a,b],[b-2],z);
 417 a*(a+1)*(1-z)^(-a-2)*z^2/((b-2)*(b-1))+2*a*(1-z)^(-a-1)*z/(b-2)+1/(1-z)^a;
 418
 419 /* ---------- 2F1(a,b; b+1; z) ---------------------------------------------- */
 420
 421 hgfred([a,b],[b+1],z);            /* No soution, a noun form */
 422 %f[2,1]([a,b],[b+1],z);
 423
 424 /* ---------- 2F1(a,b; a+b-1/2; z) ------------------------------------------ */
 425
 426 /* This function is reduced to 2F1(a,b; a+b+1/2).
 427  * This is a test of the routines gered1 and legf20 called from legfun.*/
 428
 429 result:hgfred([a,b],[a+b-1/2],z);
 430 2^(b+a-3/2)*assoc_legendre_p(b-a-1/2,-b-a+3/2,sqrt(1-z))*gamma(b+a-1/2)
 431            *z^((-b-a+3/2)/2)
 432  /sqrt(1-z)$
 433
 434 /* Check for two special cases a=-1/2, b=2 and b=3 */
 435
 436 result - hgfred([-1/2,2],[-1/2+2-1/2],z), a=-1/2,b=2,ratsimp;
 437 0;
 438
 439 result - hgfred([-1/2,3],[-1/2+3-1/2],z), a=-1/2,b=3,ratsimp;
 440 0;
 441
 442 /* ---------- 2F1(a,b; a+b+1/2; z) ------------------------------------------ */
 443
 444 /* This is a test of routine legf20 called from legfun */
 445
 446 result:hgfred([a,b],[a+b+1/2],z);
 447 assoc_legendre_p(b-a-1/2,-b-a+1/2,sqrt(1-z))
 448  *2^(b+a-1/2)*gamma(b+a+1/2)*z^((-b-a+1/2)/2);
 449
 450 /* Check for two special cases a=2, b=-3/2 and a=3, b=-5/2 */
 451
 452 result - hgfred([2,-3/2],[2-3/2+1/2],z), a=2, b=-3/2, ratsimp;
 453 0;
 454
 455 result - hgfred([3,-5/2],[3-5/2+1/2],z), a=3, b=-5/2, ratsimp;
 456 0;
 457
 458 /* ---------- 2F1(a,b; a+b+3/2; z) ------------------------------------------ */
 459
 460 /* This is reduced to 2F1(a,b; a+b+1/2).
 461  * The implementation does not work as expected because the associated
 462  * legendre function is not fully implemented.
 463  */
 464
 465 result:hgfred([a,b],[a+b+3/2],z);
 466 (b+a+1/2)*(1-z)^(3/2)
 467          *(2^(b+a-3/2)*assoc_legendre_p(b-a-1/2,-b-a+1/2,sqrt(1-z))
 468                       *gamma(b+a+1/2)*z^((-b-a+1/2)/2)
 469           /(1-z)^(3/2)
 470           +(-b-a+1/2)*2^(b+a-3/2)*assoc_legendre_p(b-a-1/2,-b-a+1/2,sqrt(1-z))
 471                      *gamma(b+a+1/2)*z^((-b-a+1/2)/2-1)
 472            /sqrt(1-z)
 473           +2^(b+a-3/2)*gamma(b+a+1/2)
 474                       *(2*b*assoc_legendre_p(b-a+1/2,-b-a+1/2,sqrt(1-z))
 475                        -(b-a+1/2)*assoc_legendre_p(b-a-1/2,-b-a+1/2,sqrt(1-z))
 476                                  *sqrt(1-z))*z^((-b-a+1/2)/2-1)
 477            /(1-z))
 478  /((a+1/2)*(b+1/2));
 479
 480 /* --- Check for specific values --- */
 481
 482 /* These tests fail because the current implementation of the
 483  * associated legendre function is not complete.
 484  */
 485
 486 result - hgfred([-1,3/2],[-1+3/2+3/2],z), a=-1, b=3/2, ratsimp;
 487 0;
 488
 489 result - hgfred([-1,5/2],[-1+5/2+3/2],z), a=-1, b=5/2, ratsimp;
 490 0;
 491
 492 result - hgfred([-2,5/2],[-2+5/2+3/2],z), a=-2, b=5/2, ratsimp;
 493 0;
 494
 495 /* ---------- 2F1(a,b; a-b+1; z) -------------------------------------------- */
 496
 497 /* Test of the routine legf16 called from legfun.
 498  * Maxima distinguishes a solution for x<0 from a solution from x>0.
 499  * This has to be checked in more detail.
 500  */
 501
 502 assume(x>0);
 503 [x>0];
 504
 505 result:hgfred([a,b],[a-b+1],x);
 506 gamma(-b+a+1)*assoc_legendre_p(b-1,b-a,(x+1)/(1-x))*x^((b-a)/2)/(1-x)^b;
 507
 508 /* Check for some values */
 509
 510 result - hgfred([1,1],[1-1+1],x), a=1, b=1, ratsimp;
 511 0;
 512
 513 /* Fails because the associated legendre function is not fully implemented */
 514 result - hgfred([2,1],[2-1+1],x), a=2, b=1, ratsimp;
 515 0;
 516
 517 /* Fails because of different phase factor, we expect 1/(1-x)^3 */
 518 result - hgfred([3,2],[3-2+1],x), a=3, b=2, radcan;
 519 0;
 520
 521 /* Now the second result for a negative argument */
 522
 523 hgfred([a,b],[a-b+1],-x);
 524 gamma(-b+a+1)*assoc_legendre_p(b-1,b-a,(1-x)/(x+1))*x^((b-a)/2)/(x+1)^b;
 525
 526 forget(x>0);
 527 [x>0];
 528
 529 /* ---------- 2F1(a,b; a-b+2; z) -------------------------------------------- */
 530
 531 hgfred([a,b],[a-b+2],x);          /* No solution, a noun form */
 532 %f[2,1]([a,b],[-b+a+2],x);
 533
 534 /* ---------- 2F1(a,b; (a+b)/2; z) ------------------------------------------ */
 535
 536 hgfred([a,b],[(a+b)/2],x);        /* No solution, a noun form */
 537 %f[2,1]([a,b],[b/2+a/2],x);
 538
 539 /* ---------- 2F1(a,b; (a+b+1)/2; z) ---------------------------------------- */
 540
 541 /* Test of the routine legf14 and gered1 called from legfun.
 542  * Maxima distinguishes 0 < x < 1 from x < 0 and x > 1.
 543  * This has to be checked in more detail.
 544  * For this case we have no examples for specific values of a and b, because
 545  * the order and degree of the legendre function can not be an integer at
 546  * the same time, the legendre function does not simplify to a simple result.
 547  */
 548
 549 assume(x>0, x<1);
 550 [x>0, x<1];
 551
 552 hgfred([a,b],[(a+b+1)/2],x);
 553 assoc_legendre_p(b/2-a/2-1/2,-b/2-a/2+1/2,1-2*x)
 554  *gamma(b/2+a/2+1/2)*(1-x)^(-b/2+(b/2+a/2-1/2)/2-a/2+1/2)*x^((-b/2-a/2+1/2)/2);
 555
 556 forget(x>0, x<1);
 557 [x>0, x<1];
 558
 559 /* ---------- 2F1(a,b; (a+b)/2+1; z) ---------------------------------------- */
 560
 561 hgfred([a,b],[(a+b)/2+1],x);      /* No solution, noun form */
 562 %f[2,1]([a,b],[b/2+a/2+1],x);
 563
 564 /* ---------- 2F1(a,b; 2*b; z) ---------------------------------------------- */
 565
 566 /* Test of routine legf36 called from legfun
 567  * For this case Maxima does not distinguish 0 < x< 1 from x<0 and x>1.
 568  * This might be not correct. Further tests are needed.
 569  */
 570
 571 result:hgfred([a,b],[2*b],x);
 572 %e^-(%i*%pi*(b-a))*2^(1-(b-a)/2)*assoc_legendre_q(b-1,b-a,2/x-1)*gamma(2*b)
 573                   *(2/x-2)^((b-a)/2)*x^((b-a)/2-b)
 574  /(gamma(b)*gamma(2*b-a));
 575
 576 /* Check for specific values */
 577
 578 /* Fails because of different phase, we expect -log(1-x)/x */
 579 result - hgfred([1,1],[2*1],x), a=1, b=1, ratsimp;
 580 0;
 581
 582 /* Fails, we can not reproduce the expected result -3*(3*log(3)-4)/4 */
 583 result - hgfred([1,2],[2*2],-2), a=1, b=2, x=-2, ratsimp;
 584 0;
 585
 586 /* ---------- 2F1(a,b; 1/2; z) ---------------------------------------------- */
 587
 588 hgfred([a,b],[1/2],x);            /* No solution, noun form */
 589 %f[2,1]([a,b],[1/2],x);
 590
 591 /* ---------- 2F1(a,a+1/2; c; z) -------------------------------------------- */
 592
 593 /* Test of routine legf24 called from legfun.
 594  * Maxima distingushes x<0 and x>0. We get a problem with the phase of the
 595  * result. More tests are needed.
 596  */
 597
 598 assume(x>0);
 599 [x>0];
 600
 601 result:hgfred([a,a+1/2],[c],x);
 602 2^(c-1)*gamma(c)*assoc_legendre_p(c-2*a-1,1-c,1/sqrt(1-x))*(1-x)^((1-c)/-2-a)
 603        *x^((1-c)/2);
 604
 605 /* Fails because of different sign. The expected result is (x+2)/2 */
 606 result - hgfred([-3/2,-3/2+1/2],[3],x), a=-3/2, c=3, ratsimp;
 607 0;
 608
 609 forget(x>0);
 610 [x>0];
 611
 612 /* ---------- 2F1(a,1-a; c; z) ---------------------------------------------- */
 613
 614 /* Test for legf14 called from legfun.
 615  * Maxima distinguishes 0 < x < 1 from x>1 and x<1.
 616  * We look only at the result for 0 < x < 1. More tests are needed.
 617  */
 618
 619 assume(x>0, x<1);
 620 [x>0, x<1];
 621
 622 result:hgfred([a,1-a],[c],x);
 623 assoc_legendre_p(a-1,1-c,1-2*x)*gamma(c)*(1-x)^((c-1)/2)*x^((1-c)/2);
 624
 625 /* We can verify the general result for some specific values */
 626
 627 result - hgfred([2,1-2],[1],x), a=2, c=1, ratsimp;
 628 0;
 629
 630 result - hgfred([3,1-3],[1],x), a=3, c=1, ratsimp;
 631 0;
 632
 633 forget(x>0, x<1);
 634 [x>0, x<1];
 635
 636 /* ---------- 2F1(a,2-a; c; z) ---------------------------------------------- */
 637
 638 hgfred([a,2-a],[c],x);            /* No solution, noun form */
 639 %f[2,1]([a,2-a],[c],x);
 640
 641 /* ---------- 2F1(a,3-a; c; z) ---------------------------------------------- */
 642
 643 hgfred([a,3-a],[c],x);            /* No solution, noun form */
 644 %f[2,1]([a,3-a],[c],x);
 645
 646 /* ---------- 2F1(-n,b; c; z) ----------------------------------------------- */
 647
 648 /* We have extended 2f1polys to handle symbolic negative integers.
 649  * Furthermore, we can suppress the expansion of polynomials with the user
 650  * flag expand_polynomials.
 651  */
 652
 653 declare(n,integer);
 654 done;
 655 assume(n>0);
 656 [n>0];
 657
 658 hgfred([-n,b],[c],z),expand_polynomials:false;
 659 n!*jacobi_p(n,c-1,-n-c+b,1-2*z)/pochhammer(c,n);
 660
 661 remove(n,integer);
 662 done;
 663 forget(n>0);
 664 [n>0];
 665
 666 /* ---------- 2F1(a,a/2+1; a/2; z) ------------------------------------------ */
 667
 668 /* Test of splitpfq */
 669
 670 hgfred([a,a/2+1],[a/2],z);
 671 2*(1-z)^(-a-1)*z+1/(1-z)^a;
 672
 673 /* ---------- 2F1(a,a+1/2; 2*a-1; z) ---------------------------------------- */
 674
 675 /* Test of step7 */
 676
 677 result:hgfred([a,a+1/2],[2*a-1],z),ratsimp;
 678 ((2-2*a)*2^(2*a)*z^4+sqrt(1-z)*((12*a-11)*2^(2*a)*z^3+(58-76*a)*2^(2*a)*z^2
 679                                                      +(128*a-80)*2^(2*a)*z
 680                                                      +(32-64*a)*2^(2*a))
 681                     +(38*a-32)*2^(2*a)*z^3+(94-132*a)*2^(2*a)*z^2
 682                     +(160*a-96)*2^(2*a)*z+(32-64*a)*2^(2*a))
 683  /((4-8*a)*(sqrt(1-z)+1)^(2*a)*z^4+sqrt(1-z)
 684                                    *((32*a-16)*(sqrt(1-z)+1)^(2*a)*z^3
 685                                     +(64-128*a)*(sqrt(1-z)+1)^(2*a)*z^2
 686                                     +(160*a-80)*(sqrt(1-z)+1)^(2*a)*z
 687                                     +(32-64*a)*(sqrt(1-z)+1)^(2*a))
 688                                   +(80*a-40)*(sqrt(1-z)+1)^(2*a)*z^3
 689                                   +(100-200*a)*(sqrt(1-z)+1)^(2*a)*z^2
 690                                   +(192*a-96)*(sqrt(1-z)+1)^(2*a)*z
 691                                   +(32-64*a)*(sqrt(1-z)+1)^(2*a));
 692
 693 /* A check for specific values is correct, but we do not get the simple expected
 694  * result 2^(-3/2). The answer from Maxima is
 695  *    (3465*2^(5/2)+19601)/(19601*2^(3/2)+55440),
 696  * which looks much more complicate. */
 697 result - hgfred([1,1+1/2],[2*1-1],-1), a=1, z=-1, ratsimp;
 698 0;
 699
 700 /* ---------- 2F1(a,a+1/2; 2*a; z) ------------------------------------------ */
 701
 702 /* Test of routine hyp-cos */
 703
 704 hgfred([a,a+1/2],[2*a],z);
 705 2^(2*a-1)*(sqrt(1-z)+1)^(1-2*a)/sqrt(1-z);
 706
 707 /* ---------- 2F1(a,a+1/2; 2*a+1; z) ---------------------------------------- */
 708
 709 /* Test of routine hyp-cos */
 710
 711 hgfred([a,a+1/2],[2*a+1],z);
 712 2^(2*a)/(sqrt(1-z)+1)^(2*a);
 713
 714 /* ---------- 2F1(a,a+1/2; 1/2; z) ------------------------------------------ */
 715
 716 /* Test of routine trig-log */
 717
 718 assume(x>0);
 719 [x>0];
 720
 721 /* A&S 15.1.9 */
 722 hgfred([a,a+1/2],[1/2],z^2);
 723 1/2*((z+1)^(-2*a)+(1-z)^(-2*a));
 724
 725 hgfred([a,a+1/2],[1/2],x);
 726 (1/(sqrt(x)+1)^(2*a)+1/(1-sqrt(x))^(2*a))/2;
 727
 728 hgfred([a,a+1/2],[1/2],-x);
 729 cos(2*a*atan(sqrt(x)))/(x+1)^a;
 730
 731 forget(x>0);
 732 [x>0];
 733
 734 /* ---------- 2F1(a,a+1/2; 3/2; z) ------------------------------------------ */
 735
 736 /* A&S 15.1.10 */
 737 hgfred([a,a+1/2],[3/2],z^2);
 738 ((z+1)^(1-2*a)-(1-z)^(1-2*a))/2/(1-2*a)/z;
 739
 740 /* ---------- 2F1(a,-a; 1/2; z) --------------------------------------------- */
 741
 742 /* Test of trig-log */
 743
 744 assume(x>0);
 745 [x>0];
 746
 747 hgfred([a,-a],[1/2],x);
 748 cos(2*a*asin(sqrt(x)));
 749
 750 forget(x>0);
 751 [x>0];
 752
 753 /* For a negative integer we get a polynomial */
 754
 755 declare(n,integer);
 756 done;
 757 assume(n>0);
 758 [n>0];
 759
 760 hgfred([n,-n],[1/2],z);
 761 chebyshev_t(n,1-2*z);
 762
 763 remove(n,integer);
 764 done;
 765 forget(n>0);
 766 [n>0];
 767
 768
 769
 770
 771
 772
 773 /* ---------- 2F1(n,m; k; z), n,m,k integers -------------------------------- */
 774
 775 declare(n,integer);
 776 done;
 777
 778 /* For three integers result in terms of log(1-x) */
 779
 780 hgfred([1,2],[3],x);
 781 -2*(log(1-x)/x+1/x+log(1-x)*(1-x)/x^2);
 782
 783 hgfred([1,2],[4],x);
 784 -3*(-1/((1-x)*x)-2*log(1-x)/x^2-2/x^2-2*log(1-x)*(1-x)/x^3)*(1-x);
 785
 786 /* a or b can be a symbol declared to be an integer */
 787
 788 hgfred([1,n],[3],x);
 789 2*((1/x-1)/(n-1)+1)/((2-n)*x)+2*(1-x)^(2-n)/((n-2)*(n-1)*x^2);
 790
 791 hgfred([n,2],[3],x);
 792 2*((1/x-1)/(2-n)+1)*(1-x)^(1-n)/((n-1)*x)+2/((1-n)*(2-n)*x^2);
 793
 794 /* a,b,c integer and z=-1 */
 795
 796 hgfred([1,1],[2],-1);
 797 log(2);
 798
 799 hgfred([1,1],[3],-1);
 800 4*(log(2)-1/2);
 801
 802 hgfred([1,1],[4], -1);
 803 6*(2*log(2)-5/4);
 804
 805 hgfred([1,1],[5],-1);
 806 16*(6*log(2)-4)/3;
 807
 808 hgfred([1,2],[3], -1);
 809 -2*(log(2)-1);
 810
 811 /* a,b,c integer and z=1/2 */
 812
 813 hgfred([1,3],[2],1/2);
 814 3;
 815
 816 hgfred([1,4],[2],1/2);
 817 14/3;
 818
 819 hgfred([1,4],[3],1/2);
 820 8/3;
 821
 822 hgfred([1,5],[2],1/2);
 823 15/2;
 824
 825 /* a,b,c integer and z=1 */
 826
 827 hgfred([1,2],[4],1);
 828 3;
 829
 830 hgfred([1,2],[5],1);
 831 2;
 832
 833 /* Bug ID: 1869296 - hgfred([1,2],[6],1) bogus*/
 834 hgfred([1,2],[6],1);
 835 5/3;
 836
 837
 838 /* ---------- More tests ---------------------------------------------------- */
 839
 840 /* Tests of confluent hypergeometric function */
 841
 842 /* See A&S 13.6 */
 843 /* Bessel functions */
 844
 845 /* A&S 13.6.13 */
 846 besselexpand:true;
 847 true$
 848
 849 expand(hgfred([1],[2],-2*%i*z));
 850 exp(-%i*z)*sin(z)/z$
 851
 852 /* A&S 13.6.14 */
 853 expand(hgfred([1],[2],2*z));
 854 exp(z)*sinh(z)/z$
 855
 856 /* A&S 13.6.19 */
 857 hgfred([1/2],[3/2],-z^2);
 858 sqrt(%pi)/2/abs(z)*erf(abs(z))$
 859
 860 /* A&S 7.1.21 */
 861 hgfred([1],[3/2],z^2);
 862 sqrt(%pi)*exp(z^2)*erf(abs(z))/2/abs(z)$
 863
 864 /* Hypergeometric function 2F1 tests
 865 *
 866 * See A&S 15.1
 867 *
 868 */
 869 /* A&S 15.1.3 */
 870 hgfred([1,1],[2],z);
 871 -log(1-z)/z;
 872
 873 /* A&S 15.1.4 */
 874 hgfred([1/2,1], [3/2], z^2);
 875 log((1+z)/(1-z))/(2*z)$
 876
 877 /* A&S 15.1.5 */
 878 hgfred([1/2,1], [3/2], -z^2);
 879 atan(z)/z$
 880
 881 /* A&S 15.1.6 */
 882 hgfred([1/2,1/2], [3/2], z^2);
 883 asin(z)/z$
 884
 885 hgfred([1,1], [3/2], z^2);
 886 asin(z)/z/sqrt(1-z^2)$
 887
 888 /* A&S 15.1.7 */
 889 hgfred([1/2,1/2], [3/2], -z^2);
 890 log(z+sqrt(1+z^2))/z$
 891
 892 hgfred([1,1], [3/2], -z^2);
 893 log(z+sqrt(1+z^2))/z/sqrt(1+z^2)$
 894
 895 /* A&S 15.1.8 */
 896 hgfred([a,b],[b],z);
 897 (1-z)^(-a)$
 898
 899 /* A&S 15.1.11 */
 900 hgfred([-a,a],[1/2],-z^2);
 901 ((z+sqrt(1+z^2))^(2*a) + (sqrt(1+z^2)-z)^(2*a))/2$
 902
 903 /* A&S 15.1.12 */
 904 hgfred([a,1-a],[1/2],-z^2);
 905 ((sqrt(1+z^2)+z)^(2*a-1)+(sqrt(1+z^2)-z)^(2*a-1))/2/sqrt(1+z^2)$
 906
 907 /* A&S 15.1.13 */
 908 hgfred([a,a+1/2],[1+2*a],z);
 909 2^(2*a)/(1+sqrt(1-z))^(2*a)$
 910
 911 hgfred([a+1,a+1/2],[1+2*a],z);
 912 2^(2*a)/(1+sqrt(1-z))^(2*a)/sqrt(1-z)$
 913
 914 /* A&S 15.1.14 */
 915 hgfred([a,a+1/2],[2*a],z);
 916 2^(2*a-1)/sqrt(1-z)*(1+sqrt(1-z))^(1-2*a)$
 917
 918 /* A&S 15.1.15 */
 919 hgfred([a,1-a],[3/2],sin(z)^2),triginverses:all;
 920 sin((1-2*a)*z)/sin(z)/(1-2*a)$
 921
 922 /* A&S 15.1.16 */
 923 hgfred([a,2-a],[3/2],sin(z)^2),triginverses:all;
 924 sin((2-2*a)*z)/(2-2*a)/(sin(z)*cos(z))$
 925
 926 /* A&S 15.1.17 */
 927 assume(not(equal(sin(z),0)));
 928 [not equal(sin(z),0)]$
 929
 930 hgfred([-a,a],[1/2],sin(z)^2),triginverses:all;
 931 cos(2*a*z)$
 932
 933 /* A&S 15.1.18 */
 934 hgfred([a,1-a],[1/2],sin(z)^2),triginverses:all;
 935 cos((2*a-1)*z)/cos(z)$
 936
 937 /* A&S 15.1.19 */
 938 assume(not(equal(tan(z),0)));
 939 [not equal(tan(z), 0)]$
 940
 941 hgfred([a,a+1/2],[1/2],-tan(z)^2),triginverses:all;
 942 cos(z)^(2*a)*cos(2*a*z)$
 943
 944 /* A&S 15.4.3 */
 945 assume(x > 0);
 946 [x > 0]$
 947
 948 hgfred([-n,n],[1/2],x);
 949 cos(2*n*asin(sqrt(x)))$
 950
 951 hgfred([-10,10],[1/2],x),expand_polynomials:false;
 952 chebyshev_t(10,1-2*x)$
 953
 954 /* A&S 15.4.4
 955  * (And A&S 8.2.1)
 956  */
 957 assume(x1>0,x1<1);
 958 [x1 > 0, x1 < 1]$
 959
 960 hgfred([-n,n+1],[1],x1);
 961 legendre_p(n,1-2*x1);
 962
 963 /* A&S 8.1.2
 964  * Is there a bug in that formula?
 965  * http://functions.wolfram.com/HypergeometricFunctions/LegendreP2General/26/01/02/
 966  * shows the formula with (1+z) and (1-z).
 967  *
 968  * But this differs from formula 14 in Higher Transcendental
 969  * Functions, which agrees with A&S.
 970  *
 971  * Formula 14 says:
 972  *
 973  * P(v,u,z) = (z+1)^(u/2)*(z-1)^(-u/2)/gamma(1-u)*F(-v, 1+v; 1-u,(1-z)/2)
 974  *
 975  * Now what?
 976  *
 977  * This tests the legf14 function.
 978  */
 979 assume(z1>1);
 980 [z1 > 1]$
 981 radcan(((z1+1)/(z1-1))^(u/2)/gamma(1-u)*hgfred([-v,v+1],[1-u],(1-z1)/2));
 982 assoc_legendre_p(v,u,z1)$
 983
 984 /* A&S 15.4.16
 985  *
 986  * F(a,1-a;c;x) = gamma(c)*(-x)^(1/2-c/2)*(1-x)^(c/2-1/2)*P(-a,1-c,1-2*x)
 987  *
 988  * This tests the legf14 function.
 989  */
 990 radcan(hgfred([a,1-a],[c],z1)*(-z1)^(c/2-1/2)*(1-z1)^(1/2-c/2)/gamma(c));
 991 assoc_legendre_p(a-1, 1-c, 1-2*z1)$
 992
 993 /* A&S 15.4.10
 994  *
 995  * F(a,a+1/2;c;z) = 2^(c-1)*gamma(c)*z^(1/2-c/2)*(1-z)^(c/2-a-1/2)*P(2*a-c,1-c,1/sqrt(1-z))
 996  *
 997  * Test for legf24 function.
 998  */
 999 assume(zg1 > 1);
1000 [zg1 > 1];
1001 radcan(hgfred([a,a+1/2],[c],zg1)*2^(1-c)*zg1^(c/2-1/2)*(1-zg1)^(1/2+a-c/2)/gamma(c));
1002 /* Use A&S 8.2.1 to get this */
1003 assoc_legendre_p(c-2*a-1,1-c,1/sqrt(1-zg1))$
1004
1005 radcan(2^u*(zg1^2-1)^(-u/2)*zg1^(v+u)/gamma(1-u)*hgfred([-v/2-u/2,1/2-v/2-u/2],[1-u],1-1/zg1^2));
1006 assoc_legendre_p(v,u,zg1)$
1007
1008 /* A&S 15.4.12
1009  *
1010  * F(a,b;a+b+1/2;z) = 2^(a+b-1/2)*gamma(a+b+1/2)*(-z)^(1/2-a-b)*P(a-b-1/2,1/2-a-b,sqrt(1-z))
1011  * Test for legf20 function.
1012  */
1013 assume(x2>1);
1014 [x2 > 1]$
1015 expand(ratsimp(2^(1/2-a-b)*(x2)^(1/2*(a+b-1/2))/gamma(a+b+1/2)*hgfred([a,b],[a+b+1/2],x2)));
1016 assoc_legendre_p(b-a-1/2,1/2-a-b,sqrt(1-x2))$
1017
1018 2^m*(1-z1^2)^(-m/2)/gamma(1-m)*hgfred([1/2+n/2-m/2, -n/2-m/2], [1-m], 1-z1^2);
1019 assoc_legendre_p(n,m,z1)$
1020
1021 /*
1022  * Formula 36:
1023  * exp(-%i*m*%pi)*Q(n,m,z) =
1024  *     2^n*gamma(1+n)*gamma(1+n+m)*(z+1)^(m/2-n-1)*(z-1)^(-m/2)/gamma(2+2*n)
1025  *     *hgfred([1+n-m,1+n],[2+2*n],2/(1+z))
1026  *
1027  * Test for legf36 function.
1028  */
1029 expand(radcan(hgfred([a,b],[2*b],z)*2^(-2*b)*sqrt(%pi)/gamma(b+1/2)*gamma(2*b-a)*z^b*(1-z)^((a-b)/2)*exp(%i*%pi*(b-a))));
1030 /*
1031  * Maxima doesn't expand gamma(2*b) which would cancel the other
1032  * terms, so we leave it in.
1033  */
1034 2*sqrt(%pi)*assoc_legendre_q(b-1,b-a,2/z-1)*gamma(2*b)/(2^(2*b)*gamma(b)*gamma(b+1/2))$
1035
1036 hgfred([a,b],[2*b],z) - hgfred([b,a],[2*b],z);
1037 0$
1038
1039 exp(%i*%pi*u)*2^v*gamma(1+v)*gamma(1+v+u)*(z+1)^(u/2-v-1)*(z-1)^(-u/2)/gamma(2+2*v)*hgfred([1+v-u,1+v],[2+2*v],2/(1+z));
1040 assoc_legendre_q(v, u, z)$
1041
1042 /*
1043  * A&S 15.4.14
1044  *
1045  * F(a,b;a-b+1;z) = gamma(a-b+1)*z^(b/2-a/2)*(1-z)^(-b)*P(-b,b-a,(1+z)/(1-z))
1046  *
1047  * Test for legf16 function
1048  */
1049 assume(zp1 > 1);
1050 [zp1 > 1]$
1051 hgfred([a,b],[a-b+1],zp1);
1052 gamma(-b+a+1)*assoc_legendre_p(b-1,b-a,(1+zp1)/(1-zp1))*zp1^((b-a)/2)*(1-zp1)^(-b)$
1053
1054 /*
1055  * Test for legf16 function
1056  */
1057 2^(-v)*(zp1+1)^(u/2+v)*(zp1-1)^(-u/2)/gamma(1-u)*ratsimp(hgfred([-v,-v-u],[1-u],(zp1-1)/(zp1+1)));
1058 assoc_legendre_p(v, u, zp1)$
1059
1060 /* Regression tests */
1061
1062 /* simpg tests */
1063
1064 /* F([...,0,...],[non-zero], z) is 1 */
1065 hgfred([a,0],[c], z);
1066 1$
1067
1068 /* F([non-zero],[...,zero or neg integer,...],z) is undefined */
1069 hgfred([a,b],[0],z);
1070 %f[2,1]([a,b],[0],z)$
1071
1072 hgfred([a,b],[c,4,-3],z);
1073 %f[2,3]([a,b],[c,4,-3],z)$
1074
1075 /* F([...,neg integer,...],[...],z) is a polynomial */
1076
1077
1078 /* Test 2F1 polynomials */
1079
1080 hgfred([-2,1],[1/2],z);
1081 /*8*jacobi_p(2,-1/2,-3/2,1-2*z)/3$*/
1082 8*z^2/3-4*z+1$
1083
1084 /* F(a,b;c;z) = F(b,a;c;z) */
1085 hgfred([1,2],[3],z);
1086 -2*(log(1-z)/z+1/z+log(1-z)*(1-z)/z^2)$
1087
1088 hgfred([1,2],[3],z)-hgfred([2,1],[3],z);
1089 0$
1090
1091 /*
1092  * Test F(a,b;c;z) for some integral values of a,b,c.
1093  * These can all be written in terms of F(1,1;2;z).
1094  */
1095
1096 /*
1097  * A&S 15.2.6:
1098  * diff((1-z)^(a+b-c)*F(a,b;c;z),z,n) = poch(c-a,n)*poch(c-b,n)/poch(c,n)*
1099  *    (1-z)^(a+b-c-n)*F(a,b;c+n;z)
1100  *
1101  * With a=b=1 and c = 2, we have
1102  *
1103  * diff(F(1,1;2;z),z,n) = poch(1,n)*poch(1,n)/poch(2,n)*(1-z)^(-n)*F(1,1;2+n,z)
1104  *
1105  * And A&S 15.1.3 says F(1,1;2;z) = -log(1-z)/z.
1106  *
1107  * So F(1,1;2+n;z) = diff(F(1,1;2;z),z,n)*(1-z)^n*poch(2,n)/poch(1,n)/poch(1,n)
1108  *
1109  * F(1,1;3;z) = diff(F(1,1;2;z),z,1)*(1-z)*2/1!/1!
1110  *            = 2*(1-z)*diff(F(1,1;2;z),z,1)
1111  * F(1,1;4;z) = diff(F(1,1;2;z),z,2)*(1-z)^2*(2*3)/2!/2!
1112  *            = 3/2*(1-z)^2*diff(F(1,1;2;z),z,2);
1113  */
1114
1115 hgfred([1,1],[3],z);
1116 2*(1/((1-z)*z)+log(1-z)/z^2)*(1-z)$
1117
1118 hgfred([1,1],[4],z);
1119 3*(1/((1-z)^2*z)-2/((1-z)*z^2)-2*log(1-z)/z^3)*(1-z)^2/2$
1120
1121 /*
1122  * A&S 15.2.7:
1123  * diff((1-z)^(a+n-1)*F(a,b;c;z),z,n)
1124  *    = (-1)^n*poch(a,n)*poch(c-b,n)/poch(c,n)*(1-z)^(a-1)*F(a+n,b;c+n;z)
1125  *
1126  * So
1127  * F(a+n,b;c+n;z) = (-1)^n*poch(c,n)/poch(a,n)/poch(c-b,n)*(1-z)^(1-a)
1128  *                  *diff((1-z)^(a+n-1)*F(a,b;c;z),z,n)
1129  *
1130  * Thus,
1131  * F(1+n,1;2+n;z) = (-1)^n*poch(2,n)/n!/n!*diff((1-z)^n*F(1,1;2;z),z,n)
1132  *
1133  * F(2,1;3;z) = -1*2*diff((1-z)*F(1,1;2;z),z,1)
1134  * F(3,1;4;z) = (2*3)/2!/2!*diff((1-z)^2*F(1,1;2;z),z,2)
1135  */
1136 hgfred([2,1],[3],z);
1137 -2*(log(1-z)/z+1/z+log(1-z)*(1-z)/z^2)$
1138
1139 hgfred([3,1],[4],z);
1140 3*(-2*log(1-z)/z-3/z-4*log(1-z)*(1-z)/z^2-2*(1-z)/z^2
1141                        -2*log(1-z)*(1-z)^2/z^3)
1142         /2$
1143
1144 /*
1145  * A&S 15.2.2:
1146  * diff(F(a,b;c;z),z,n) = poch(a,n)*poch(b,n)/poch(c,n)*F(a+n,b+n;c+n;z)
1147  *
1148  * So
1149  *
1150  * F(a+n,b+n;c+n;z) = poch(c,n)/poch(a,n)/poch(b,n)*diff(F(a,b;c;z),z,n)
1151  *
1152  * F(1+n,1+n;2+n;z) = poch(2,n)/n!/n!*diff(F(1,1;2;z),z,n)
1153  *
1154  * F(2,2;3;z) = 2*diff(F(1,1;2;z),z,1);
1155  * F(3,3;4;z) = 2*3/2!/2!*diff(f(1,1;2;z),z,2)
1156  */
1157
1158 hgfred([2,2],[3],z);
1159 2*(1/((1-z)*z)+log(1-z)/z^2)$
1160
1161 hgfred([3,3],[4],z);
1162 3*(1/((1-z)^2*z)-2/((1-z)*z^2)-2*log(1-z)/z^3)/2$
1163
1164 /*
1165  * A&S 15.3.9:
1166  * F(a,b;c;z) = A*z^(-a)*F(a,a-c+1;a+b-c+1;1-1/z)
1167  *              + B*(1-z)^(c-a-b)*z^(a-c)*F(c-a,1-a;c-a-b+1,1-1/z)
1168  *
1169  * where A = gamma(c)*gamma(c-a-b)/gamma(c-a)/gamma(c-b)
1170  *       B = gamma(c)*gamma(a+b-c)/gamma(a)/gamma(b)
1171  *
1172  * F(2,b;3;z) = A*z^(-2)*F(2,0;b;1-1/z) + B*(1-z)^(1-b)*z^(-1)*F(1,-1;2-b;1-1/z)
1173  * where A = gamma(3)*gamma(1-b)/gamma(1)/gamma(3-b) = 2*gamma(1-b)/gamma(3-b)
1174  *       B = gamma(3)*gamma(b-1)/gamma(2)/gamma(b) = 2*gamma(b-1)/gamma(b)
1175  *
1176  * And we know that F(2,0;b;1-1/z) = 1 and
1177  * F(1,-1;2-b;1-1/z) = jacobi_p(1,1-b,b-2,1-2*(1-1/z))/(2-b).
1178  *
1179  * This tests geredf in hyp.lisp.
1180  *
1181  */
1182 declare(ia,integer,ib,integer,ic,integer);
1183 done$
1184
1185 hgfred([2,ib],[3],z),gamma_expand:true; /*Expand gamma function */
1186 /*2*jacobi_p(1,1-ib,ib-2,1-2*(1-1/z))*gamma(ib-1)*(1-z)^(1-ib)/((2-ib)*gamma(ib)*z)
1187         +2*gamma(1-ib)/(gamma(3-ib)*z^2)$*/
1188 2*((1/z-1)/(2-ib)+1)*(1-z)^(1-ib)/((ib-1)*z)+2/((1-ib)*(2-ib)*z^2);
1189
1190 /*
1191  * Test of step4, which handles F(a,b;c;z) when a+b is an integer and
1192  * c + 1/2 is an integer.
1193  *
1194  */
1195 assume(zn < 0);
1196 [zn < 0]$
1197
1198 /*
1199  * Tests f81
1200  * First test.  F(a,-a;3/2;z)
1201  * A&S 15.2.6 tells us how to compute this:
1202  *
1203  * diff((1-z)^(a+b-c)*F(a,b;c;z),z,n)
1204  *     = poch(c-a,n)*poch(c-b,n)/poch(c,n)*(1-z)^(a+b-c-n)*F(a,b;c+n;z)
1205  *
1206  * So
1207  *
1208  * diff((1-z)^(-1/2)*F(a,-a;1/2;z),z,1)
1209  *     = poch(1/2-a,1)*poch(1/2+a,n)/poch(1/2,n)*(1-z)^(-3/2)*F(a,-a;3/2;z)
1210  *     = (1/2-a)*(1/2+a)/(1/2)*(1-z)^(-3/2)*F(a,-a;3/2;z)
1211  *
1212  * So
1213  * F(a,-a;3/2;z) = (1-z)^(3/2)*2/(1/2-a)/(1/2+a)*diff((1-z)^(-1/2)*F(a,-a;1/2;z),z,1)
1214  *
1215  * A&S 15.1.11 and 15.1.17 tells us what F(a,-a;1/2;z) is.
1216  *
1217  */
1218 hgfred([-a,a],[3/2],zn)-(1-zn)^(3/2)*(1/2)/(1/2-a)/(1/2+a)*diff((1-zn)^(-1/2)*hgfred([-a,a],[1/2],zn),zn,1);
1219 0$
1220
1221 /*
1222  * Test f85
1223  * F(a+2,-a;1/2;zn)
1224  *
1225  * A&S 15.2.3:
1226  * diff(z^(a+n-1)*F(a,b;c;z),z,n) = poch(a,n)*z^(a-1)*F(a+n,b;c;z)
1227  *
1228  * or
1229  *
1230  * F(a+n,b;c;z) = z^(1-a)/poch(a,n)*diff(z^(a+n-1)*F(a,b;c;z),z,n)
1231  *
1232  * In this case, we have
1233  *
1234  * F(a+2,-a;1/2;z) = z^(1-a)/poch(a,2)*diff(z^(a+1)*F(a,-a;c;z),z,2)
1235  */
1236 hgfred([a+2,-a],[1/2],zn)-zn^(1-a)/(a*(a+1))*diff(zn^(a+1)*hgfred([a,-a],[1/2],zn),zn,2);
1237 0$
1238
1239 /*
1240  * Tests f83
1241  * F(a-1,-a;1/2;z)
1242  *
1243  * F(a-1,-a;1/2;z)
1244  *
1245  * A&S 15.2.5:
1246  * diff(z^(c-a+m-1)*(1-z)^(a+b-c)*F(a,b;c;z),z,m)
1247  *   = poch(c-a,m)*z^(c-a-1)*(1-z)^(a+b-c-m)*F(a-m,b;c;z)
1248  *
1249  * So
1250  *
1251  * diff(z^(1/2-a)*(1-z)^(-1/2)*F(a,-a;1/2;z),z,1)
1252  *    = poch(1/2-a,1)*z^(-a-1/2)*(1-z)^(-3/2)*F(a-1,-a;1/2;z)
1253  *
1254  * F(a-1,-a;1/2;z)
1255  *    = (1-z)^(3/2)*z^(a+1/2)/(1/2-a)*diff(z^(1/2-a)*(1-z)^(-1/2)*F(a,-a;1/2;z),z,1)
1256  *
1257  */
1258 ratsimp(hgfred([a-1,-a],[1/2],zn)-(1-zn)^(3/2)*zn^(a+1/2)/(1/2-a)*diff(zn^(1/2-a)*(1-zn)^(-1/2)*hgfred([a,-a],[1/2],zn),zn,1));
1259 0$
1260
1261 /*
1262  * Tests f84
1263  * F(a,-a;-1/2;z)
1264  *
1265  * A&S 15.2.4:
1266  * diff(z^(c-1)*F(a,b;c;z),z,n) = poch(c-n,n)*z^(c-n-1)*F(a,b;c-n;z)
1267  *
1268  * A&S 15.2.3:
1269  * diff(z^(a+m-1)*F(a,b;c;z),z,m) = poch(a,n)*z^(a-1)*F(a+n,b;c;z)
1270  *
1271  * So:
1272  * diff(z^(-1/2)*F(a,-a;1/2;z),z,1) = poch(1/2-1,1)*z^(1/2-1-1)*F(a,-a;-1/2;z)
1273  *
1274  * F(a,-a;-1/2;z) = z^(3/2)/poch(-1/2,1)*diff(z^(-1/2)*F(a,-a;1/2;z),z,1)
1275  *
1276  */
1277 ratsimp(hgfred([a,-a],[-1/2],zn) - zn^(3/2)/(-1/2)*diff(zn^(-1/2)*hgfred([a,-a],[1/2],zn),zn,1));
1278 0$
1279
1280 /*
1281  * Tests f82
1282  * F(a-1,-a;-1/2;z)
1283  *
1284  * A&S 15.2.4
1285  * diff(z^(c-1)*F(a,b;c;z),z,n)
1286  *     = poch(c-n,n)*z^(c-n-1)*F(a;b;c-n;z)
1287  *
1288  * A&S 15.2.8:
1289  * diff(z^(c-1)*(1-z)^(b-c+n)*F(a,b;c;z),z,n)
1290  *     - poch(c-n,n)*z^(c-n-1)*(1-z)^(b-c)*F(a-n,b;c-n;z)
1291  *
1292  * So:
1293  * diff(z^(-1/2)*(1-z)^(-a+1/2)*F(a,-a;1/2;z),z,1)
1294  *      = poch(-1/2,1)*z^(-3/2)*(1-z)^(-a-1/2)*F(a-1,-a;-1/2;z)
1295  *
1296  * F(a-1,-a;-1/2;z) = z^(3/2)*(1-z)^(a+1/2)/poch(-1/2,1)*diff(z^(-1/2)*(1-z)^(-a+1/2)*F(a,-a;1/2;z),z,1)
1297  *
1298  */
1299 ratsimp(hgfred([a-1,-a],[-1/2],zn) - zn^(3/2)*(1-zn)^(a+1/2)/(-1/2)*diff(zn^(-1/2)*(1-zn)^(-a+1/2)*hgfred([a,-a],[1/2],zn),zn,1));
1300 0$
1301
1302 /*
1303  * Tests f86
1304  * F(a-2,-a;-1/2;z)
1305  *
1306  * A&S 15.2.5
1307  * diff(z^(c-a+n-1)*(1-z)^(a+b-c)*F(a,b;c;z),z,n)
1308  *     = poch(c-a,n)*z^(c-a-1)*(1-z)^(a+b-c-n)*F(a-n,b;c;z)
1309  *
1310  * A&S 15.2.8:
1311  * diff(z^(c-1)*(1-z)^(b-c+n)*F(a,b;c;z),z,n)
1312  *     = poch(c-n,n)*z^(c-n-1)*(1-z)^(b-c)*F(a-n,b;c-n;z)
1313  *
1314  * For our problem:
1315  *
1316  * diff(z^(-a+1-1/2)*(1-z)^(-1/2)*F(a,-a;1/2;z),z,1)
1317  *     = poch(1/2-a,1)*z^(-a-1/2)*(1-z)^(-1/2-1)*F(a-1,-a;1/2;z)
1318  *     = (1/2-a)*z^(-a-1/2)*(1-z)^(-3/2)*F(a-1,-a;-1/2;z)
1319  *
1320  * diff(z^(-1/2)*(1-z)^(-a-1/2+1)*F(a-1,-a;1/2;z),z,1)
1321  *     = poch(1/2-1,1)*z^(-1-1/2)*(1-z)^(-a-1/2)*F(a-2,-a;-1/2;z)
1322  *     = -1/2*z^(-3/2)*(1-z)^(-a-1/2)*F(a-2,-a;-1/2;z)
1323  *
1324  * Let
1325  * G(z) = z^(-a-1/2)*(1-z)^(-3/2)*F(a-1,-a;1/2;z)
1326  *      = (1/2-a)^(-1)*diff(z^(-a+1/2)*(1-z)^(-1/2)*F(a,-a;1/2;z),z,1)
1327  *
1328  * F(a-2,-a;-1/2;z)
1329  *      = z^(3/2)*(1-z)^(a+1/2)/(-1/2)
1330  *         *diff(z^(-1/2)*(1-z)^(-a+1/2)*F(a-1,-a;1/2;z),z,1)
1331  *      = z^(3/2)*(1-z)^(a+1/2)/(-1/2)
1332  *         *diff(z^a*(1-z)^(2-a)*G(z),z,n)
1333  *
1334  */
1335 ratsimp(hgfred([a-2,-a],[-1/2],zn) -
1336    zn^(3/2)*(1-zn)^(a+1/2)/(-1/2)*(1/2-a)^(-1)*
1337    diff(zn^a*(1-zn)^(2-a)
1338         *diff(zn^(-a+1/2)*(1-zn)^(-1/2)*hgfred([a,-a],[1/2],zn),zn,1)
1339         ,zn,1));
1340 0$
1341
1342 /*
1343  * F(1,2;3/2;z)
1344  *
1345  * This is easily derived from F(1,1;3/2;z) via A&S 15.2.3:
1346  *
1347  * diff(z^(a+n-1)*F(a,b;c;z),z,n) = poch(a,n)*z^(a-1)*F(a+n,b;c;z)
1348  *
1349  * And from A&S 15.1.6, F(1,1;3/2;z^2) = asin(z)/z/sqrt(1-z^2).
1350  *
1351  * So
1352  * diff(z*F(1,1;3/2;z),z,1) = poch(1,1)*F(2,1;3/2;z)
1353  *                          = F(1,2;3/2;z)
1354  *
1355  */
1356
1357 /* Tests step4-int */
1358 hgfred([1,2],[3/2],zn) - diff(zn*hgfred([1,1],[3/2],zn),zn,1);
1359 0$
1360
1361 /*
1362  * Some tests of splitpfq
1363  */
1364
1365 /*
1366  * From A&S 15.2.2
1367  *
1368  * diff(z^(a+n-1)*F(a,b;c;z),z,n) = poch(a,n)*z^(a-1)*F(a+n,b;c;z)
1369  *
1370  * So
1371  *
1372  * F(a+n,b;c;z) = z^(1-a)/poch(a,n)*diff(z^(a+n-1)*F(a,b;c;z),z,n)
1373  */
1374
1375 /*
1376  * F(a+3,b;a+2;z) = z^(1-a-2)/poch(a+2,1)*diff(z^(a+2+1-1)*F(a+2;b;a+2;z),z,1)
1377  *                = z^(-1-a)/(a+2)*diff(z^(a+2)*F(b;;z),z,1)
1378  * F(b;;z) = (1-z)^(-b)
1379  *
1380  * F(a+3,b;a+2,z) = z^(-1-a)/(a+2)*diff(z^(a+2)*(1-z)^(-b),z,1)
1381  */
1382 multthru(hgfred([a+3,b],[a+2],z));
1383 b*(1-z)^(-b-1)*z/(a+2)+1/(1-z)^b$
1384
1385 /*
1386  * F(a+4,b;a+2;z) = z^(1-a-2)/poch(a+2,2)*diff(z^(a+2+2-1)*F(a+2;b;a+2;z),z,2)
1387  *                = z^(-1-a)/(a+2)/(a+3)*diff(z^(a+3)*F(b;;z),z,2)
1388  *
1389  */
1390 ratsimp(hgfred([b,a+4],[a+2],z)-
1391 (b*(b+1)*(1-z)^(-b-2)*z^2/((a+2)*(a+3))
1392         +2*b*(1-z)^(-b-1)*z/(a+2)+1/(1-z)^b));
1393 0$
1394
1395 /*
1396  * Bug 1155241:  F(1/2,-1/2;3/2,z)
1397  *
1398  * Applying A&S 15.3.3 gives
1399  *
1400  *  F(1/2,-1/2;3/2,z) = (1-z)^(3/2)*F(1,2;3/2;z)
1401  *
1402  * And we know what F(1,2;3/2;z) is.
1403  *
1404  * Alternatively, apply A&S 15.2.5 to get
1405  *
1406  * diff(z*(1-z)^(-1/2)*F(1/2,1/2;3/2;z),z,1)
1407  *   = (1-z)^(-3/2)*F(-1/2,1/2;3/2;z)
1408  *
1409  * Both of these give the same answer:
1410  *
1411  * (sqrt(1-z)*sqrt(z)+asin(sqrt(z)))/(2*sqrt(z))
1412  *
1413  * when z > 0.
1414  */
1415 (assume(zp > 0), true);
1416 true$
1417
1418 ratsimp(hgfred([1/2,-1/2],[3/2],zp));
1419 (sqrt(1-zp)*sqrt(zp)+asin(sqrt(zp)))/(2*sqrt(zp))$
1420
1421 /*
1422  * Tests for step7
1423  */
1424
1425 /*
1426  * F(s,s+1/2;2*s+2;z) can be derived from F(s,s+1/2;2*s+1;z) via
1427  * A&S 15.2.6:
1428  *
1429  * F(a,b;c+n;z) = poch(c,n)/poch(c-a,n)/poch(c-b,n)*(1-z)^(c+n-a-b)
1430  *                 *diff((1-z)^(a+b-c)*F(a,b;c;z)),z,n)
1431  *
1432  * F(s,s+1/2;2*s+2;z)
1433  *    =  (2*s+1)/(s+1)/(s+1/2)*(1-z)^(3/2)
1434  *         *diff((1-z)^(-1/2)*F(s,s+1/2;2*s+1;z),z,1)
1435  *
1436  * (Tests step-7-pm)
1437  */
1438 hgfred([s,s+1/2],[2*s+2],z)
1439   - (2*s+1)/(s+1)/(s+1/2)*(1-z)^(3/2)
1440      *diff((1-z)^(-1/2)*hgfred([s,s+1/2],[2*s+1],z),z,1);
1441 0$
1442
1443 hgfred([s,s+1/2],[2*s+2],z) - hgfred([s+1/2,s],[2*s+2],z);
1444 0$
1445
1446 /*
1447  * F(s+2,s+1/2;2*s+1;z) can be derived from F(s,s+1/2;2*s+1;z) via
1448  * A&S 15.2.3:
1449  *
1450  * F(a+n,b;c;z) = z^(1-a)/poch(a,n)*diff(z^(a+n-1)*fun,z,n)
1451  *
1452  * F(s+2,s+1/2;2*s+1;z)
1453  *     = z^(1-s)/(s*(s+1))*diff(z^(s+1)*F(s,s+1/2;2*s+1;z),z,2)
1454  *
1455  * (Tests step-7-pp)
1456  */
1457 hgfred([s+2,s+1/2],[2*s+1],z)
1458   - z^(1-s)/s/(s+1)*diff(z^(s+1)*hgfred([s,s+1/2],[2*s+1],z),z,2);
1459 0$
1460
1461 hgfred([s+2,s+1/2],[2*s+1],z)-hgfred([s+1/2,s+2],[2*s+1],z);
1462 0$
1463
1464 /*
1465  * F(s-1,s+1/2;2*s+1;z) can be derived from F(s,s+1/2;2*s+1;z) via
1466  * A&S 15.2.5:
1467  *
1468  * F(a-n,b;c;z) = 1/poch(c-a,n)*z^(1+a-c)*(1-z)^(c+n-a-b)
1469  *                 *diff(z^(c-a+n-1)*(1-z)^(a+b-c)*F(a,b;c;z),z,n)
1470  *
1471  * F(s-1,s+1/2;2*s+1;z)
1472  *     = 1/(s+1)*z^(-s)*(1-z)^(3/2)
1473  *        *diff(z^(s+1)*(1-z)^(-1/2)*F(s,s+1/2;2*s+1;z),z,1)
1474  *
1475  * (Tests step-7-mp)
1476  */
1477 hgfred([s-1,s+1/2],[2*s+1],z)
1478  - 1/(s+1)*z^(-s)*(1-z)^(3/2)
1479     *diff(z^(s+1)*(1-z)^(-1/2)*hgfred([s,s+1/2],[2*s+1],z),z,1);
1480 0$
1481
1482 hgfred([s-1,s+1/2],[2*s+1],z) - hgfred([s+1/2,s-1],[2*s+1],z);
1483 0$
1484
1485 /*
1486  * F(s-1,s+1/2;2*s+2;z) can be derived from F(s,s+1/2;2*s+1;z) via
1487  * A&S 15.2.5:
1488  *
1489  * F(a-n,b;c;z) = 1/poch(c-a,n)*z^(1+a-c)*(1-z)^(c+n-a-b)
1490  *                 *diff(z^(c-a+n-1)*(1-z)^(a+b-c)*F(a,b;c;z),z,n)
1491  *
1492  * and A&S 15.2.6:
1493  *
1494  * F(a,b;c+n;z) = poch(c,n)/poch(c-a,n)/poch(c-b,n)*(1-z)^(c+n-a-b)
1495  *                 *diff((1-z)^(a+b-c)*F(a,b;c;z)),z,n)
1496  *
1497  *
1498  * F(s-1,s+1/2;2*s+2;z)
1499  *    = 1/(s+2)*z^(-s-1)*(1-z)^(5/2)
1500  *       *diff(z^(s+2)*(1-z)^(-3/2)*F(s,s+1/2;2*s+2;z)),z,1)
1501  *
1502  * F(s,s+1/2;2*s+2;z)
1503  *    =  (2*s+1)/(s+1)/(s+1/2)*(1-z)^(3/2)
1504  *       *diff((1-z)^(-1/2)*F(s,s+1/2;2*s+1;z),z,1)
1505  * (Tests step-7-mm)
1506  */
1507 hgfred([s-1,s+1/2],[2*s+2],z)
1508  - 1/(s+2)*z^(-s-1)*(1-z)^(5/2)
1509    *diff(z^(s+2)*(1-z)^(-3/2)
1510          *(2*s+1)/(s+1)/(s+1/2)*(1-z)^(3/2)
1511           *diff((1-z)^(-1/2)*hgfred([s,s+1/2],[2*s+1],z),z,1)
1512          ,z,1);
1513 0$
1514
1515 hgfred([s-1,s+1/2],[2*s+2],z) - hgfred([s+1/2,s-1],[2*s+2],z);
1516 0$
1517
1518 (assume(p>0),true);
1519 true$
1520 /* Table of Integral Transforms */
1521 /* Formula 36, p 147 */
1522 specint(exp(-a*exp(-t))*exp(-p*t),t);
1523 a^(-p)*gamma_incomplete_lower(p,a)$
1524
1525 /* Formula 37, p 147 */
1526 specint(exp(-a*exp(t))*exp(-p*t),t);
1527 a^p*gamma_incomplete(-p,a)$
1528
1529 /* Formula 16, p 217 */
1530 (assume(2*v+2*u+1>0),true);
1531 true$
1532 (assume(2*v-2*u+1>0),true);
1533 true$
1534 specint(t^(v-1)*%w[k,u](a*t)*exp(-p*t),t);
1535 a^(u+1/2)*(p+a/2)^(-v-u-1/2)*gamma(v-u+1/2)*gamma(v+u+1/2)
1536                            *%f[2,1]([v+u+1/2,u-k+1/2],[v-k+1],(p-a/2)/(p+a/2))
1537          /gamma(v-k+1)$
1538
1539 /*
1540  * Bug 1376865
1541  * Formula 24, p. 178.
1542  * The symbol %ei has gone. The new function is expintegral_ei.
1543  */
1544 specint(expintegral_ei(t)*exp(-p*t),t);
1545 -log(p-1)/p;
1546
1547 /*
1548  * Formula 25, p. 178
1549  *
1550  * t^(-1/2)*expintegral_ei(-t)
1551  * -> -2*%pi/sqrt(p)*log(sqrt(p)+sqrt(p+1))
1552  * This is equivalent to:
1553  * -> -2*%pi/sqrt(p)*asinh(sqrt(p))
1554  *
1555  * The result of Maxima is equivalent to
1556  * -> -2*sqrt(%pi)/sqrt(p)*asinh(sqrt(p))
1557  *
1558  * This differs by a factor sqrt(%pi) from the result above.
1559  * I (DK) think Maxima is correct.
1560  *
1561  * Maxima will ask for the sign of 2*psey-1.
1562  * We add this fact to the database and get the a result.
1563  */
1564 (assume(2*?psey-1>0),done);
1565 done;
1566 ratsimp(specint(t^(-1/2)*expintegral_ei(-t)*exp(-p*t),t));
1567 -sqrt(%pi)*(log(sqrt(p+1)+sqrt(p))-log(sqrt(p+1)-sqrt(p)))/sqrt(p);
1568 (forget(2*?psey-1>0),done);
1569 done;
1570
1571 /*
1572  * Formula 26, p. 178
1573  *
1574  * sin(a*t)*%ei(-t)
1575  * -> -(p^2+a^2)^(-1)*(a/2*log((p+1)^2+a^2) - p*atan(a/(p+1)))
1576  */
1577 ratsimp(rectform(specint(sin(a*t)*expintegral_ei(-t)*exp(-p*t),t)));
1578 (2*p*atan(a/(p+1))-a*log(p^2+2*p+a^2+1))/(2*p^2+2*a^2);
1579
1580 /*
1581  * Formula 27, p. 178
1582  *
1583  * cos(a*t)*%ei(-t)
1584  * -> -(p^2+a^2)^(-1)*(p/2*log((p+1)^2+a^2) + a*atan(a/(p+1)))
1585  */
1586 ratsimp(rectform(specint(cos(a*t)*expintegral_ei(-t)*exp(-p*t),t)));
1587 -((2*a*atan(a/(p+1))+p*log(p^2+2*p+a^2+1))/(2*p^2+2*a^2));
1588
1589 /*
1590  * Formula 34, p. 179
1591  *
1592  * exp(b*t)*gamma_incomplete_lower(v,a*t)
1593  * -> gamma(v)*a^v*(p-b)^(-1)*(p+a-b)^(-v)
1594  */
1595 (assume(p>b,v>-1),0);
1596 0;
1597 radcan(specint(exp(b*t)*gamma_incomplete_lower(v,a*t)*exp(-p*t),t));
1598 a^v*gamma(v+1)/((p-b)*(p-b+a)^v*v);
1599
1600 /*
1601  * Formula 1, p. 214
1602  *
1603  * kbateman[0](t) -> 1/(p+1)
1604  *
1605  */
1606 radcan(specint(kbateman[0](t)*exp(-p*t),t));
1607 1/(p+1);
1608
1609 /*
1610  * More tests, based on special values from
1611  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/
1612  *
1613  * These tests mostly provided by Edmond Orignac
1614  */
1615 ratsimp(hgfred([-1/2,1/2,2],[1,3/2],x)-
1616                (3/4*sqrt(1-x)+1/4*asin(sqrt(x))/sqrt(x)));
1617 0;
1618
1619 ratsimp(hgfred([-1/2,1/2,2],[1,5/2],x)-
1620                (3*(sqrt(x)*sqrt(1-x)*(6*x-1)+(4*x+1)*asin(sqrt(x)))/(32*x**(3/2))));
1621 0;
1622
1623 ratsimp(hgfred([-1/2,1/2,5/2],[3/2,3/2],x)-
1624                (1/3*asin(sqrt(x))/sqrt(x)+2/3*sqrt(1-x)));
1625 0;
1626
1627 ratsimp(hgfred([-1/2,1,3/2],[1/2,2],x)-
1628                (2*((1+2*x)*sqrt(1-x)-1)/(3*x)));
1629 0;
1630
1631 ratsimp(factor(hgfred([-1/2,1,3/2],[1/2,3],x))-
1632                (4*(-2*(2*x+3)*(1-x)**(3/2)-5*x+6)/(15*x**2)));
1633 0;
1634
1635 ratsimp(factor(hgfred([-1/2,1,5/2],[1/2,1/2],x))-
1636         1/3*((8*x^2-13*x+3)/(1-x)^2-8*sqrt(x)*atanh(sqrt(x))));
1637 0;
1638
1639 ratsimp(factor(hgfred([-1/2,1,5/2],[1/2,3/2],x))-
1640         1/3*((3-4*x)/(1-x)-4*sqrt(x)*atanh(sqrt(x))));
1641 0;
1642
1643 ratsimp(factor(hgfred([-1/2,1,5/2],[3/2,3/2],x))-
1644                1/3*((1-2*x)*atanh(sqrt(x))/sqrt(x)+2));
1645 0;
1646
1647 ratsimp(factor(hgfred([-1/2,1,3],[2,2],x))-
1648                ((2-(2-5*x)*sqrt(1-x))/(6*x)));
1649 0;
1650
1651
1652 ratsimp(hgfred([-1/2,3/2,3/2],[1/2,5/2],x)-
1653                (-3/4*asin(sqrt(x))/x**(3/2)+3*(1+x-2*x**2)/(4*sqrt(1-x)*x)));
1654 0;
1655
1656
1657 ratsimp(hgfred([-1/2,3/2,3],[5/2,2],x)-
1658                ((3*sqrt(x)*sqrt(1-x)*(10*x-1)+3*asin(sqrt(x)))/(32*x**(3/2))));
1659 0;
1660
1661 /*
1662  * Tests for 2F1.  From
1663  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/06/01/
1664  */
1665
1666 /* Fixed a, b, z */
1667
1668 hgfred([a,b],[a],z);
1669 (1-z)^(-b);
1670
1671 hgfred([a,b],[b],z);
1672 (1-z)^(-a);
1673
1674 ratsimp(hgfred([a,b],[b-1],z)
1675         -((1-z)^(-a-1)*(b+(a-b+1)*z-1)/(b-1)));
1676 0;
1677
1678 /*
1679  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/06/04/
1680  */
1681
1682 ratsimp(hgfred([a,a/2+1],[a/2],z)
1683         - (1-z)^(-a-1)*(z+1));
1684 0;
1685
1686 hgfred([a,a+1/2],[2*a],z);
1687 2^(2*a-1)*(sqrt(1-z)+1)^(1-2*a)/sqrt(1-z);
1688
1689
1690 hgfred([a,a+1/2],[2*a+1],z);
1691 2^(2*a)*(sqrt(1-z)+1)^(-2*a);
1692
1693 hgfred([a,a+1/2],[1/2],zp);
1694 ((1+sqrt(zp))^(-2*a)+(1-sqrt(zp))^(-2*a))/2;
1695
1696 ratsimp(hgfred([a,a+1/2],[3/2],z)
1697               -((1-sqrt(z))^(1-2*a)-(sqrt(z)+1)^(1-2*a))/(2*(2*a-1)*sqrt(z)));
1698 0;
1699
1700 hgfred([a,-a],[1/2],zp);
1701 cos(2*a*asin(sqrt(zp)));
1702
1703 hgfred([a,1-a],[1/2],zp);
1704 cos((2*a-1)*asin(sqrt(zp)))/sqrt(1-zp);
1705
1706 hgfred([a,1-a],[1],x1);
1707 legendre_p(a-1,1-2*x1);
1708
1709 hgfred([a,1-a],[3/2],z);
1710 sin((1-2*a)*asin(sqrt(z)))/((1-2*a)*sqrt(z));
1711
1712 hgfred([a,2-a],[3/2],z);
1713 sin((2-2*a)*asin(sqrt(z)))/((2-2*a)*sqrt(1-z)*sqrt(z));
1714
1715 /*
1716  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/01/
1717  */
1718
1719 factor(hgfred([-1/2,-1/2],[1/2],zp));
1720 asin(sqrt(zp))*sqrt(zp) + sqrt(1 - zp);
1721
1722 ratsimp(factor(hgfred([-1/2,-1/2],[3/2],zp))-
1723 1/4*((2*zp+1)*asin(sqrt(zp))/sqrt(zp)+3*sqrt(1-zp)));
1724 0;
1725
1726 ratsimp(factor(hgfred([-1/2,-1/2],[5/2],zp))-
1727 3/64/zp^(3/2)*(sqrt(1-zp)*sqrt(zp)*(1+14*zp)+(8*zp^2+8*zp-1)*asin(sqrt(zp)))
1728 );
1729 0;
1730
1731 ratsimp(factor(hgfred([-1/2,-1/2],[7/2],zp))-
1732 5/768/zp^(5/2)*(sqrt(1-zp)*sqrt(zp)*(-3+16*zp+92*zp^2)+3*(16*zp^3+24*zp^2-6*zp+1)*asin(sqrt(zp)))
1733 );
1734 0;
1735
1736 /*
1737  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/02/
1738  */
1739
1740 ratsimp(factor(hgfred([-1/2,1/2],[3/2],zp))-
1741 (asin(sqrt(zp))+sqrt(zp)*sqrt(1-zp))/2/sqrt(zp)
1742 );
1743 0;
1744
1745 ratsimp(factor(hgfred([-1/2,1/2],[5/2],zp))-
1746 3/16/zp^(3/2)*(sqrt(1-zp)*sqrt(zp)*(1+2*zp)+(4*zp-1)*asin(sqrt(zp)))
1747 );
1748 0;
1749
1750 ratsimp(factor(hgfred([-1/2,1/2],[7/2],zp))-
1751 5/128/zp^(5/2)*(sqrt(1-zp)*sqrt(zp)*(3+2*zp)*(4*zp-1)+3*(8*zp^2-4*zp+1)*asin(sqrt(zp)))
1752 );
1753 0;
1754
1755 /*
1756  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/03/
1757  */
1758
1759 hgfred([-1/2,1],[1/2],zp);
1760 1-sqrt(zp)*atanh(sqrt(zp));
1761
1762 ratsimp(factor(hgfred([-1/2,1],[3/2],zp))-
1763 1/2*(1-(zp-1)/sqrt(zp)*atanh(sqrt(zp)))
1764 );
1765 0;
1766
1767 ratsimp(factor(hgfred([-1/2,1],[2],zp))-
1768 2/3/zp*(1-(1-zp)^(3/2))
1769 );
1770 0;
1771
1772 ratsimp(factor(hgfred([-1/2,1],[5/2],zp))-
1773 3/8/zp^(3/2)*(sqrt(zp)*(zp+1)-(zp-1)^2*atanh(sqrt(zp)))
1774 );
1775 0;
1776
1777 ratsimp(factor(hgfred([-1/2,1],[3],zp))-
1778 4/15/zp^2*(-2+5*zp+2*(1-zp)^(5/2))
1779 );
1780 0;
1781
1782 ratsimp(factor(hgfred([-1/2,1],[7/2],zp))-
1783 (-5/48/zp^(5/2)*(3*atanh(sqrt(zp))*(zp-1)^3+sqrt(zp)*(-3*zp^2-8*zp+3)))
1784 );
1785 0;
1786
1787 ratsimp(factor(hgfred([-1/2,1],[4],zp))-
1788 (-(2*(8*(1-zp)^(7/2)-35*zp^2+28*zp-8))/(35*zp^3))
1789 );
1790 0;
1791
1792 /*
1793  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/04/
1794  */
1795
1796 ratsimp(factor(hgfred([-1/2,3/2],[1/2],zp))-
1797 (1-2*zp)/sqrt(1-zp)
1798 );
1799 0;
1800
1801 ratsimp(factor(hgfred([-1/2,3/2],[5/2],zp))-
1802 3/8/zp^(3/2)*(sqrt(1-zp)*sqrt(zp)*(-1+2*zp)+asin(sqrt(zp)))
1803 );
1804 0;
1805
1806 ratsimp(factor(hgfred([-1/2,3/2],[7/2],zp))-
1807 5/32/zp^(5/2)*(sqrt(1-zp)*sqrt(zp)*(3-4*zp+4*zp^2)+(6*zp-3)*asin(sqrt(zp)))
1808 );
1809 0;
1810
1811 /*
1812  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/05/
1813  */
1814
1815 ratsimp(factor(hgfred([-1/2,2],[1/2],zp))-
1816 (3*zp-3*(zp-1)*atanh(sqrt(zp))*sqrt(zp)-2)/2/(zp-1)
1817 );
1818 0;
1819
1820 ratsimp(factor(hgfred([-1/2,2],[1],zp))-
1821 (2-3*zp)/2/sqrt(1-zp)
1822 );
1823 0;
1824
1825 ratsimp(factor(hgfred([-1/2,2],[3/2],zp))-
1826 1/4*((1-3*zp)/sqrt(zp)*atanh(sqrt(zp))+3)
1827 );
1828 0;
1829
1830 ratsimp(factor(hgfred([-1/2,2],[3],zp))-
1831 4/15/zp^2*(sqrt(1-zp)*(zp-1)*(2+3*zp)+2)
1832 );
1833 0;
1834
1835 ratsimp(factor(hgfred([-1/2,2],[5/2],zp))-
1836 1/16/zp^(3/2)*(3*sqrt(zp)*(3*zp-1)+(-9*zp^2+6*zp+3)*atanh(sqrt(zp)))
1837 );
1838 0;
1839
1840 ratsimp(factor(hgfred([-1/2,2],[4],zp))-
1841 8/35/zp^3*(-4+7*zp+sqrt(1-zp)*(zp-1)^2*(4+3*zp))
1842 );
1843 0;
1844
1845 ratsimp(factor(hgfred([-1/2,2],[7/2],zp))-
1846 (-5/32/zp^(5/2)*(3*(zp+1)*atanh(sqrt(zp))*(zp-1)^2+sqrt(zp)*(-3*zp^2+2*zp-3)))
1847 );
1848 0;
1849
1850 /*
1851  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/06/
1852  */
1853
1854 ratsimp(factor(hgfred([-1/2,5/2],[1/2],zp))-
1855 (8*zp^2-12*zp+3)/3/(1-zp)^(3/2)
1856 );
1857 0;
1858
1859 ratsimp(factor(hgfred([-1/2,5/2],[3/2],zp))-
1860 (3-4*zp)/3/sqrt(1-zp)
1861 );
1862 0;
1863
1864 ratsimp(factor(hgfred([-1/2,5/2],[7/2],zp))-
1865 5/48/zp^(5/2)*(sqrt(1-zp)*sqrt(zp)*(1+2*zp)*(-3+4*zp)+3*asin(sqrt(zp)))
1866 );
1867 0;
1868
1869 /*
1870  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/07/
1871  */
1872
1873 ratsimp(factor(hgfred([-1/2,3],[1/2],zp))-
1874 (-15*sqrt(zp)*atanh(sqrt(zp))*(zp-1)^2+15*zp^2-25*zp+8)/8/(zp-1)^2
1875 );
1876 0;
1877
1878 ratsimp(factor(hgfred([-1/2,3],[1],zp))-
1879 (15*zp^2-24*zp+8)/8/(1-zp)^(3/2)
1880 );
1881 0;
1882
1883 ratsimp(factor(hgfred([-1/2,3],[3/2],zp))-
1884 (sqrt(zp)*(15*zp-13)-3*(5*zp^2-6*zp+1)*atanh(sqrt(zp)))/16/(zp-1)/sqrt(zp)
1885 );
1886 0;
1887
1888 ratsimp(factor(hgfred([-1/2,3],[2],zp))-
1889 (4-5*zp)/4/sqrt(1-zp)
1890 );
1891 0;
1892
1893 ratsimp(factor(hgfred([-1/2,3],[5/2],zp))-
1894 3/64/zp^(3/2)*(sqrt(zp)*(15*zp-1)+(-15*zp^2+6*zp+1)*atanh(sqrt(zp)))
1895 );
1896 0;
1897
1898 ratsimp(factor(hgfred([-1/2,3],[7/2],zp))-
1899 (-5/128/zp^(5/2)*(sqrt(zp)*(-15*zp^2+4*zp+3)+3*(5*zp^3-3*zp^2-zp-1)*atanh(sqrt(zp))))
1900 );
1901 0;
1902
1903 ratsimp(factor(hgfred([-1/2,3],[4],zp))-
1904 2/35/zp^3*(8-(1-zp)^(3/2)*(8+12*zp+15*zp^2))
1905 );
1906 0;
1907
1908 /*
1909  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/08/
1910  */
1911
1912 ratsimp(factor(hgfred([-1/2,7/2],[1/2],zp))-
1913 (5-30*zp+40*zp^2-16*zp^3)/5/(1-zp)^(5/2)
1914 );
1915 0;
1916
1917 ratsimp(factor(hgfred([-1/2,7/2],[3/2],zp))-
1918 (24*zp^2-40*zp+15)/15/(1-zp)^(3/2)
1919 );
1920 0;
1921
1922 ratsimp(factor(hgfred([-1/2,7/2],[5/2],zp))-
1923 (5-6*zp)/5/sqrt(1-zp)
1924 );
1925 0;
1926
1927 /*
1928  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/09/
1929  */
1930
1931 ratsimp(factor(hgfred([-1/2,4],[1/2],zp))-
1932 1/48/(zp-1)^3*(-105*sqrt(zp)*atanh(sqrt(zp))*(zp-1)^3+105*zp^3-280*zp^2+231*zp-48)
1933 );
1934 0;
1935
1936 ratsimp(factor(hgfred([-1/2,4],[1],zp))-
1937 (16-72*zp+90*zp^2-35*zp^3)/16/(1-zp)^(5/2)
1938 );
1939 0;
1940
1941 ratsimp(factor(hgfred([-1/2,4],[3/2],zp))-
1942 (sqrt(zp)*(105*zp^2-190*zp+81)-15*(zp-1)^2*(7*zp-1)*atanh(sqrt(zp)))/96/(zp-1)^2/sqrt(zp)
1943 );
1944 0;
1945
1946 ratsimp(factor(hgfred([-1/2,4],[2],zp))-
1947 (35*zp^2-60*zp+24)/24/(1-zp)^(3/2)
1948 );
1949 0;
1950
1951 ratsimp(factor(hgfred([-1/2,4],[5/2],zp))-
1952 (sqrt(zp)*(105*zp^2-100*zp+3)-3*(35*zp^3-45*zp^2+9*zp+1)*atanh(sqrt(zp)))/128/(zp-1)/zp^(3/2)
1953 );
1954 0;
1955
1956 ratsimp(factor(hgfred([-1/2,4],[3],zp))-
1957 (6-7*zp)/6/sqrt(1-zp)
1958 );
1959 0;
1960
1961 ratsimp(factor(hgfred([-1/2,4],[7/2],zp))-
1962 5/768/zp^(5/2)*(sqrt(zp)*(105*zp^2-10*zp-3)-3*(35*zp^3-15*zp^2-3*zp-1)*atanh(sqrt(zp)))
1963 );
1964 0;
1965
1966 /*
1967  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/20/
1968  */
1969
1970 ratsimp(factor(hgfred([1/2,1],[-1/2],zp))-
1971 (1-3*zp)/(zp-1)^2
1972 );
1973 0;
1974
1975 hgfred([1/2,1],[3/2],zp);
1976 log((sqrt(zp)+1)/(1-sqrt(zp)))/(2*sqrt(zp));
1977
1978 ratsimp(factor(hgfred([1/2,1],[2],zp))-
1979 2*(1-sqrt(1-zp))/zp
1980 );
1981 0;
1982
1983 ratsimp(factor(hgfred([1/2,1],[5/2],zp))-
1984 3/2/zp^(3/2)*((zp-1)*atanh(sqrt(zp))+sqrt(zp))
1985 );
1986 0;
1987
1988 ratsimp(factor(hgfred([1/2,1],[3],zp))-
1989 (-4/3/zp^2*((2*sqrt(1-zp)-3)*zp-2*sqrt(1-zp)+2))
1990 );
1991 0;
1992
1993 ratsimp(factor(hgfred([1/2,1],[7/2],zp))-
1994 5/8/zp^(5/2)*(3*atanh(sqrt(zp))*(zp-1)^2+sqrt(zp)*(5*zp-3))
1995 );
1996 0;
1997
1998 ratsimp(factor(hgfred([1/2,1],[4],zp))-
1999 (2*(-8*(1-zp)^(5/2)+15*zp^2-20*zp+8))/5/zp^3
2000 );
2001 0;
2002
2003 /*
2004  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/22/
2005  */
2006
2007 factor(hgfred([1/2,2],[-1/2],z));
2008 (5*z-1)/(z-1)^3;
2009
2010 ratsimp(factor(hgfred([1/2,2],[1],zp))-
2011 (2-zp)/2/(1-zp)^(3/2)
2012 );
2013 0;
2014
2015 ratsimp(factor(hgfred([1/2,2],[3/2],zp))-
2016 1/2*(atanh(sqrt(zp))/sqrt(zp)+1/(1-zp))
2017 );
2018 0;
2019
2020 ratsimp(factor(hgfred([1/2,2],[5/2],zp))-
2021 ((3*(zp+1)/4/zp^(3/2)*atanh(sqrt(zp))-3/4/zp))
2022 );
2023 0;
2024
2025 ratsimp(factor(hgfred([1/2,2],[3],zp))-
2026 4/3/zp^2*(2-(zp+2)*sqrt(1-zp))
2027 );
2028 0;
2029
2030 ratsimp(factor(hgfred([1/2,2],[7/2],zp))-
2031 15/16/zp^(5/2)*((zp^2+2*zp-3)*atanh(sqrt(zp))-(zp-3)*sqrt(zp))
2032 );
2033 0;
2034
2035 ratsimp(factor(hgfred([1/2,2],[4],zp))-
2036 8/5/zp^3*(-sqrt(1-zp)*zp^2+(5-3*sqrt(1-zp))*zp+4*sqrt(1-zp)-4)
2037 );
2038 0;
2039
2040 /*
2041  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/03/07/24/
2042  */
2043
2044 ratsimp(hgfred([1/2,3],[-1/2],z) - (1-7*z)/(z-1)^4);
2045 0$
2046
2047 ratsimp(factor(hgfred([1/2,3],[1],z))-
2048 (3*z^2-8*z+8)/8/(1-z)^(5/2)
2049 );
2050 0;
2051
2052 ratsimp(factor(hgfred([1/2,3],[3/2],zp))-
2053 (3*atanh(sqrt(zp))*(zp-1)^2+sqrt(zp)*(5-3*zp))/8/(zp-1)^2/sqrt(zp)
2054 );
2055 0;
2056
2057 ratsimp(factor(hgfred([1/2,3],[2],zp))-
2058 (4-3*zp)/4/(1-zp)^(3/2)
2059 );
2060 0;
2061
2062 ratsimp(factor(hgfred([1/2,3],[5/2],zp))-
2063 3/16/(zp-1)/zp^(3/2)*(sqrt(zp)*(1-3*zp)+(3*zp^2-2*zp-1)*atanh(sqrt(zp)))
2064 );
2065 0;
2066
2067 ratsimp(factor(hgfred([1/2,3],[7/2],zp))-
2068 15/64/zp^(5/2)*((3*zp^2+2*zp+3)*atanh(sqrt(zp))-3*sqrt(zp)*(zp+1))
2069 );
2070 0;
2071
2072 /*
2073  * Tests for 2F2
2074  */
2075 hgfred([a,b],[a,c],z);
2076 %f[1,1]([b],[c],z)$
2077
2078 /* test with exponentials */
2079
2080 /* Part 1: exponential times polynomial */
2081
2082 ratsimp(factor(hgfred([3/2,3/2],[1/2,1/2],x))-
2083         exp(x)*(4*x^2+8*x+1));
2084 0;
2085
2086 ratsimp(factor(hgfred([3/2,2],[1/2,1],x))-
2087         exp(x)*(2*x^2+5*x+1));
2088 0;
2089
2090 ratsimp(factor(hgfred([3/2,3],[1/2,2],x))-
2091         exp(x)*(2*x^2+7*x+2)/2);
2092 0;
2093
2094 ratsimp(factor(hgfred([2,2],[1,1],x))-
2095         exp(x)*(x^2+3*x+1));
2096 0;
2097
2098 ratsimp(factor(hgfred([3,3],[2,2],x))-
2099         exp(x)*(x^2+5*x+4)/4);
2100 0;
2101
2102 ratsimp(factor(hgfred([3,3],[1,2],x))-
2103         exp(x)*(x^3+8*x^2+14*x+4)/4);
2104 0;
2105
2106 /* Part 2: exponential times rational fraction */
2107
2108 ratsimp(factor(exponentialize(expand(block([gamma_expand:true], hgfred([1,3/2],[1/2,2],x)))))-
2109         (exp(x)*(2*x-1)+1)/x);
2110 0;
2111
2112 ratsimp(factor(exponentialize(expand(block([gamma_expand:true], hgfred([1,5/2],[3/2,2],x)))))-
2113         (exp(x)*(2*x+1)-1)/(3*x));
2114 0;
2115
2116 ratsimp(factor(exponentialize(expand(block([gamma_expand:true], hgfred([1,3],[2,2],x)))))-
2117         (exp(x)*(x+1)-1)/(2*x));
2118 0;
2119
2120
2121 /* test with incomplete gamma */
2122 /*
2123 hgfred([1/2,b],[3/2,b+1],x);
2124 b*(sqrt(%pi/x)*erfi(sqrt(x))-(-x)^(-b)*(gamma(b)-gamma_incomplete(b,-z)))/(2*b-1)$
2125 */
2126
2127 /* Note: I am unsure about Maxima notation for the incomplete gamma
2128 function */
2129
2130 /* Tests with erfi
2131    erfi is implemented as a simplifying function. The following definition no
2132    longer works. When erf_representation is ERF the simplifier transform erfi to
2133    a representation in terms of erf.
2134 (erfi(x) := -%i*erf(%i*x),0);
2135 0;
2136 */
2137
2138 erf_representation : erf;
2139 erf;
2140
2141 ratsimp(factor(hgfred([-1/2,3/2],[1/2,1/2],x))-
2142         (exp(x)-2*sqrt(%pi)*sqrt(x)*erfi(sqrt(x))));
2143 0;
2144
2145
2146 ratsimp(factor(hgfred([-1/2,2],[1/2,1],x))-
2147         (exp(x)-3/2*sqrt(%pi)*sqrt(x)*erfi(sqrt(x))));
2148 0;
2149
2150
2151 ratsimp(factor(hgfred([-1/2,5/2],[1/2,3/2],x))-
2152         (exp(x)-4/3*sqrt(%pi)*sqrt(x)*erfi(sqrt(x))));
2153 0;
2154
2155 ratsimp(factor(hgfred([-1/2,3],[1/2,2],x))-
2156         (exp(x)-5/4*sqrt(%pi)*sqrt(x)*erfi(sqrt(x))));
2157 0;
2158
2159 erf_representation : false;
2160 false;
2161
2162 /*
2163 ratsimp(factor(hgfred([(-1)/2,1],[1/2,2],x))-
2164         (1/(3*x)*(exp(x)*(2*x+1)-2*sqrt(%pi)*erfi(sqrt(x))*x^(3/2)-1)));
2165 0;
2166 */
2167
2168 /*
2169 hgfred([(-1)/2,1],[1/2,3],x);
2170 1/(15*x^2)*(exp(x)*(8*x^2+4*x+2)-8*sqrt(%pi)*erfi(sqrt(x))*x^(5/2)-2*(5*x+3))$
2171 */
2172
2173 /*
2174 hgfred([(-1)/2,1],[3/2,2],x);
2175 1/(6*x)*(2*exp(x)*(x-1)+sqrt(%pi)*sqrt(x)*erfi(sqrt(x))*(3-2*x)+2)$
2176 */
2177
2178 /*
2179 hgfred([(-1)/2,1],[3/2,3],x);
2180 2/(15*x^2)*(exp(x)*(2*x^2-4*x-1)+sqrt(%pi)*x^(3/2)*erfi(sqrt(x))*(5-2*x)+5*z+1)$
2181 */
2182
2183 /* Tests with erf */
2184 ratsimp(factor(block([gamma_expand:true], hgfred([1,3],[1/2,2],x)))-
2185         (2*x+4+sqrt(%pi)*sqrt(x)*(2*x+5)*exp(x)*erf(sqrt(x)))/4);
2186 0;
2187
2188
2189
2190
2191 /* Tests with modified bessel  functions I_0 and I_1 */
2192
2193 /*
2194 hgfred([1,3/2],[2,2],x);
2195 2*(exp(x/2)*bessel_i(0,x/2)-1)/x$
2196 */
2197
2198 /*
2199 hgfred([1/2,5/2],[3/2,2],x);
2200 exp(x/2)*(3*bessel_i(0,x/2)-bessel_i(1,x/2))/3$
2201 */
2202
2203 /*
2204 hgfred([1/2,3],[2,2],x);
2205 exp(x/2)*(2*bessel_i(0,x/2)-bessel_i(1,x/2))/2$
2206 */
2207
2208 /*
2209 hgfred([-1/2,3/2],[1/2,1],x);
2210 exp(x/2)*((1-2*x)*bessel_i(0,x/2)+*x*bessel_i(1,x/2))$
2211 */
2212
2213 /*
2214 hgfred([(-1)/2,1],[2,2],x);
2215 1/(9*x)*(4*exp(x/2)*x*(x-2)*bessel_i(1,x/2)-2*exp(x/2)*(2*x^2-6*x+3)*bessel_i(0,x/2)+6)$
2216 */
2217
2218 /*
2219 hgfred([(-1)/2,1],[2,3],x);
2220 4/(45*x)*(exp(x/2)*(4*x^2-14*x+3)*bessel_i(1,x/2)-exp(x/2)*(4*x^2-18*x+15)*bessel_i(0,x/2)+15)$
2221 */
2222
2223
2224 /* Test with exponential integral */
2225
2226 /*
2227 hgfred([1,1],[2,2],x);
2228 (2*%ei(x)+log(1/x)-log(x)-2%gamma)/(2*x)$
2229 */
2230 /* The above formula may be simplificated by noting
2231 log(1/x)=-log(x) for x>0; the expression might be in terms of E1 instead
2232 of Ei */
2233
2234
2235 /*
2236  * Tests for 3F2
2237  */
2238
2239 /*
2240  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/05/
2241  */
2242
2243 factor(hgfred([-1/2,1/2,2],[1,3/2],zp));
2244 (3*sqrt(1-zp)*sqrt(zp)+asin(sqrt(zp)))/(4*sqrt(zp));
2245
2246 ratsimp(factor(hgfred([-1/2,1/2,2],[1,5/2],zp))-
2247 3*(sqrt(1-zp)*sqrt(zp)*(6*zp-1)+(4*zp+1)*asin(sqrt(zp)))/(32*zp^(3/2)));
2248 0;
2249
2250 /*
2251  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/06/
2252  */
2253
2254 factor(hgfred([-1/2,1/2,5/2],[3/2,3/2],zp));
2255 (2*sqrt(1-zp)*sqrt(zp)+asin(sqrt(zp)))/(3*sqrt(zp));
2256
2257 /*
2258  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/07/
2259  */
2260
2261 ratsimp(factor(hgfred([-1/2,1/2,3],[1,3/2],zp))
2262 -(sqrt(zp)*(13-15*zp)+3*sqrt(1-zp)*asin(sqrt(zp)))/(16*sqrt(zp)*sqrt(1-zp)));
2263 0;
2264
2265 ratsimp(factor(hgfred([-1/2,1/2,3],[1,5/2],zp))
2266 -(3*sqrt(1-zp)*sqrt(zp)*(30*zp-1)+3*(12*zp+1)*asin(sqrt(zp)))/(128*zp^(3/2)));
2267 0;
2268
2269 factor(hgfred([-1/2,1/2,3],[3/2,2],zp));
2270 (5*sqrt(1-zp)*sqrt(zp)+3*asin(sqrt(zp)))/(8*sqrt(zp));
2271
2272 ratsimp(factor(hgfred([-1/2,1/2,3],[5/2,2],zp))
2273 -3*(sqrt(1-zp)*sqrt(zp)*(10*zp+1)+(12*zp-1)*asin(sqrt(zp)))/(64*zp^(3/2)));
2274 0;
2275
2276 /*
2277  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/09/
2278  */
2279
2280 ratsimp(factor(hgfred([-1/2,1,3/2],[1/2,2],z))
2281 -2*((1+2*z)*sqrt(1-z)-1)/(3*z));
2282 0;
2283
2284 ratsimp(factor(hgfred([-1/2,1,3/2],[1/2,3],z))
2285 -4*(-2*(2*z+3)*(1-z)^(3/2)-5*z+6)/(15*z^2));
2286 0;
2287
2288 /*
2289  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/11/
2290  */
2291
2292 ratsimp(factor(hgfred([-1/2,1,5/2],[1/2,2],z))
2293 -2*(4*z*(1-2*z)-sqrt(1-z)+1)/(9*sqrt(1-z)*z));
2294 0;
2295
2296 ratsimp(factor(hgfred([-1/2,1,5/2],[1/2,3],z))
2297 -4*(-5*z+2*sqrt(1-z)*(8*z^2+4*z+3)-6)/(45*z^2));
2298
2299 0;
2300
2301 ratsimp(factor(hgfred([-1/2,1,5/2],[3/2,2],z))
2302 -2*(1-(1-4*z)*sqrt(1-z))/(9*z));
2303 0;
2304
2305 ratsimp(factor(hgfred([-1/2,1,5/2],[3/2,3],z))
2306 -4*(5*z+2-2*(4*z+1)*(1-z)^(3/2))/(45*z^2));
2307 0;
2308
2309 /*
2310  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/13/
2311  */
2312 ratsimp(factor(hgfred([-1/2,3/2,3/2],[1/2,1/2],z))
2313 -(4*z^2-6*z+1)/(1-z)^(3/2));
2314 0;
2315
2316 ratsimp(factor(hgfred([-1/2,3/2,3/2],[1/2,5/2],zp))
2317 -3*(sqrt(zp)*(-2*zp^2+zp+1)-sqrt(1-zp)*asin(sqrt(zp)))/(4*sqrt(1-zp)*zp^(3/2)));
2318 0;
2319
2320 /*
2321  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/14/
2322  */
2323 ratsimp(factor(hgfred([-1/2,3/2,2],[1/2,1],z))
2324 -(6*z^2-9*z+2)/(2*(1-z)^(3/2)));
2325 0;
2326
2327 ratsimp(factor(hgfred([-1/2,3/2,2],[1/2,3],z))
2328 -4*((2*z^2+z+2)*sqrt(1-z)-2)/(5*z^2));
2329 0;
2330
2331 ratsimp(factor(hgfred([-1/2,3/2,2],[1,5/2],zp))
2332 -3*(sqrt(zp)*sqrt(1-zp)*(6*zp+1)-asin(sqrt(zp)))/(16*zp^(3/2)));
2333 0;
2334
2335 /*
2336  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/15/
2337  */
2338
2339 ratsimp(factor(hgfred([-1/2,3/2,5/2],[1/2,1/2],z))
2340 -(-16*z^3+40*z^2-30*z+3)/(3*(1-z)^(5/2)));
2341 0;
2342
2343 /*
2344  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/16/
2345  */
2346
2347 ratsimp(factor(hgfred([-1/2,3/2,3],[1/2,1],z))
2348 -(-30*z^3+75*z^2-56*z+8)/(8*(1-z)^(5/2)));
2349 0;
2350
2351
2352 ratsimp(factor(hgfred([-1/2,3/2,3],[1/2,2],z))
2353 -(10*z^2-15*z+4)/(4*(1-z)^(3/2)));
2354 0;
2355
2356 /*
2357  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/17/
2358  */
2359
2360 ratsimp(factor(hgfred([-1/2,2,2],[1,1],z))
2361 -(9*z^2-14*z+4)/(4*(1-z)^(3/2)));
2362 0;
2363
2364 ratsimp(factor(hgfred([-1/2,2,2],[1,3],z))
2365 -2*((9*z^2+2*z+4)*sqrt(1-z)-4)/(15*z^2));
2366 0;
2367
2368 /*
2369  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/18/
2370  */
2371
2372 ratsimp(factor(hgfred([-1/2,2,5/2],[1/2,1],z))
2373 -(2-15*z+20*z^2-8*z^3)/(2*(1-z)^(5/2)));
2374 0;
2375
2376 ratsimp(factor(hgfred([-1/2,2,5/2],[1/2,3],z))
2377 -4*(-8*z^3+4*z^2+z+2*sqrt(1-z)-2)/(15*z^2*sqrt(1-z)));
2378 0;
2379
2380 ratsimp(factor(hgfred([-1/2,2,5/2],[1,3/2],z))
2381 -(12*z^2-19*z+6)/(6*(1-z)^(3/2)));
2382 0;
2383
2384 ratsimp(factor(hgfred([-1/2,2,5/2],[3,3/2],z))
2385 -4/(45*z^2)*((12*z^2+z+2)*sqrt(1-z)-2));
2386 0;
2387
2388 /*
2389  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/19/
2390  */
2391
2392 ratsimp(factor(hgfred([-1/2,2,3],[1,1],z))
2393 -(-45*z^3+114*z^2-88*z+16)/16/(1-z)^(5/2));
2394 0;
2395
2396 /*
2397  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/20/
2398  */
2399
2400 ratsimp(factor(hgfred([-1/2,5/2,5/2],[1/2,1/2],z))
2401 -(64*z^4-224*z^3+280*z^2-144*z+9)/9/(1-z)^(7/2));
2402 0;
2403
2404 ratsimp(factor(hgfred([-1/2,5/2,5/2],[1/2,3/2],z))
2405 -(-32*z^3+80*z^2-60*z+9)/9/(1-z)^(5/2));
2406 0;
2407
2408 ratsimp(factor(hgfred([-1/2,5/2,5/2],[3/2,3/2],z))
2409 -(16*z^2-26*z+9)/9/(1-z)^(3/2));
2410 0;
2411
2412 /*
2413  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/21/
2414  */
2415
2416 ratsimp(factor(hgfred([-1/2,5/2,3],[1/2,1],z))
2417 -(40*z^4-140*z^3+175*z^2-88*z+8)/8/(1-z)^(7/2));
2418 0;
2419
2420 ratsimp(factor(hgfred([-1/2,5/2,3],[1/2,2],z))
2421 -(-40*z^3+100*z^2-75*z+12)/12/(1-z)^(5/2));
2422 0;
2423
2424 ratsimp(factor(hgfred([-1/2,5/2,3],[3/2,1],z))
2425 -(-20*z^3+51*z^2-40*z+8)/8/(1-z)^(5/2));
2426 0;
2427
2428 ratsimp(factor(hgfred([-1/2,5/2,3],[3/2,2],z))
2429 -(20*z^2-33*z+12)/12/(1-z)^(3/2));
2430 0;
2431
2432 /*
2433  * http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric3F2/03/08/22/
2434  */
2435
2436 ratsimp(factor(hgfred([-1/2,3,3],[1,1],z))
2437 -(225*z^4-792*z^3+1000*z^2-512*z+64)/64/(1-z)^(7/2));
2438 0;
2439
2440 ratsimp(factor(hgfred([-1/2,3,3],[1,2],z))
2441 -(-75*z^3+192*z^2-152*z+32)/32/(1-z)^(5/2));
2442 0;
2443
2444 ratsimp(factor(hgfred([-1/2,3,3],[2,2],z))
2445 -(25*z^2-42*z+16)/16/(1-z)^(3/2));
2446 0;
2447
2448 /*******************************************************************************
2449  * Bug ID: 2847387 - hgfred([3/2,-b],[5/2],-1) bogus and
2450  * Bug ID: 2876277 - hgfred([3/2,-2],[5/2],-x) not fully simplified
2451  *
2452  * Some examples to show that we get correct and simple results.
2453  * Reference: wolframalpha.com
2454  */
2455
2456 hgfred([3/2,3],[5/2],-1);
2457 3*%pi/32;
2458
2459 hgfred([3/2,-2],[5/2],-1);
2460 92/35;
2461
2462 hgfred([3/2,-2],[5/2],1);
2463 8/35;
2464
2465 declare(a,noninteger);
2466 done;
2467
2468 ratsimp(hgfred([a,-2],[5/2],-1));
2469 1/35*(4*a^2+32*a+35);
2470
2471 ratsimp(hgfred([a,-2],[5/2],1));
2472 1/35*(4*a^2-24*a+35);
2473
2474 ratsimp(hgfred([-2,1],[1/2],z));
2475 (8*z^2-12*z+3)/3;
2476
2477 /******************************************************************************/
2478
2479 /* Bug 3495182: hgfred([-n/2, (n-1)/2],[1],x) --> error
2480  *
2481  * Just check that we don't get an error.
2482  */
2483
2484 hgfred([-n/2, (n-1)/2],[1],x);
2485 %f[2,1]([-n/2, n/2-1/2],[1],x);
2486
2487 /*
2488  * Bug 3578373: hgfred([-1/2,a+1],[a+2],x);
2489  *
2490  * Just check that we don't get an error about  0 to a negative exponent.
2491  */
2492 hgfred([-1/2,a+1],[a+2],x);
2493 %f[2,1]([-1/2, a+1], [a+2], x);