1module ghyper; % Generalized Hypergeometric Functions. 2 3% Author : Victor Adamchik, Byelorussian University Minsk, Byelorussia. 4 5% Redistribution and use in source and binary forms, with or without 6% modification, are permitted provided that the following conditions are met: 7% 8% * Redistributions of source code must retain the relevant copyright 9% notice, this list of conditions and the following disclaimer. 10% * Redistributions in binary form must reproduce the above copyright 11% notice, this list of conditions and the following disclaimer in the 12% documentation and/or other materials provided with the distribution. 13% 14% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 15% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 16% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 17% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 18% CONTRIBUTORS 19% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 20% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 21% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 22% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 23% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 24% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 25% POSSIBILITY OF SUCH DAMAGE. 26% 27 28 29% Major modifications by: Winfried Neun, ZIB Berlin. 30 31% Oct 22, 2001, hypergeometric({a,b},{0},z) returns 32% unevaluated (without error) as requested by Francis Wright 33 34 35% The next 2 declarations enable better checking of number of arguments 36% by simpiden 37 38flag('(hypergeometic), 'specfn); 39put('hypergeomtric, 'number!-of!-args, 3); 40 41 put('ghf,'simpfn,'simpghf)$ 42 43 symbolic procedure simpghf u; 44 if null cddr u then 45 rerror('specialf,125, 46 "WRONG NUMBER OF ARGUMENTS TO GHF-FUNCTION") 47 else 48 if or(not numberp car u,not numberp cadr u) then 49 rerror('specialf,126,"INVALID AS INTEGER") 50 else 51 begin scalar vv,v; 52 v:=redpar1(cddr u,car u); 53 vv:=redpar1(cdr v,cadr u); 54 if null cddr vv then return 55 ghfsq(list(car u,cadr u),listsq car v, 56 listsq car vv, simp cadr vv); 57 return rerror ('specialf,127, 58 "WRONG NUMBER OF ARGUMENTS TO GHF-FUNCTION"); 59 end$ 60 61 symbolic procedure ghfexit(a,b,z); 62 begin scalar aa,bb; 63 aa:= 'list . listprepsq a; 64 bb:= 'list . listprepsq b; 65 return mksqnew({'hypergeometric,aa,bb,prepsq z}); 66 end; 67 68%*********************************************************************** 69%* GHF as a polynomial * 70%*********************************************************************** 71 72symbolic procedure listmaxsq u; 73 % u - list of numbers of SQ. 74 % return - max value. 75 if null cdr u then car u else 76 if null caar u then car u else 77 if null caadr u then cadr u else 78 if caar u > caadr u or car u = cadr u then 79 listmaxsq((car u) . cddr u) else 80 listmaxsq((cadr u) . cddr u)$ 81 82 symbolic procedure ghfpolynomp (u,a); 83 begin scalar w1,w2; 84m1: 85 if null u then 86 if null w1 then <<u:=w2; return (nil . a) >> 87 else <<u:=listmaxsq(w1); 88 a:=u . append(delete(u,w1),w2); 89 return (t . a)>> 90 else 91 if parfool(car u) then (w1:=(car u) . w1) 92 else (w2:=(car u) . w2); 93 u:=cdr u; 94 goto m1; 95 end$ 96 97 symbolic procedure polynom(u,a,b,z); 98 % u - list of SQ. 99 begin scalar s; integer k; 100 if null caar(u) then return '(1 . 1) else 101 s := ghfpolynomp (b,a); 102 a := cdr s; 103 if car s then 104 if null caar a or greaterp(caar a,caar u) then 105 <<%rerror('special,124, 106 % "zero in the denominator of the GHF-function"); 107 b:=a; a:=u; 108 return ghfexit(a,b,z); 109 >> 110 else b:=a; 111 k:=1; s:=1 . 1; 112m: 113 s:=addsq(s,quotsq(multsq(multpochh(u,simp k),exptsq(z,k)), 114 multpochh(append(list('(1 . 1)),b),simp k))); 115 k:=k+1; 116 if greaterp(k,numr negsq(car u)) then return s else goto m; 117 end$ 118 119%*********************************************************************** 120%* Lowering of the order GHF * 121%*********************************************************************** 122 123% This assigns to a and b! 124 symbolic smacro procedure ghflowering1p; 125 begin scalar sa,sb,w1,w2; 126 sa:=a; sb:=b; 127 m1: if null b then << a:=sa; b:=sb; return nil 128 >>; 129 m2: if null a then << w2:= (car b) . w2; 130 b:=cdr b; 131 a:=sa; w1:=nil; 132 goto m1 133 >> 134 else 135 if numberp(prepsq subtrsq(car a,car b)) and 136 greaterp(car(subtrsq(car a,car b)),0) then 137 << 138 b:=car b . append(w2,cdr b); 139 a:=subtrsq(car a,car b) . append(w1,cdr a); 140 return t 141 >> else 142 << 143 w1:=(car a) . w1; 144 a:=cdr a; 145 goto m2 146 >>; 147 end$ 148 149 symbolic procedure lowering1(x,y,u,z); 150 % x -- (m . a). 151 % y -- (g . b). 152 addsq(ghfsq(u,append(list(subtrsq(addsq(car x,car y),'(1 . 1))), 153 cdr x), 154 append(list(car y),cdr y),z), 155 multsq(ghfsq(u,append(list(addsq(car x,car y)),specfn!-listplus(cdr x, 156 '(1 . 1))), 157 append(list(addsq(car y,'(1 . 1))), 158 specfn!-listplus(cdr y,'(1 . 1))),z), 159 quotsq(multsq(z,multlist(cdr x)), 160 multsq(car y,multlist(cdr y)))))$ 161 162% This assigns to a and b 163 symbolic smacro procedure ghflowering2p; 164 begin scalar sa,sb,w1,wa,fl; 165 if equal(z,'(1 . 1)) then return nil; 166 sa:=a; sb:=b; 167 m1: if null b then 168 << b:=sb; 169 if fl then a:=wa . sa else a:= sa; 170 return nil 171 >>; 172 m2: if null a then 173 << b:=cdr b; 174 a:=sa; w1:=nil; 175 goto m1 176 >> 177 else 178 if numberp(prepsq subtrsq(car b,car a)) and 179 lessp(car(subtrsq(car a,car b)),0) then 180 if fl then 181 if not equal(wa,car a) then 182 << b:=sb; 183 a:=list(wa,car a) . append(w1,cdr a); 184 return t 185 >> 186 else 187 << 188 w1:=(car a) . w1; 189 a:=cdr a; 190 goto m2 191 >> 192 else 193 << fl:=t; 194 sa:=append(w1,cdr a); 195 wa:=car a; 196 b:=cdr b; a:=sa; w1:=nil; 197 goto m1 198 >> 199 else 200 << w1:= (car a) .w1; 201 a:=cdr a; 202 goto m2 203 >>; 204 end$ 205 206 symbolic procedure lowering2(x,b,u,z); 207 % x -- (r s).(a). 208 subtrsq(multsq(ghfsq(u,append(list(caar x,addsq('(1 . 1),cadar x)), 209 cdr x),b,z), 210 quotsq(cadar x,subtrsq(cadar x,caar x))), 211 multsq(ghfsq(u,append(list(addsq('(1 . 1),caar x),cadar x), 212 cdr x),b,z), 213 quotsq(caar x,subtrsq(cadar x,caar x))))$ 214 215% This assigns to a and b 216 symbolic smacro procedure ghflowering3p; 217 %return a = (mmm . a1). 218 begin scalar sa,w,mmm; % MM used in SPDE as a global. 219 sa:=a; 220m1: if null a then << a:=sa; return nil >> 221 else 222 if not numberp(prepsq car a) then 223 <<w:= (car a) . w; a:=cdr a; goto m1 >>; 224 225 if member ('(1 . 1), a) then <<mmm := '(1 . 1); 226 a:= delete('(1 . 1),a)>> 227 else << mmm:= car a; a:=cdr a >>; % WN 2.2 94 228 229m2: if null a then 230 if listnumberp b then << a:=mmm . w; return t >> 231 else << a:=sa; return nil>> 232 else 233 if equal(car a,'(1 . 1)) then 234 <<a:=sa; return nil>> 235 else 236 <<w:=(car a) . w; 237 a:=cdr a; 238 goto m2 239 >>; 240 end$ 241 242 symbolic procedure listnumberp(v); 243 % v -- list of SQ. 244 % value is T if numberp exist in (v). 245 if null v then nil 246 else 247 if numberp(prepsq car v) then t 248 else listnumberp(cdr v)$ 249 250 symbolic procedure lowering3(a,b,u,z); 251 multsq(quotsq(multlist(difflist(b,'(1 . 1))),multsq(z,multlist( 252 difflist(cdr a,'(1 . 1))))), 253 subtrsq(ghfsq(u, (car a) . difflist(cdr a,'(1 . 1)), 254 difflist(b,'(1 . 1)),z), 255 ghfsq(u,append(list(subtrsq(car a,'(1 . 1))), 256 difflist(cdr a,'(1 . 1))),difflist(b,'(1 . 1)),z)))$ 257 258%*********************************************************************** 259%* GHFsq, main entry * 260%*********************************************************************** 261 262 symbolic procedure ghfsq(u,a,b,z); 263 % u -- (p q) PF. 264 % a,b -- lists of SQ. 265 % z -- SQ. 266 begin scalar c,aaa; 267 u:=redpar(a,b); 268 a:=car u;b:=cadr u;u:=list(length(a), length(b)); 269 if null numr z then return '(1 . 1) else 270 if listparfool(b,(nil .1)) and not listparfool(a,(nil . 1)) then 271 % return rerror('specialf,128, 272 %"zero in the denominator of the GHF-function") 273 return ghfexit(a,b,z) 274 else 275 aaa := ghfpolynomp(a,a); 276 a := cdr aaa; 277 if car aaa then return polynom(a,a,b,z) else 278 if ghflowering1p() then return lowering1(a,b,u,z) else 279 if ghflowering2p() then return lowering2(a,b,u,z) else 280 if ghflowering3p() then return lowering3(a,b,u,z) else 281 if car u = 0 and cadr u = 0 then return expdeg(simp 'e,z) else 282 if car u = 0 and cadr u = 1 then return ghf01(a,b,z) else 283 if car u = 1 and cadr u = 0 then 284 if z='(1 . 1) then return ghfexit(a,b,z) 285 else 286 return expdeg(subtrsq('(1 . 1),z),if null a then '(nil . 1) 287 else negsq(car a)) 288 else 289 if car u = 1 and cadr u = 1 then return ghf11(a,b,z) else 290 if car u = 1 and cadr u = 2 then return ghf12(a,b,z) else 291 if car u = 2 and cadr u = 1 then return ghf21(a,b,z) else 292 if car u = cadr u + 1 then 293 if (c:=ghfmid(a,b,z)) = 'fail 294 then return ghfexit(a,b,z) 295 else return c; 296 297 if car u <= cadr u then return ghfexit(a,b,z); 298 return rerror('specialf,131,"hypergeometric series diverges"); 299 end$ 300 301 302%*********************************************************************** 303% p = q+1 * 304%*********************************************************************** 305 306 symbolic procedure ghfmid(a,b,z); 307 begin scalar c; 308 c:= redpar(a,difflist(b, '(1 . 1))); 309 if length(cadr c) > 0 or length(car c) > 1 then return 'fail 310 else 311 return formulaformidcase(length(b), caar c, 312 subtrsq(car b,'(1 . 1)), z); 313 end$ 314 315 symbolic procedure formulaformidcase(p,b,a,z); 316 if not(p = 1) and b = '(1 . 1) and z = '(1 . 1) then 317 multsq(simpx1(prepsq(multsq('(-1 . 1),a)),p,1), 318 quotsq(dfpsisq(a,simp(p-1)),gamsq(simp p))) 319 else 320 if b = '(1 . 1) and z='(-1 . 1) then 321 quotsq(multsq(simpx1(prepsq(multsq('(-1 . 2),a)),p,1), 322 subtrsq(dfpsisq(multsq(a, '(1 . 2)),simp(p-1)), 323 dfpsisq(multsq(addsq('(1 . 1),a),'(1 . 2)), 324 simp(p-1)))), 325 gamsq(simp p)) 326 else 327 if z = '(1 . 1) and not numberp(prepsq b) then 328 multsq( 329 subsqnew( 330 derivativesq( 331 quotsq(gamsq(simp 'r),gamsq(addsq(simp 'r,subtrsq('(1 . 1),b)))), 332 'r,simp(p-1)), 333 a,'r), 334 quotsq( 335 multsq(multsq(simpx1(prepsq(multsq('(-1 . 1),a)),p,1), '(-1 . 1)), 336 gamsq(subtrsq('(1 . 1),b))), 337 gamsq(simp p))) 338 else 339 if z='(-1 . 1) and numberp prepsq(b) then 340 begin scalar c; integer k; 341 return multsq( 342 subsqnew( 343 derivativesq( addsq( 344 << 345 k:=prepsq(b) - 1; c:='(nil . 1); 346 while prepsq(k)>0 do 347 << 348 c:=addsq(c, multsq(gamsq(b), 349 simppochh(subtrsq(simp(1+k),simp 'r), 350 simp(prepsq(b)-1-k)))); 351 k:=k-1 352 >>; 353 c 354 >>, 355 quotsq( 356 multsq(gamsq(subtrsq(b,simp 'r)), 357 subtrsq(psisq(multsq(addsq(simp 'r,'(1 . 1)),'(1 . 2))), 358 psisq(multsq(simp 'r,'(1 . 2))))), 359 multsq((2 . 1), gamsq(subtrsq('(1 . 1),simp 'r))))), 360 'r,p-1), a, 'r), 361 quotsq( 362 multsq(simpx1(prepsq(multsq('(-1 . 1),a)),p,1), '(-1 . 1)), 363 multsq(gamsq(simp p),gamsq(simp b)))) 364 end 365 else 'fail$ 366 367%*********************************************************************** 368%* Particular cases * 369%*********************************************************************** 370 371 symbolic procedure ghf01(a,b,z); 372 if znak z then 373 multsq(gamsq(car b),multsq(bessmsq(subtrsq(car b,'(1 . 1)), 374 multsq('(2 . 1),simpx1(prepsq z,1,2))), 375 expdeg(z,multsq(subtrsq('(1 . 1),car b),'(1 . 2))))) 376 else 377 multsq(gamsq(car b),multsq(besssq(subtrsq(car b,'(1 . 1)), 378 multsq('(2 . 1),simpx1(prepsq(negsq z),1,2))),expdeg(negsq z, 379 multsq(subtrsq('(1 . 1),car b),'(1 . 2))))) $ 380 381 symbolic procedure ghf11(a,b,z); 382 if equal(car b,multsq('(2 . 1),car a)) then 383 multsq(multsq(gamsq(addsq('(1 . 2),car a)),expdeg(simp 'e, 384 multsq(z,'(1 . 2)))), 385 multsq(expdeg(multsq(z,'(1 . 4)),subtrsq('(1 . 2),car a)), 386 bessmsq(subtrsq(car a,'(1 . 2)),multsq(z,'(1 . 2))))) 387 else 388 if equal(car a,'(1 . 2)) and equal(car b,'(3 . 2)) then 389 multsq(multsq(simpx1('pi,1,2),'(1 . 2)), 390 if znak z then 391 quotsq(simpfunc('erfi,simpx1(prepsq z,1,2)),simpx1(prepsq z,1,2)) 392 else 393 quotsq(simpfunc('erf,simpx1(prepsq(negsq z),1,2)), 394 simpx1(prepsq(negsq z),1,2))) 395 else 396 if equal(car a,'(1 . 1)) and equal(car b,'(3 . 2)) and znak z then 397 multsq(multsq('(1 . 2),expdeg(simp 'e,z)), 398 multsq(simpfunc('erf,simpx1(prepsq z,1,2)),simpx1(prepsq 399 quotsq(simp('pi),z),1,2))) 400 else ghfexit(a,b,z)$ 401 402 symbolic procedure ghf21(a,b,z); 403 if and(equal(car a,'(1 . 2)),equal(cadr a,'(1 . 2)), 404 equal(car b,'(3 . 2)),znak(z)) 405 then 406 quotsq(simpfunc('asin,simpx1(prepsq(z),1,2)), 407 simpx1(prepsq(z),1,2)) 408 else 409 if ((equal(car a,'(1 . 2)) and equal(cadr a,'(1 . 1))) or 410 (equal(car a,'(1 . 1)) and equal(cadr a,'(1 . 2)))) and 411 equal(car b,'(3 . 2)) 412 then 413 << 414 if not znak(z) then 415 quotsq(simpfunc('atan,simpx1(prepsq(negsq z),1,2)), 416 simpx1(prepsq(negsq z),1,2)) else 417% if not equal(z,'(1 . 1)) then 418% quotsq(simpfunc('log,addsq('(1 . 1),simpx1(prepsq z,1,2))), 419% multsq(simpfunc('log,subtrsq('(1 . 1),simpx1(prepsq z,1,2))), 420% multsq('(2 . 1),simpx1(prepsq z,1,2)))) else 421 if not equal(z,'(1 . 1)) then 422 multsq(simpfunc('log,quotsq(addsq('(1 . 1),simpx1(prepsq z,1,2)), 423 subtrsq('(1 . 1),simpx1(prepsq z,1,2)))), 424 invsq(multsq('(2 . 1),simpx1(prepsq z,1,2)))) else 425 ghfexit(a,b,z) 426 >> 427 else 428 if and(equal(car a,'(1 . 1)),equal(cadr a,'(1 . 1)), 429 equal(car b,'(2 . 1)),not equal(z,'(1 . 1))) 430 then 431 quotsq(simpfunc('log,addsq('(1 . 1),negsq z)),negsq z) 432 else 433 if equal(subtrsq(addsq(car a,cadr a),car b),'(-1 . 2)) and 434 (equal(multsq('(2 . 1),car a),car b) or 435 equal(multsq('(2 . 1),cadr a),car b)) 436 then 437 multsq(expdeg(addsq('(1 . 1), 438 simpx1(prepsq(subtrsq('(1 . 1),z)),1,2)), 439 subtrsq('(1 . 1),car b)),expdeg('(2 . 1),addsq(car b,'(-1 . 1)))) 440 else 441 if z='(1 . 1) 442 and (not numberp prepsq subtrsq(car b,addsq(car a, cadr a)) 443 or prepsq(subtrsq(car b,addsq(car a, cadr a))) > 0 ) 444 then quotsq(multsq(gamsq(car b), 445 gamsq(subtrsq(car b,addsq(car a,cadr a))) ), 446 multsq(gamsq(subtrsq(car b,car a)), 447 gamsq(subtrsq(car b,cadr a)))) 448 else 449 if car a='(1 . 1) and cadr a='(1 . 1) and numberp prepsq car b and 450 prepsq car(b) > 0 and not(z='(1 . 1)) then 451 formula136(prepsq car b,z) else 452ghfexit(a,b,z)$ 453 454 symbolic procedure formula136(m,z); 455 begin scalar c; integer k; 456 c:='(nil . 1); k:=2; 457 while k<=m-1 do 458 << c:=addsq(c,quotsq(exptsq(subtrsq(z,'(1 . 1)),k), 459 multsq(exptsq(z,k),simp(m-k)))); 460 k:=k+1 461 >>; 462 c:=subtrsq(c,multsq(exptsq(quotsq(subtrsq(z,'(1 . 1)),z),m), 463 simpfunc('log,subtrsq('(1 . 1),z)))); 464 return 465 multsq(c, 466 quotsq(multsq(simp(m-1),z),exptsq(subtrsq(z,'(1 . 1)),2))); 467 end$ 468 469 symbolic procedure ghf12(a,b,z); 470 if equal(car a,'(3 . 4)) and (equal(car b,'(3 . 2)) and equal(cadr b, 471 '(7 . 4)) or equal(car b,'(7 . 4)) and equal(cadr b,'(3 . 2))) 472 and not znak z then 473 <<z:=multsq('(2 . 1),simpx1(prepsq(negsq z),1,2)); 474 multsq(quotsq(multsq('(3 . 1),simpx1('pi,1,2)), 475 multsq(simpx1(2,1,2),simpx1(prepsq z,3,2))), 476 simpfunc('intfs,z)) >> 477 else 478 if equal(car a,'(1 . 4)) and (equal(car b,'(1 . 2)) and equal(cadr b, 479 '(5 . 4)) or equal(car b,'(5 . 4)) and equal(cadr b,'(1 . 2))) 480 and not znak z then 481 <<z:=multsq((2 . 1),simpx1(prepsq(negsq z),1,2)); 482 multsq(quotsq(simpx1('pi,1,2),multsq(simpx1(2,1,2), 483 simpx1(prepsq z,1,2))),simpfunc('intfc,z)) >> 484 else ghfexit(a,b,z)$ 485 486symbolic inline procedure ghyper_fehlerf(); 487 rerror('specialf,139,"Wrong arguments to hypergeometric"); 488 489symbolic procedure hypergeom(u); 490 491begin scalar list1,list2,res,res1; 492 493if not (length(u) = 3) then ghyper_fehlerf(); 494if pairp u then list1 :=car u else ghyper_fehlerf(); 495if pairp cdr u then list2 := cadr u else ghyper_fehlerf(); 496if not pairp cddr u then ghyper_fehlerf(); 497 498if not eqcar(list1,'list) then ghyper_fehlerf(); 499if not eqcar(list2,'list) then ghyper_fehlerf(); 500 501list1 := for each x in cdr list1 collect simp reval x; 502list2 := for each x in cdr list2 collect simp reval x; 503res := ghfsq(list (length list1,length list2), 504 list1,list2,simp caddr u); 505res1 := prepsq res; 506return if eqcar(res1,'hypergeometric) then res else simp res1; 507 end; 508 509remflag('(hypergeometric),'full); 510put('hypergeometric,'simpfn,'hypergeom); 511 512% differentiation of hypergeometric function 513 514symbolic procedure dfform_hypergeometric(ghfform,dfvar,n); 515 begin scalar a,b,var,fct,result; 516 a:= cdr cadr ghfform; 517 b:= cdr caddr ghfform; 518 var := cadddr ghfform; 519 % diff. w.r.t. one of indizes --> return unchanged 520 if depends(a,dfvar) or depends(b,dfvar) 521 then result := !*kk2q {'df,ghfform,dfvar} 522 else << fct := simp!* {'quotient,retimes a,retimes b}; 523 result := simp!* {'hypergeometric, 524 'list . for each el in a collect {'plus, el, 1}, 525 'list . for each el in b collect {'plus, el, 1}, 526 var}; 527 result := multsq(fct,result); 528 if dfvar neq var then result := multsq(result,simp!*{'df,var,dfvar}) >>; 529 if n neq 1 530 then result := multsq(!*t2q((ghfform .** (n-1)) .* n),result); 531 return result; 532 end; 533 534put('hypergeometric,'dfform,'dfform_hypergeometric); 535 536 537% something is missing: 538 539algebraic let {hypergeometric({1/2,1/2},{3/2},-(~x)^2) => asinh(x)/x }; 540 541algebraic let hypergeometric({~a,~b},{~c},-(~z/(1-~z))) => 542 hypergeometric({a,c-b},{c},z) * (1-z)^a; 543 % Pfaff's reflection law 544 545flag ('(permutationof),'boolean); 546 547symbolic procedure permutationof(set1,set2); 548 length set1 = length set2 549 and not setdiff(set1,set2); 550 551algebraic let { 552 hypergeometric({},~lowerind,~z) => 553 3/(32*sqrt(2)*(-z)^(3/4))* 554 (cosh(2*(-z*4)^(1/4))*sin(2*(-z*4)^(1/4)) - 555 sinh(2*(-z*4)^(1/4))*cos(2*(-z*4)^(1/4))) 556 when permutationof(lowerind,{5/4,3/2,7/4}) 557 and numberp z and z < 0, 558 559 hypergeometric({},~lowerind,~z) => 560 1/(4*(-4*z)^(1/4))* 561 (sinh(2*(-z*4)^(1/4))*cos(2*(-z*4)^(1/4)) + 562 cosh(2*(-z*4)^(1/4))*sin(2*(-z*4)^(1/4))) 563 when permutationof(lowerind,{5/4,1/2,3/4}) 564 and numberp z and z < 0, 565 566 hypergeometric({},~lowerind,~z) => 567 1/(8*(-z)^(1/2))*sinh(2*(-z*4)^(1/4))*sin(2*(-z*4)^(1/4)) 568 when permutationof(lowerind,{3/4,5/4,3/2}) 569 and numberp z and z < 0, 570 571 hypergeometric({},~lowerind,~z) => 572 cosh(2*(-z*4)^(1/4))*cos(2*(-z*4)^(1/4)) 573 when permutationof(lowerind,{1/4,1/2,3/4}) 574 and numberp z and z < 0, 575 576 hypergeometric({},~lowerind,~z) => 577 3/(64*z^(3/4))*(sinh(4*z^(1/4)) -sin(4*z^(1/4))) 578 when permutationof(lowerind,{5/4,3/2,7/4}), 579 580 hypergeometric({},~lowerind,~z) => 581 1/(8*z^(1/4))*(sinh(4*z^(1/4)) +sin(4*z^(1/4))) 582 when permutationof(lowerind,{5/4,1/2,3/4}), 583 584 hypergeometric({},~lowerind,~z) => 585 1/(16*z^(1/2))*(cosh(4*z^(1/4)) -cos(4*z^(1/4))) 586 when permutationof(lowerind,{3/4,5/4,3/2}), 587 588 hypergeometric({},~lowerind,~z) => 589 1/2*(cosh(4*z^(1/4)) + cos(4*z^(1/4))) 590 when permutationof(lowerind,{1/4,1/2,3/4}) 591}; 592 593 594algebraic 595<< hypergeometric_rules:= 596 597{ hypergeometric({~a},{},~x) => (1-x)^(-a) when not(numberp x and x=1), 598 599% F(a;b;z) 600 601 hypergeometric({1/2},{5/2},~x) => 602 3/(4*x)*((1+2*x)/2*sqrt(pi/x)*erfi(sqrt(x))-e^x), 603 604 hypergeometric({1},{1/2},~x) => 605 1+sqrt(pi*x)*e^x*erf(sqrt(x)), 606 607 hypergeometric({1},{3/2},~x) => 608 1/2*sqrt(pi/x)*e^x*erf(sqrt(x)), 609 610 hypergeometric({1},{5/2},~x) => 611 3/(2*x)*(1/2*sqrt(pi/x)*e^(x)*erf(sqrt(x))-1), 612 613 hypergeometric({1},{7/2},~x) => 614 5/(4*x^2)*(3/2*sqrt(pi/x)*e^x*erf(sqrt(x))-3-2*x), 615 616 hypergeometric({3/2},{5/2},-~x) => 617 e^(-x)*hypergeometric({1},{5/2},x), 618 619 hypergeometric({3/2},{5/2},~x) => 620 3/(2*x)*(e^x-1/2*sqrt(pi/x)*erfi(sqrt(x))), 621 622 hypergeometric({5/2},{7/2},-~x) => 623 e^(-x)*hypergeometric({1},{7/2},x), 624 625 hypergeometric({7/2},{9/2},-~x) => 626 e^(-x)*hypergeometric({1},{9/2},x), 627 628 hypergeometric({~a},{~b},~x) => 629 a*(-x)^(-a)*m_gamma(a,-x) when b = a + 1, 630 631 hypergeometric({~a},{~b},~x) => 632 (a+1)*(e^(x)+(-x-a)*(-x)^(-a-1)*m_gamma(a+1,-x)) 633 when b = a + 2, 634 635% F(a,b;c;z) 636 637 hypergeometric({-1/2,1},{3/2},-~x) => 638 (1/2)*(1+(1+x)*(atan(sqrt(x))/sqrt(x))), 639 640 hypergeometric({-1/2,1},{3/2},~x) => 641 (1/2)*(1+(1-x)*(atanh(sqrt(x))/sqrt(x))), 642 643 hypergeometric({1/2,1},{5/2},-~x) => 644 (3/2*-x)*(1-(1+x)*(atan(sqrt(x))/sqrt(x))), 645 646 hypergeometric({1/2,1},{5/2},~x) => 647 (3/2*x)*(1-(1-x)*(atanh(sqrt(x))/sqrt(x))), 648 649 hypergeometric({~a + 1/2,~a},{1/2},~x) => 650 (1-x)^(-a)*cos(2*a*atan(sqrt(-x))), 651 652 hypergeometric({5/4,3/4},{1/2},~x) => 653 (1-x)^(-3/4)*cos(3/2*atan(sqrt(-x))), 654 655 hypergeometric({(~n+1)/2 + 1/2,(~n+1)/2},{1/2},~x) => 656 (1-x)^(-(n+1)/2)*cos((n+1)*atan(sqrt(-x))), 657 658 659 hypergeometric({7/4,5/4},{3/2},~x) => 660 2/3*(1-x)^(-3/4)/sqrt(-x)*sin(3/2*atan(sqrt(-x))), 661 662 hypergeometric({~a + 1/2,~a},{3/2},~x) => 663 ((1-x)^(1/2-a))/((2*a-1)*sqrt(-x))*sin((2*a-1)*atan(sqrt(-x))), 664 665 hypergeometric({(~n+2)/2 + 1/2,(~n+2)/2},{3/2},~x) => 666 ((1-x)^(1/2-(n+2)/2))/((2*(n+2)/2-1)*sqrt(-x))*sin((2*(n+2)/2-1) 667 *atan(sqrt(-x))), 668 669% F(a;b,c;z); 670 671 hypergeometric({-1/2},{1/2,1/2},-~x) => 672 cos(2*sqrt(x))+2*sqrt(x)*si(2*sqrt(x)), 673 674 hypergeometric({-1/2},{1/2,1/2},~x) => 675 cosh(2*sqrt(x))-2*x*shi(2*sqrt(x)), 676 677 hypergeometric({-1/2},{1/2,3/2},-~x) => 678 (1/2)*(cos(2*sqrt(x))+(sin(2*sqrt(x)))/(2*sqrt(x)) 679 +2*sqrt(x)*si(2*sqrt(x))), 680 681 hypergeometric({-1/2},{1/2,3/2},~x) => 682 (1/2)*(cosh(2*sqrt(x))+(sinh(2*sqrt(x)))/(2*sqrt(x)) 683 -2*sqrt(x)*shi(2*sqrt(x))), 684 685 hypergeometric({-1/2},{3/2,3/2},-~x) => 686 (1/4)*(cos(2*sqrt(x))+(sin(2*sqrt(x)))/(2*sqrt(x)) 687 +(1+2*x)*(si(2*sqrt(x)))/sqrt(x)), 688 689 hypergeometric({-1/2},{3/2,3/2},~x) => 690 (1/4)*(cosh(2*sqrt(x)) +(sinh(2*sqrt(x)))/(2*sqrt(x)) 691 +(1-2*x)*(shi(2*sqrt(x)))/sqrt(x)), 692 693 694 hypergeometric({1/2},{3/2,3/2},-~x) => 695 si(2*sqrt(x))/(2*sqrt(x)), 696 697 hypergeometric({1/2},{3/2,3/2},~x) => 698 shi(2*sqrt(x))/(2*sqrt(x)), 699 700 701 hypergeometric({1/2},{5/2,3/2},-~x) => 702 3/(8*-x)*(2*sqrt(x)*si(2*sqrt(x))-cos(2*sqrt(x))+ 703 (sin(2*sqrt(x)))/(2*sqrt(x))), 704 705 hypergeometric({1/2},{5/2,3/2},~x) => 706 3/(8*x)*(2*sqrt(x)*shi(2*sqrt(x))-cosh(2*sqrt(x))+ 707 (sinh(2*sqrt(x)))/(2*sqrt(x))), 708 709 710 hypergeometric({1},{3/4,5/4},~x) => 711 1/2*sqrt(pi/sqrt(-x))*(cos(2*sqrt(-x))*Fresnel_C(2*sqrt(-x)) 712 + sin(2*sqrt(-x))*Fresnel_S(2*sqrt(-x))), 713 714 hypergeometric({1},{5/4,7/4},~x) => 715 3*sqrt(pi)/(8*(sqrt(-x))^(3/2))*(sin(2*sqrt(-x)) 716 *Fresnel_C(2*sqrt(-x))-cos(2*sqrt(-x))*Fresnel_S(2*sqrt(-x))), 717 718 hypergeometric({5/2},{7/2,7/2},-~x) => 719 (75/(16*x^2))*(3*si(2*sqrt(x))/(2*sqrt(x)) 720 - 2*sin(2*sqrt(x))/sqrt(x) + cos(2*sqrt(x))), 721 722 hypergeometric({5/2},{7/2,7/2},~x) => 723 (75/(16*x^2))*(3*shi(2*sqrt(x))/(2*sqrt(x)) 724 - 2*sinh(2*sqrt(x))/sqrt(x) + cosh(2*sqrt(x))), 725 726 727 hypergeometric({~a},{~b,3/2},~x) => 728 -2^(1-2*a)*a*(sqrt(-x))^(-2*a)* 729 (gamma(2*a-1)*cos(a*pi)+Fresnel_S(2*sqrt(-x),2*a-1)) 730 731 when b = a + 1, 732 733 hypergeometric({~a},{~b,1/2},~x) => 734 2^(1-2*a)*a*(sqrt(-x))^(-2*a)* 735 (gamma(2*a)*cos(a*pi)-Fresnel_C(2*sqrt(-x),2*a)) 736 737 when b = a + 1 738}; 739 740let hypergeometric_rules; 741 742operator poisson!-charlier, toronto; 743 744let { toronto(~m,~n,~x) => 745 gamma(m/2 + 1/2)/factorial n * x^(2*n-2*m+1)*exp(-x^2) * 746 KummerM(m/2+1/2,1+n,x^2), 747 748 poisson!-charlier(~n,~nu,~x) => 749 Pochhammer(1 + nu-n,n)/(sqrt factorial n * x^(n/2))* 750 sum(Pochhammer(-n,i)*x^i/ 751 (Pochhammer(1+nu-n,i) * factorial i) 752 ,i,0,n) 753 }; 754>>; 755 756 757endmodule; 758 759end; 760 761 762 763 764