share/hypergeometric/rtest_hg.mac

   1 (kill(all),
   2    fpprec : 25,
   3    load("hypergeometric"),
   4    expand_hypergeometric : false,
   5    max_fpprec : 1000,
   6    start : absolute_real_time (),
   7    0);
   8 0$
   9 errcatch(max_fpprec : %pi);
  10 [false]$
  11
  12 max_fpprec;
  13 1000$
  14
  15 errcatch(expand_hypergeometric : 42);
  16 [false]$
  17
  18 expand_hypergeometric;
  19 false$
  20
  21 /* argument length checks */
  22
  23 errcatch(hypergeometric([3],[10]));
  24 []$
  25
  26 errcatch(hypergeometric());
  27 []$
  28
  29 errcatch(hypergeometric([]));
  30 []$
  31
  32 /* 0F0 simplification */
  33 hypergeometric([],[],x);
  34 exp(x)$
  35
  36 /* 1F0 simplification */
  37 hypergeometric([a],[],x);
  38 1/(1-x)^a$
  39
  40 hypergeometric([-5],[],x);
  41 (1-x)^5$
  42
  43 hypergeometric([%i],[],x);
  44 1/(1-x)^%i$
  45
  46 /* sort and remove parameters */
  47
  48 hypergeometric([a,a],[b],x);
  49 hypergeometric([a,a],[b],x)$
  50
  51 hypergeometric([a,a],[a],x);
  52 1/(1-x)^a$
  53
  54 hypergeometric([a,a],[a,b],x);
  55 hypergeometric([a],[b],x)$
  56
  57 hypergeometric([a,b],[c,d],0);
  58 1$
  59
  60 /* Special values */
  61 hypergeometric([5/3, 42],[%i, 23],0);
  62 1$
  63
  64 hypergeometric([5/3, 7.892],[%i, 23],0.0);
  65 1.0$
  66
  67 hypergeometric([5/3, 7.892],[%i, 23],0.0b0);
  68 1.0b0$
  69
  70 hypergeometric([],[],0.5);
  71 ''(exp(0.5))$
  72
  73 hypergeometric([],[],0.5b0);
  74 ''(exp(0.5b0))$
  75
  76 hypergeometric([],[],0.5 + %i);
  77 ''(exp(0.5 + %i))$
  78
  79 hypergeometric([],[],0.5b0 + %i);
  80 ''(exp(0.5b0 + %i))$
  81
  82 hypergeometric([0],[],1);
  83 1$
  84
  85 hypergeometric([0,a],[b],x);
  86 1$
  87
  88 /* undefined cases */
  89 hypergeometric([7],[-9],x);
  90 'hypergeometric([7],[-9],x)$
  91
  92 hypergeometric([inf],[-9],x);
  93 'hypergeometric([inf],[-9],x);
  94
  95 hypergeometric([7],[8,9,minf],0);
  96 'hypergeometric([7],[8,9,minf],0)$
  97
  98 /* 1F1 reflection */
  99 hypergeometric([a],[b],-x) - hypergeometric([b-a],[b],x)*%e^(-x);
 100 0$
 101
 102 hypergeometric([a],[b],p-q) - hypergeometric([b-a],[b],q-p)*%e^(p-q);
 103 0$
 104
 105 /* polynomial cases (include taylor cases)*/
 106 (fake_hg(a,b,x,n) := block([listarith : true],
 107     sum(product(xreduce("*",a + l) * x /((1 + l) * xreduce("*", b + l)),l,0,k-1), k, 0,n)), 0)$
 108 0$
 109
 110 (expand_hypergeometric : true,0);
 111 0$
 112
 113 ratsimp(hypergeometric([-5,a],[b],x) - fake_hg([-5,a], [b],x, 6));
 114 0$
 115
 116 ratsimp(hypergeometric([-5,3/5],[],x) - fake_hg([-5,3/5], [],x, 6));
 117 0$
 118
 119 ratsimp(hypergeometric([-5,3/5 + %i],[1/12],x) - fake_hg([-5,3/5 + %i], [1/12],x, 6));
 120 0$
 121
 122 ratsimp(hypergeometric([-5,3/5 + %i],[1/12], 4/3) - fake_hg([-5,3/5 + %i], [1/12], 4/3, 6));
 123 0$
 124
 125 (expand_hypergeometric : false,0);
 126 0$
 127
 128 ratsimp(hypergeometric([-5,3/5 + %i],[1/12], 4/3 + %i/9) - fake_hg([-5,3/5 + %i], [1/12], 4/3 + %i/9, 6));
 129 0$
 130
 131 hypergeometric([9,3],[1/2], taylor(x,x,0,3)) - fake_hg([9,3],[1/2], taylor(x,x,0,3),3);
 132 ''(taylor(0,x,0,3))$
 133
 134 hypergeometric([2,3,5],[8,8,74], taylor(x,x,0,4)) - fake_hg([2,3,5],[8,8,74], taylor(x,x,0,4),5);
 135 ''(taylor(0,x,0,4))$
 136
 137 /* Exp doesn't taylorize; this isn't a bug in hypergeometric */
 138 hypergeometric([],[], taylor(x,x,0,5)) - exp(x);
 139 ''(taylor(0,x,0,5))$
 140
 141 hypergeometric([3/8],[], taylor(x,x,0,5)) - 1/(1-x)^(3/8);
 142 ''(taylor(0,x,0,5))$
 143
 144 /* Chu-Vandermonde identity */
 145 sublis([a = 1/2, b = -2/3, c = 12, x = 1],
 146   hypergeometric([a,b],[c],x) -  gamma(c) * gamma(c - a - b) / (gamma(c-a) * gamma (c-b)));
 147 0$
 148
 149 sublis([a = 1/2, b = -2/3, c = 121/10, x = 1],
 150   hypergeometric([a,b],[c],x) -  gamma(c) * gamma(c - a - b) / (gamma(c-a) * gamma (c-b)));
 151 0$
 152
 153 /* numerical */
 154
 155 (fpprec : 20, close_p(a,b, dig) := is(abs(a-b) < 10^-dig * abs(a)),0);
 156 0$
 157
 158 (check_id(id, xo, dig) := close_p(rectform(sublis(xo, rhs(id))), rectform(sublis(xo, lhs(id))), dig),0);
 159 0$
 160
 161 close_p(hypergeometric([-30,1/2], [4,2 + %i], 1/2), hypergeometric([-30,1/2], [4,2 + %i],0.5b0), 20);
 162 true$
 163
 164 close_p(hypergeometric([-75,1/2], [1/2,2 + %i], 1/2),hypergeometric([-75,1/2], [1/2,2 + %i], 0.5b0), 17);
 165 true$
 166
 167 (id : hypergeometric([1/2,1/2],[3/2],z^2) = asin(z)/z,0);
 168 0$
 169
 170 check_id(id, [z= 0.1], 13);
 171 true$
 172
 173 check_id(id, [z=0.1 + 0.2 * %i], 13);
 174 true$
 175
 176 check_id(id, [z=-0.99b0], 18);
 177 true$
 178
 179 /* G&R 9.121, 2 */
 180 (id : hypergeometric([-n/2,-(n-1)/2],[1/2], z^2) = ((1 + z)^n + (1- z)^n)/2,0);
 181 0$
 182
 183 check_id(id, [z=0.42, n = 11], 15);
 184 true$
 185
 186 check_id(id, [z=0.42b0 * %i, n = 21], 15);
 187 true$
 188
 189 check_id(id, [z=0.42b0 * %i, n = 37], 15);
 190 true$
 191
 192 /* G&R 9.121, 6 */
 193
 194 (id : hypergeometric([1,1],[2],-z) = log(1 + z) / z,0);
 195 0$
 196
 197 check_id(id, [z = -0.99b0], 20);
 198 true$
 199
 200 check_id(id, [z = 0.99b0], 20);
 201 true$
 202
 203 check_id(id, [z = 0.056b0 - 0.1b0 * %i], 18);
 204 true$
 205
 206 check_id(id, [z = 0.042b0 - 0.1b0 * %i], 18);
 207 true$
 208
 209 (id : hypergeometric([100],[2],x) = hgfred([100],[2],x),0);
 210 0$
 211
 212 check_id(id,[x = 0.9],13);
 213 true$
 214
 215 abramowitz_id(abramowitz_id(hypergeometric([a,b],[c],x), "15.3.3"), "15.3.3");
 216 hypergeometric([a,b],[c],x)$
 217
 218 abramowitz_id(abramowitz_id(hypergeometric([6,2],[9],x), "15.3.3"), "15.3.3");
 219 hypergeometric([6,2],[9],x)$
 220
 221 /*
 222 hypergeom([2308094809027379/562949953421312,6980579422424269/562949953421312],
 223           [2871044762448691/562949953421312],-5685794529555251/562949953421312);
 224 */
 225
 226 close_p(hypergeometric([4.1, 12.4],[5.1], -10.1), float(1.8514804991703517693765213431502813b-7), 14);
 227 true$
 228
 229 /* expand */
 230
 231 block([expand_hypergeometric : false], hypergeometric([-4,5],[7],x));
 232 hypergeometric([-4,5],[7],x)$
 233
 234 block([expand_hypergeometric : true], hypergeometric([-4,5],[7],x));
 235 x^4/3-(5*x^3)/3+(45*x^2)/14-(20*x)/7+1$
 236
 237 hypergeometric([-4,5],[7],2/3);
 238 163/1701$
 239
 240 hypergeometric([-4,5],[7],2/(1 + %i));
 241 1/7-(5*%i)/21$
 242
 243 hypergeometric([-6],[2], rat(x));
 244 ''(rat((x^6-42*x^5+630*x^4-4200*x^3+12600*x^2-15120*x+5040)/5040))$
 245
 246 hypergeometric_simp(hypergeometric([2,1],[3],x));
 247 -2*(log(1-x)/x+1/x+(log(1-x)*(1-x))/x^2)$
 248
 249 hypergeometric_simp(hypergeometric([a,b],[c],x));
 250 hypergeometric([a,b],[c],x)$
 251
 252 /* conjugates */
 253
 254 conjugate(hypergeometric([a],[c],x + %i*y));
 255 hypergeometric([a],[c],x - %i*y)$
 256
 257 conjugate(hypergeometric([%i],[c,d],x + %i*y));
 258 hypergeometric([-%i],[c,d],x - %i*y)$
 259
 260 conjugate(hypergeometric([%i],[%i-5,d],x + %i*y));
 261 hypergeometric([-%i],[-%i-5,d],x - %i*y)$
 262
 263 /* Kummer transformation */
 264 hypergeometric([a],[b],-x);
 265 hypergeometric([b-a],[b],x)*%e^(-x)$
 266
 267 hypergeometric([b-a],[b],x)*%e^(-x) - hypergeometric([a],[b],-x);
 268 0$
 269
 270 /* derivatives */
 271
 272 diff(hypergeometric([a,b],[c],x),x);
 273 (a*hypergeometric([a+1,b+1],[c+1],x)*b)/c$
 274
 275 diff(hypergeometric([a,b],[],x),x);
 276 a*hypergeometric([a+1,b+1],[],x)*b$
 277
 278 errcatch(diff(hypergeometric([a,b],[c],x),a));
 279 []$
 280
 281 errcatch(diff(hypergeometric([x,b],[c],x),x));
 282 []$
 283
 284 rationalize(hypergeometric([3,5],[12],1.0));
 285 ''(limit(hgfred([3,5],[12],x),x,1))$
 286
 287 /* http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=Hypergeometric2F1 */
 288
 289 close_p(hypergeometric([10,1/2],[4/5],11.34), 1.3549350983566885E-11 - 0.0378552663235813 * %i, 14);
 290 true$
 291
 292 hypergeometric([-10, 4.1],[%i], -10.1 + %i);
 293 1.045909655864634E13 + 3.8059389422734617E13 * %i$
 294
 295 hypergeometric([3.4, 5.6],[0.9], 1.01);
 296 5.5055743253266464E17 -1.7888695373728154E17 * %i$
 297
 298 (fpprec : 20, hypergeometric([2,4],[7.1b0],5.6b0 + %i));
 299 0.08245572273577160061b0 - 0.34771832247673618857b0*%i$
 300
 301 (fpprec : 20, hypergeometric([2,4],[7.1b0],5.6b0));
 302 0.26986105349383515333b0 + 0.35425557684972100084b0*%i$
 303
 304 z : float(hypergeometric([-100,1/2],[4/5 + %i], 13/10));
 305 0.05884530173873348*%i+0.062826546005002537$
 306
 307 hypergeometric([-100,1/2],[4/5 + %i], 1.3);
 308 0.05884530173873348*%i+0.062826546005002537$
 309
 310 z : float(hypergeometric([-10, 41/10],[%i], -101/10 + %i));
 311 1.0459096558646348E13 + %i * 3.8059389422734648E13$
 312
 313 hypergeometric([-10, 4.1],[%i], -10.1 + %i);
 314 1.0459096558646348E13 + %i * 3.8059389422734648E13$
 315
 316 abramowitz_id(hypergeometric([a,b],[c],x),"15.3.4");
 317 hypergeometric([a,c-b],[c],x/(x-1))/(1-x)^a$
 318
 319 /* A&S 15.3.15 I think this is an identity off [1,inf) */
 320
 321 (remvalue(z), domain : complex, fpprec : 20, 0);
 322 0$
 323
 324 (id : hypergeometric([a,b],[2*b], x) = (1-x)^(-a/2) * hypergeometric([a/2, b-a/2],[b+1/2], x^2/(4*x - 4)),0);
 325 0$
 326
 327 check_id(id, [a = 1.56, b = 12.8, x = 1.067 + 0.001 * %i], 13);
 328 true$
 329
 330 check_id(id, [a = 1.56, b = 12.8, x = -100.42 - %i], 13);
 331 true$
 332
 333 check_id(id, [a = 1.56, b = 12.8, x = cos(42.0) + %i * sin(42.0)], 15);
 334 true$
 335
 336 check_id(id, [a = 1.56b0, b = 12.8b0, x = 1.067b0 + 0.001b0 * %i], 18);
 337 true$
 338
 339 check_id(id, [a = 1.56b0, b = 12.8b0, x = -100.42b0 - %i], 10);
 340 true$
 341
 342 check_id(id, [a = 1.56b0, b = 12.8b0, x = cos(42.0b0) + %i * sin(42.0b0)], 10);
 343 true$
 344
 345 /* A&S 15.3.16 */
 346
 347 (id : hypergeometric([a,b],[2*b],x)=(1-x/2)^(-a) * hypergeometric([a/2, 1/2 + a/2],[b + 1/2], x^2/(2-x)^2),0);
 348 0$
 349
 350 check_id(id, [a = -1.056b0, b = -1.8b0, x = -10.66b0 - %i], 18);
 351 true$
 352
 353 check_id(id, [a = -1.056b0, b = -1.8b0, x = cos(2.1b0) + %i * sin(2.1b0)], 18);
 354 true$
 355
 356 /* A&S 15.3.31 */
 357
 358 (id : hypergeometric([a,1-a],[c],x) = (1-x)^(c-1) * hypergeometric([(c-a)/2, (c + a -1)/2],[c],4*x-4*x^2),0);
 359 0$
 360
 361 check_id(id, [a = 0.056b0, c = -0.8b0, x = -101.7b0 + %i], 18);
 362 true$
 363
 364 check_id(id, [a = 0.056b0, c = -0.8b0, x = cos(2.1b0) + %i * sin(2.1b0)], 18);
 365 true$
 366
 367 /* A&S 15.2.19 */
 368
 369 (id : (b-a)*(1-x) * hypergeometric([a,b],[c],x) - (c-a) * hypergeometric([a-1,b],[c],x) +
 370      (c-b) * hypergeometric([a,b-1],[c],x),0);
 371 0$
 372
 373 close_p(hypergeometric([1,1],[2],0.998b0), -log(1 - 0.998b0)/0.998b0, fpprec-1);
 374 true$
 375
 376 close_p(hypergeometric([1,1],[2],0.998), (-log(1 - 0.998)/0.998), 15);
 377 true$
 378
 379 /* Maple uses base ten floats, so comparing Maple to Maxima is tricky. To
 380 compare a Maple result to Maxima, do:
 381
 382  (1) In Maxima, set fpprec and use rationalize to convert floats to rationals. For
 383      big floats, the value of rationalize depends on fpprec, so be sure to set fpprec
 384      to the desired value.
 385
 386  (2) In Maple, set Digits. Evaluate convert(hypergeom(...),float), where the arguments are
 387      rationals determined in step 1.
 388
 389 Don't get me wrong: Maple does allow IEEE double floats, but Maple big floats (that is
 390 non-doubles) are base ten, I think.
 391 */
 392
 393 close_p(hypergeometric([1,1],[2], 0.999b0), -log(1-0.999b0)/0.999b0, fpprec-1);
 394 true$
 395
 396 close_p(hypergeometric([1,1],[2],1.3), -log(1-1.3)/1.3, 15);
 397 true$
 398
 399 close_p(hypergeometric([1,1],[2],1.3b0), -log(1-1.3b0)/1.3b0, fpprec-1);
 400 true$
 401
 402 hypergeometric([100],[2],0.9b0);
 403 ''(nfloat(hgfred([100],[2],x),[x = 0.9b0], fpprec))$
 404
 405 is(abs(nfloat(id, [a = 2/3, b= 42, c = 27/7, x=-1776.0], 18)) < 1.0e-14);
 406 true$
 407
 408 /* actual regresssions */
 409
 410 diff(hypergeometric([a,b],[c,d],-x),x);
 411 -(a*hypergeometric([a+1,b+1],[c+1,d+1],-x)*b)/(c*d)$
 412
 413 close_p(nfloat(hypergeometric([3.4],[7.8],'z), ['z = 5.67]), 18.61833264303484, 15);
 414 true$
 415
 416 block([max_fpprec : 2000], hypergeometric([101],[11], -1181.0b0));
 417 ''(bfloat(sublis([x=-1181], hgfred([101],[11],x))))$
 418
 419 block([expand_hypergeometric : true], hypergeometric([-2],[2/3 + %i], x));
 420 (9/442-189*%i/442)*x^2+(18*%i/13-12/13)*x+1$
 421
 422 /* Using doubles, this test causes an overflow. The code should switch to big floats.
 423 Yes, the test is slow :(. */
 424
 425 hypergeometric([3],[8.7],-1000.1);
 426 2.8986069810127266E-7$
 427
 428 block([fpprec : 25], hypergeometric([-1000],[2/3],-1181.0b0));
 429 1.813981492619807238503838b731$
 430
 431 /* A&S 13.4.1 */
 432 (id : (b-a) * hypergeometric([a-1],[b],x) + (2*a-b+x)*hypergeometric([a],[b],x)
 433     - a * hypergeometric([a+1],[b],x),0);
 434 0$
 435
 436 is(abs(nfloat(id, [a=1.67, b = -56.5, x = -0.42], 15)) < 1.0e-14);
 437 true$
 438
 439 is(abs(nfloat(id, [a=1.67 + %i, b = -56.5, x = 10.42 + 11.88*%i], 15)) < 5.0e-14);
 440 true$
 441
 442 close_p(subst([a=1.67, b=56.1, x = -1.42], hypergeometric([a], [b], x) * (x - b + 2 * a)), -51.96358410181497,13);
 443 true$
 444
 445 close_p(subst([a=1.67, b=56.1, x = -1.42], hypergeometric([a], [b], x) * (x - b)), -55.16695012064225,13);
 446 true$
 447
 448 /* A&S 15.1.7 */
 449
 450 (id : hypergeometric([1/2,1/2],[3/2],-x^2) - log(x + sqrt(1+x^2)) / x,0);
 451 0$
 452
 453 is(abs(nfloat(id, [x = cos(-47.0) + %i * sin(-47.0)], 12)) < 10^-10);
 454 true$
 455
 456 /* A&S 15.1.4 */
 457
 458 (id : hypergeometric([1/2,1],[3/2],x^2) - log((1+x)/(1-x))/(2*x),0);
 459 0$
 460
 461 is(abs(nfloat(id, [x = exp(5.0 - 2.807 * %i)], 12)) < 10^-10);
 462 true$
 463
 464 is(abs(nfloat(id, [x = -10.7 + %i], 20)) < 10^-19);
 465 true$
 466
 467 is(abs(nfloat(id, [x = 28.07 + 0.01 * %i], 20)) < 10^-19);
 468 true$
 469
 470
 471 close_p(nfloat(hypergeometric([1/2,1],[3/2], x), [x =cos(23.0) + %i * sin(23.0)], 18),
 472   bfloat(rectform(hgfred([1/2,1],[3/2], cos(23) + %i * sin(23)))), 15);
 473 true$
 474
 475
 476 /* bugs!!!
 477
 478 is(abs(nfloat(id, [x = cos(-47.0b0) + %i * sin(-47.0b0)], 10)) < 10^-8);
 479 true$
 480
 481 hypergeometric([1/2,1/2],[3/2],10.0e6);
 482 hypergeometric([1/2,1/2],[3/2],10.0e2);
 483 hypergeometric([5.0b-1, 1.0b0],[1.5b0], 0.995b0);
 484 is(abs(nfloat(id, [x = -1964.1b0], 20)) < 10^-20);
 485 */
 486
 487 hypergeometric([1/2,1],[3/2], 2807.0);
 488 ''(float(3.562945408378278793890989970521468894557558424951846558009b-4
 489             - 2.964822314792439520081705392154987108698281707824080659881b-2 * %i))$
 490
 491 (load("functs"),0);
 492 0$
 493
 494 /* A&S 13.1.20 */
 495 (id : determinant(wronskian([hypergeometric([a],[b],x), x^(1-b) * hypergeometric([1+a-b],[2-b],x)],x)) - (1-b)*x^(-b)*exp(x),0);
 496 0$
 497
 498 is(abs(nfloat(id, [a = 1.45, b=2.807, x=3.14159], 18)) < 10^-18);
 499 true$
 500
 501 is(abs(nfloat(id, [a = 1.45 - %i, b=-2.807, x=3.14159 + 3.14159*%i], 18)) < 10^-18);
 502 true$
 503
 504 is(abs(nfloat(id, [a = 1.45 - %i, b=-2.807, x=-19.13], 18)) < 10^-18);
 505 true$
 506
 507 /* A&S 13.4.4 */
 508 (id : b * hypergeometric([a],[b],x) -b* hypergeometric([a-1],[b],x)-x*hypergeometric([a],[b+1],x),0);
 509 0$
 510
 511 is(abs(nfloat(id, [a = 10.45 - %i, b=-12.807, x=-29.13], 18)) < 10^-18);
 512 true$
 513
 514 /* A&S 13.6.1 */
 515 (id : hypergeometric([nu+1/2],[2*nu+1],2*%i*x) - gamma(1+nu) * exp(%i*x) * (x/2)^(-nu) * bessel_j(nu,x),0);
 516 0$
 517
 518 is(abs(rectform(sublis([nu=6.7, x = 28.07],id))) < 1/10^18);
 519 true$
 520
 521 is(abs(float(rectform(sublis([nu=-6.01, x = -2.0 + 4.0 * %i ],id)))) < 1/10^14);
 522 true$
 523
 524 f : abramowitz_id(42*elliptic_kc(x*(x-1)), "a&s17.3.9-->");
 525 21*%pi*hypergeometric([1/2,1/2],[1],(x-1)*x)$
 526
 527 abramowitz_id(f, "a&s17.3.9<--");
 528 42*elliptic_kc(x*(x-1))$
 529
 530 f : abramowitz_id(42*elliptic_ec(x*(x-1)), "a&s17.3.10-->");
 531 21*%pi*hypergeometric([-1/2,1/2],[1],(x-1)*x)$
 532
 533 abramowitz_id(f, "a&s17.3.10<--");
 534 42*elliptic_ec((x-1)*x)$
 535
 536 f : abramowitz_id(elliptic_kc(p) + elliptic_ec(q), "a&s17.3.9-->", "a&s17.3.10-->");
 537 %pi*hypergeometric([1/2,1/2],[1],p)/2+%pi*hypergeometric([-1/2,1/2],[1],q)/2$
 538
 539 abramowitz_id(f,"a&s17.3.9<--", "a&s17.3.10<--");
 540 elliptic_kc(p) + elliptic_ec(q)$
 541
 542 (id :  hypergeometric([3/2 + n/5, 2/3, -n, 2*n + 2],[n + 11/6, 4/3, n/5+1/2], 2/27) -
 543          pochhammer(5/2,n) * pochhammer(11/6,n) / (pochhammer(3/2,n) * pochhammer(7/2,n)),0);
 544 0$
 545
 546 makelist(sublis([n = k], id), k,0,10);
 547 ''(makelist(0,k,0,10))$
 548
 549 /* Value at 1/2; see http://en.wikipedia.org/wiki/List_of_hypergeometric_identities */
 550
 551 (id : hypergeometric([a,1-a],[c], 1/2) - gamma(c/2) *  gamma(c/2 + 1/2)/(gamma(c/2 + a/2) *
 552 gamma(1/2 + c/2 - a/2)),0);
 553 0$
 554
 555 is(abs(nfloat(id,[a=4/5, c = 2/3], 23)) < 1/10^23);
 556 true$
 557
 558 is(abs(nfloat(id,[a=28.07, c = 1909], 23)) < 1/10^23);
 559 true$
 560
 561 is(abs(nfloat(id,[a=-28.07 + %i, c = 19.09 - %i], 23)) < 1/10^18);
 562 true$
 563
 564 is(abs(nfloat(hypergeometric([1/6,1/2],[2/3],125/128) - 4 * 2^(1/6) / 3,[], 22)) < 1/10^22);
 565 true$
 566
 567
 568 /* and just for fun */
 569 taylor(hypergeometric([a+1,b],[c,d],x)*(((b-a-1)*x)/((c-a-1)*(d-a-1))-(c*d-2*a*d-2*d-2*a*c-2*c+3*a^2+7*a+4)/((c-a-1)*(d-a-1)))+hypergeometric([a+2,b],[c,d],x)*(((a+1)*x)/((c-a-1)*(d-a-1))-((a+1)*(d+c-3*a-5))/((c-a-1)*(d-a-1)))-((a+1)*(a+2)*hypergeometric([a+3,b],[c,d],x))/((c-a-1)*(d-a-1))+
 570 hypergeometric([a,b],[c,d],x),x,0,15);
 571 ''(taylor(0,x,0,15))$
 572
 573 /* SF bug 2855935  */
 574
 575 hypergeometric([a],[],1);
 576 hypergeometric([a],[],1)$
 577
 578 (assume(a < 0),0);
 579 0$
 580
 581 hypergeometric([a],[],1);
 582 0$
 583
 584 (forget(a < 0),0);
 585 0$
 586
 587 /* Integrals of hypergeometric functions */
 588 integrate(hypergeometric([1],[2],z),z);
 589 hypergeometric([1,1],[2,2],z)*z$
 590
 591 integrate(hypergeometric([a,b],[c],z),z);
 592 hypergeometric([a-1,b-1],[c-1],z)*(c-1)/((a-1)*(b-1))$
 593
 594 integrate(hypergeometric([a,b],[c],2*z),z);
 595 hypergeometric([a-1,b-1],[c-1],2*z)*(c-1)/(2*(a-1)*(b-1))$
 596
 597 /* next three tests regressions for Bug ID: 2949677 */
 598 hypergeometric([-4, 38/5],[13/5],1);
 599 3125/6279$
 600
 601 close_p(hypergeometric([-4, 38/5],[13/5],1.0), 3125/6279, 14);
 602 true$
 603
 604 close_p(hypergeometric([-4, 38/5],[13/5],1.0b0), 3125/6279, fpprec-2);
 605 true$
 606
 607 /* Regression for Kummer reflection - ID: 2954141 */
 608
 609 hypergeometric([-2],[3],x), expand_hypergeometric : true;
 610 x^2/12-(2*x)/3+1$
 611
 612 hypergeometric([-2],[3],-x),expand_hypergeometric : true;
 613 x^2/12+2*x/3+1$
 614
 615 hypergeometric([-1],[-2],x), expand_hypergeometric : true;
 616 x/2+1$
 617
 618
 619 /* regression: Unevaluated numerical hypergeometric() call - ID: 2996065 */
 620
 621 close_p(hypergeometric([-1/4,1/4],[1],16.0/25),0.95126303911614,15);
 622 true$
 623
 624 close_p(hypergeometric([1,1],[2], 0.998),6.227062222867928,15);
 625 true$
 626
 627 close_p(hypergeometric([1/2,-1/2],[1], 0.64), 0.812549610282201,15);
 628 true$
 629
 630 /* tests for new code for x near exp(%i * %pi / 3) */
 631
 632 close_p(hypergeometric([1/2,2/3],[6], 0.5b0+%i * sqrt(3.0b0)/2),
 633   5.56090808176771804282797055039701904767032b-2*%i
 634   +1.01996502902735430283851818152075979159011b0, fpprec-2);
 635 true$
 636
 637 close_p(hypergeometric([1/2,2/3],[6], 0.5b0 - %i * sqrt(3.0b0)/2),
 638         -5.56090808176771804282797055039701904767032b-2*%i
 639                 +1.01996502902735430283851818152075979159011b0, fpprec-2);
 640 true$
 641
 642 /* tests aimed toward increasing code coverage*/
 643
 644 (remvalue(id,a,b,c,d,x,z,f),0);
 645 0$
 646
 647 conjugate(hypergeometric([1/2, 6],[9], 1/2 + %i/8));
 648 hypergeometric([1/2, 6],[9], 1/2 - %i/8)$
 649
 650 conjugate(hypergeometric([a,b,c,d],[f], 8 + %i));
 651 '(conjugate(hypergeometric([a,b,c,d],[f], 8 + %i)))$
 652
 653 block([fpprec : 20], close_p(hypergeometric([5/8],[], 1.97b0), nfloat(1/(1-x)^(5/8),[x = 1.97b0],20), 18));
 654 true$
 655
 656 close_p(hypergeometric([-1],[3/4],0.5), 1/3, 15);
 657 true$
 658
 659 block([expand_hypergeometric : true], hypergeometric([-2, 1/2 + %i/7],[%i/9],x));
 660 (29241/16072-108783*%i/16072)*x^2+(9*%i-18/7)*x+1$
 661
 662 tex(hypergeometric([2],[9],x),false);
 663 "$$F\\left( \\left. \\begin{array}{c}2\\\\9\\end{array} \\right |,x\\right)$$
 664 "$
 665
 666 tex(hypergeometric([a,b],[c],x),false);
 667 "$$F\\left( \\left. \\begin{array}{c}a,\\;b\\\\c\\end{array} \\right |,x
 668  \\right)$$
 669 "$
 670
 671 /* float hypergeometric error | fix - ID: 3426847 */
 672 block([fpprec : 56], close_p(hypergeometric([1],[3/2,3/2], -10000.0), -4.012030270249197E-4, 15));
 673 true$
 674
 675 /* Unevaluated numerical hypergeometric() call - ID: 2996065 */
 676 close_p(hypergeometric([-1/4,1/4],[1],16.0/25), 0.95126303911614, 15);
 677 true$
 678
 679 /* clisp isn't very fast => with is(time < 40) we got build failures on launchpad. */
 680 (time : (absolute_real_time ()- start), print(time), is(time < 400));
 681 true$
 682
 683 (remvalue(start), reset(domain, fpprec, expand_hypergeometric, max_fpprec),
 684  remfunction(fake_hg, close_p, check_id),0);
 685 0$