1% Tests of PM. 2 3% TESTS OF BASIC CONSTRUCTS. 4 5% These names are used with varing numbers of arguments here in this 6% test script, so to avoid messages of complaint I tag them as 7% variadic. 8lisp flag('(f h !@), 'variadic); 9 10 11 12operator f, h$ 13 14 15 16% A "literal" template. 17 18m(f(a),f(a)); 19 20 21t 22 23 24% Not literally equal. 25 26m(f(a),f(b)); 27 28 29 30%Nested operators. 31 32m(f(a,h(b)),f(a,h(b))); 33 34 35t 36 37 38% A "generic" template. 39 40m(f(a,b),f(a,?a)); 41 42 43{?a->b} 44 45m(f(a,b),f(?a,?b)); 46 47 48{?a->a,?b->b} 49 50 51% ??a takes "rest" of arguments. 52 53m(f(a,b),f(??a)); 54 55 56{??a->[a,b]} 57 58 59% But ?a does not. 60 61m(f(a,b),f(?a)); 62 63 64 65% Conditional matches. 66 67m(f(a,b),f(?a,?b _=(?a=?b))); 68 69 70 71m(f(a,a),f(?a,?b _=(?a=?b))); 72 73 74{?a->a,?b->a} 75 76 77% "plus" is symmetric. 78 79m(a+b+c,c+?a+?b); 80 81 82{?a->a,?b->b} 83 84 85%It is also associative. 86 87m(a+b+c,c+?a); 88 89 90{?a->a + b} 91 92 93% Note the effect of using multi-generic symbol is different. 94 95m(a+b+c,c+??c); 96 97 98{??c->[a,b]} 99 100 101%Flag h as associative. 102 103flag('(h),'assoc); 104 105 106 107m(h(a,b,d,e),h(?a,d,?b)); 108 109 110{?a->h(a,b),?b->e} 111 112 113% Substitution tests. 114 115s(f(a,b),f(a,?b)->?b^2); 116 117 118 2 119b 120 121 122s(a+b,a+b->a*b); 123 124 125a*b 126 127 128% "associativity" is used to group a+b+c in to (a+b) + c. 129 130s(a+b+c,a+b->a*b); 131 132 133a*b + c 134 135 136% Only substitute top at top level. 137 138s(a+b+f(a+b),a+b->a*b,inf,0); 139 140 141f(a + b) + a*b 142 143 144 145% SIMPLE OPERATOR DEFINITIONS. 146 147% Numerical factorial. 148 149operator nfac$ 150 151 152 153s(nfac(3),{nfac(0)->1,nfac(?x)->?x*nfac(?x-1)},1); 154 155 1563*nfac(2) 157 158 159s(nfac(3),{nfac(0)->1,nfac(?x)->?x*nfac(?x-1)},2); 160 161 1626*nfac(1) 163 164 165si(nfac(3),{nfac(0)->1,nfac(?x)->?x*nfac(?x-1)}); 166 167 1686 169 170 171% General factorial. 172 173operator gamma,fac; 174 175 176 177fac(?x _=Natp(?x)) ::- ?x*fac(?x-1); 178 179 180hold(?x*fac(?x - 1)) 181 182 183fac(0) :- 1; 184 185 1861 187 188 189fac(?x) :- Gamma(?x+1); 190 191 192gamma(?x + 1) 193 194 195fac(3); 196 197 1986 199 200 201fac(3/2); 202 203 204 3*sqrt(pi) 205------------ 206 4 207 208 209% Legendre polynomials in ?x of order ?n, ?n a natural number. 210 211operator legp; 212 213 214 215legp(?x,0) :- 1; 216 217 2181 219 220 221legp(?x,1) :- ?x; 222 223 224?x 225 226 227legp(?x,?n _=natp(?n)) 228 ::- ((2*?n-1)*?x*legp(?x,?n-1)-(?n-1)*legp(?x,?n-2))/?n; 229 230 231 (2*?n - 1)*?x*legp(?x,?n - 1) - (?n - 1)*legp(?x,?n - 2) 232hold(----------------------------------------------------------) 233 ?n 234 235 236legp(z,5); 237 238 239 4 2 240 z*(63*z - 70*z + 15) 241------------------------ 242 8 243 244 245legp(a+b,3); 246 247 248 3 2 2 3 249 5*a + 15*a *b + 15*a*b - 3*a + 5*b - 3*b 250--------------------------------------------- 251 2 252 253 254legp(x,y); 255 256 257legp(x,y) 258 259 260 261% TESTS OF EXTENSIONS TO BASIC PATTERN MATCHER. 262 263comment *: MSet[?exprn,?val] or ?exprn ::: ?val 264 assigns the value ?val to the projection ?exprn in such a way 265 as to store explicitly each form of ?exprn requested. *; 266 267 268Nosimp('mset,(t t)); 269 270 271 272Newtok '((!: !: !: !-) Mset); 273 274 275 276infix :::-; 277 278 279 280precedence Mset,RSetd; 281 282 283 284?exprn :::- ?val ::- (?exprn ::- (?exprn :- ?val )); 285 286 287hold(?exprn::-(?exprn:-?val)) 288 289 290scs := sin(?x)^2 + Cos(?x)^2 -> 1; 291 292 293 2 2 294scs := cos(?x) + sin(?x) ->1 295 296 297% The following pattern substitutes the rule sin^2 + cos^2 into a sum of 298% such terms. For 2n terms (ie n sin and n cos) the pattern has a worst 299% case complexity of O(n^3). 300 301operator trig,u; 302 303 304 305trig(?i) :::- Ap(+, Ar(?i,sin(u(?1))^2+Cos(u(?1))^2)); 306 307 308 2 2 309hold(trig(?i):-ap(plus,ar(?i,sin(u(?1)) + cos(u(?1)) ))) 310 311 312if si(trig 1,scs) = 1 then write("Pm ok") else Write("PM failed"); 313 314 315Pm ok 316 317 318if si(trig 10,scs) = 10 then write("Pm ok") else Write("PM failed"); 319 320 321Pm ok 322 323 324% The next one takes about 70 seconds on an HP 9000/350, calling UNIFY 325% 1927 times. 326 327% if si(trig 50,scs) = 50 then write("Pm ok") else Write("PM failed"); 328 329% Hypergeometric Function simplification. 330 331newtok '((!#) !#); 332 333 334*** # redefined 335 336 337flag('(#), 'symmetric); 338 339 340 341operator #,@,ghg; 342 343 344 345xx := ghg(4,3,@(a,b,c,d),@(d,1+a-b,1+a-c),1); 346 347 348xx := ghg(4,3,@(a,b,c,d),@(d,a - b + 1,a - c + 1),1) 349 350 351S(xx,sghg(3)); 352 353 354*** sghg declared operator 355 356ghg(4,3,@(a,b,c,d),@(d,a - b + 1,a - c + 1),1) 357 358 359s(ws,sghg(2)); 360 361 362ghg(4,3,@(a,b,c,d),@(d,a - b + 1,a - c + 1),1) 363 364 365yy := ghg(3,2,@(a-1,b,c/2),@((a+b)/2,c),1); 366 367 368 c a + b 369yy := ghg(3,2,@(a - 1,b,---),@(-------,c),1) 370 2 2 371 372 373S(yy,sghg(1)); 374 375 376 c a + b 377ghg(3,2,@(a - 1,b,---),@(-------,c),1) 378 2 2 379 380 381yy := ghg(3,2,@(a-1,b,c/2),@(a/2+b/2,c),1); 382 383 384 c a + b 385yy := ghg(3,2,@(a - 1,b,---),@(-------,c),1) 386 2 2 387 388 389S(yy,sghg(1)); 390 391 392 c a + b 393ghg(3,2,@(a - 1,b,---),@(-------,c),1) 394 2 2 395 396 397% Some Ghg theorems. 398 399flag('(@), 'symmetric); 400 401 402 403% Watson's Theorem. 404 405SGhg(1) := Ghg(3,2,@(?a,?b,?c),@(?d _=?d=(1+?a+?b)/2,?e _=?e=2*?c),1) -> 406 Gamma(1/2)*Gamma(?c+1/2)*Gamma((1+?a+?b)/2)*Gamma((1-?a-?b)/2+?c)/ 407 (Gamma((1+?a)/2)*Gamma((1+?b)/2)*Gamma((1-?a)/2+?c) 408 *Gamma((1-?b)/2+?c)); 409 410 411 1 + ?a + ?b 412sghg(1) := ghg(3,2,@(?a,?b,?c),@(?d _= ?d=-------------,?e _= ?e=2*?c),1)->( 413 2 414 415 - ?a - ?b + 2*?c + 1 2*?c + 1 416 sqrt(pi)*gamma(-----------------------)*gamma(----------) 417 2 2 418 419 ?a + ?b + 1 - ?a + 2*?c + 1 420 *gamma(-------------))/(gamma(------------------) 421 2 2 422 423 - ?b + 2*?c + 1 ?a + 1 ?b + 1 424 *gamma(------------------)*gamma(--------)*gamma(--------)) 425 2 2 2 426 427 428% Dixon's theorem. 429 430SGhg(2) := Ghg(3,2,@(?a,?b,?c),@(?d _=?d=1+?a-?b,?e _=?e=1+?a-?c),1) -> 431 Gamma(1+?a/2)*Gamma(1+?a-?b)*Gamma(1+?a-?c)*Gamma(1+?a/2-?b-?c)/ 432 (Gamma(1+?a)*Gamma(1+?a/2-?b)*Gamma(1+?a/2-?c)*Gamma(1+?a-?b-?c)); 433 434 435sghg(2) := ghg(3,2,@(?a,?b,?c),@(?d _= ?d=1 + ?a - ?b,?e _= ?e=1 + ?a - ?c),1)-> 436 437 ?a - 2*?b - 2*?c + 2 438 (gamma(?a - ?b + 1)*gamma(?a - ?c + 1)*gamma(----------------------) 439 2 440 441 ?a + 2 442 *gamma(--------))/(gamma(?a - ?b - ?c + 1)*gamma(?a + 1) 443 2 444 445 ?a - 2*?b + 2 ?a - 2*?c + 2 446 *gamma(---------------)*gamma(---------------)) 447 2 2 448 449 450SGhg(3) := Ghg(?p,?q,@(?a,??b),@(?a,??c),?z) 451 -> Ghg(?p-1,?q-1,@(??b),@(??c),?z); 452 453 454sghg(3) := 455 456ghg(?p,?q,@(??b,?a),@(??c,?a),?z)->ghg(?p - 1,?q - 1,@(??b),@(??c),?z) 457 458 459SGhg(9) := Ghg(1,0,@(?a),?b,?z ) -> (1-?z)^(-?a); 460 461 462 1 463sghg(9) := ghg(1,0,@(?a),?b,?z)->--------------- 464 ?a 465 ( - ?z + 1) 466 467SGhg(10) := Ghg(0,0,?a,?b,?z) -> E^?z; 468 469 470 ?z 471sghg(10) := ghg(0,0,?a,?b,?z)->e 472 473SGhg(11) := Ghg(?p,?q,@(??t),@(??b),0) -> 1; 474 475 476sghg(11) := ghg(?p,?q,@(??t),@(??b),0)->1 477 478 479% If one of the bottom parameters is zero or a negative integer the 480% hypergeometric functions may be singular, so the presence of a 481% functions of this type causes a warning message to be printed. 482 483% Note it seems to have an off by one level spec., so this may need 484% changing in future. 485% 486% Reference: AS 15.1; Slater, Generalized Hypergeometric Functions, 487% Cambridge University Press,1966. 488 489s(Ghg(3,2,@(a,b,c),@(b,c),z),SGhg(3)); 490 491 492ghg(2,1,@(a,b),@(b),z) 493 494 495si(Ghg(3,2,@(a,b,c),@(b,c),z),{SGhg(3),Sghg(9)}); 496 497 498 1 499------------- 500 a 501 ( - z + 1) 502 503 504S(Ghg(3,2,@(a-1,b,c),@(a-b,a-c),1),sghg 2); 505 506 507 a - 2*b - 2*c + 1 a + 1 508 gamma(a - b)*gamma(a - c)*gamma(-------------------)*gamma(-------) 509 2 2 510--------------------------------------------------------------------- 511 a - 2*b + 1 a - 2*c + 1 512 gamma(a - b - c)*gamma(-------------)*gamma(-------------)*gamma(a) 513 2 2 514 515 516end; 517 518Tested on x86_64-pc-windows CSL 519Time (counter 1): 0 ms 520 521End of Lisp run after 0.00+0.04 seconds 522real 0.20 523user 0.03 524sys 0.03 525