1float2bf : true$ 2true; 3 4/* Rodrigue's formulae */ 5 6(jacobi_p_rod(n,a,b,x) := block([ an, rho, g ], 7 an : (-1)^n * 2^n * n!, 8 rho : (1-t)^a * (1 + t)^b, 9 g : 1-t^2, 10 rat(subst(x,t, diff(rho * g^n,t,n) / (an * rho)))), 11 12 gen_laguerre_rod(n, a, x) := block([an, rho, g], 13 an : n!, 14 rho : exp(-t) * t^a, 15 g : t, 16 rat(subst(x,t, diff(rho * g^n,t,n) / (an * rho)))), 17 18 hermite_rod(n,x) := block([an, rho, g], 19 an : (-1)^n, 20 rho : exp(-t^2), 21 g : 1, 22 rat(subst(x,t, diff(rho * g^n,t,n) / (an * rho)))), 23 24/*See A&S 10.1.25 page 439. */ 25 26 spherical_bessel_j_rod(n,x) := block([sofar,k], 27 sofar : sin(x)/x, 28 for k : 1 thru n do ( 29 sofar : -diff(sofar,x) / x 30 ), 31 x^n * sofar), 32 33 spherical_bessel_y_rod(n,x) := block([sofar,k], 34 sofar : cos(x) / x, 35 for k : 1 thru n do ( 36 sofar : -diff(sofar,x) / x 37 ), 38 -x^n * sofar), 39 40 all_functions : [jacobi_p, 41 ultraspherical, 42 assoc_legendre_p, 43 legendre_q, 44 assoc_legendre_q, 45 chebyshev_t, 46 chebyshev_u, 47 laguerre, 48 gen_laguerre, 49 legendre_p, 50 hermite, 51 spherical_hankel2, 52 spherical_hankel1, 53 spherical_bessel_j, 54 spherical_bessel_y, 55 assoc_leg_cos, 56 spherical_harmonic], 57 58 59 errors_found : [ ], 60 tests_pass : [ ], 61 62 zerop(e) := is(0=e), 63 64 check_zero_list(e) := block([ k, okay,n], 65 kill(labels), 66 okay : true, 67 k : 0, 68 n : length(e), 69 while okay and k < n do ( 70 k : k + 1, 71 if not zerop(e[ k ]) then ( 72 okay : false 73 ) 74 ), 75 if okay then ( 76 tests_pass : endcons(test_name, tests_pass), 77 print("okay: ", test_name) 78 ) else ( 79 print("error: ", test_name), 80 print("should vanish = ", e[ k ]), 81 errors_found : endcons(test_name, errors_found) 82 ) 83 ), 84 85 check_true_list(e) := block([ ], 86 if member(false, e) then ( 87 print("error: ", test_name), 88 errors_found : endcons(test_name, errors_found) 89 ) else ( 90 print("okay: ", test_name), 91 tests_pass : endcons(test_name, tests_pass) 92 ) 93 ), 94 950); 960; 97 98/* Jacobi Rodrigues test*/ 99 100(print (test_name : "Jacobi Rodrigues test"), 101 makelist(jacobi_p(k,a,b,x) - jacobi_p_rod(k,a,b,x), k,0,5), 102 rat(%%)); 103''(makelist (0, k, 0, 5)); 104 105/* gen_laguerre Rodrigues*/ 106 107(print (test_name : "gen_laguerre Rodrigues"), 108 makelist(gen_laguerre(k,a,x) - gen_laguerre_rod(k,a,x), k,0,5), 109 rat(%%)); 110''(makelist (0, k, 0, 5)); 111 112/* hermite Rodrigues*/ 113 114(print (test_name : "hermite Rodrigues"), 115 makelist(hermite(k,x) - hermite_rod(k,x), k,0,5), 116 rat(%%)); 117''(makelist (0, k, 0, 5)); 118 119/* spherical_bessel_j Rodrigues*/ 120 121(print (test_name : "spherical_bessel_j Rodrigues"), 122 makelist(spherical_bessel_j(k,x) - spherical_bessel_j_rod(k,x), k,0,5), 123 ratsimp(exponentialize(%%))); 124''(makelist (0, k, 0, 5)); 125 126/* spherical_bessel_y Rodrigues*/ 127 128(print (test_name : "spherical_bessel_y Rodrigues"), 129 makelist(spherical_bessel_y(k,x) - spherical_bessel_y_rod(k,x), k,0,5), 130 ratsimp(exponentialize(%%))); 131''(makelist (0, k, 0, 5)); 132 133/*---------------*/ 134 135(print (test_name : "spherical harmonic orthogonality"), 136 137 f(n1,m1,n2,m2) := defint(defint(trigreduce(spherical_harmonic(n1,m1,th,ph) * 138 conjugate(spherical_harmonic(n2,m2,th,ph)) * sin(th)),th,0,%pi),ph,0,2*%pi), 139 0); 1400; 141 142block ([sofar : [ ]], 143 for n1 : 0 thru 2 do ( 144 for m1 : -n1 thru n1 do ( 145 for n2 : 0 thru 2 do ( 146 for m2 : - n2 thru n2 do ( 147 sofar : cons(f(n1,m1,n2,m2) - kron_delta(n1,n2) * kron_delta(m1,m2), sofar))))), 148 if every (zerop, sofar) then true else sofar); 149true; 150 151/*----------------*/ 152 153/* See A&S 22.3.14 page 776 and 22.5.4 page 777. 154*/ 155 156(print (test_name : "A&S 22.3.14"), 157 makelist(ultraspherical(n,a,cos(x))/a - 2*cos(n*x)/n,n,1, 3), 158 limit(%%,a,0), 159 ratsimp(trigreduce(%%))), 160''(makelist (0, 3)); 161 162/* See A&S 22.3.14 page 776 163*/ 164 165(print (test_name : "A&S 22.3.14"), 166 makelist(chebyshev_t(n,cos(x)) - cos(n*x),n,0, 3), 167 ratsimp(trigreduce(%%)))$ 168''(makelist (0, 4)); 169 170/* See A&S 22.3.15 page 776 171*/ 172 173(print (test_name : "A&S 22.3.15"), 174 makelist(sin(x) * chebyshev_u(n,cos(x)) - sin((n+1)*x),n,0, 4), 175 ratsimp(trigreduce(%%)))$ 176''(makelist (0, 5)); 177 178/* See A&S 22.7.15 page 782. 179*/ 180 181(print (test_name : "A&S 22.7.15"), 182 jacobi_rec(n) := (n + a/2+b/2+1)*(1-x)*jacobi_p(n,a+1,b,x) 183 - (n+a+1)*jacobi_p(n,a,b,x) + (n+1)*jacobi_p(n+1,a,b,x), 184 makelist(ratsimp(jacobi_rec(n)),n,0, 7)); 185''(makelist (0, 8)); 186 187/* See A&S 22.7.16 page 782 188*/ 189 190(print (test_name : "A&S 22.7.16"), 191 jacobi_rec(n) := (n + a/2+b/2+1)*(1+x)*jacobi_p(n,a,b+1,x) 192 - (n+b+1)*jacobi_p(n,a,b,x) - (n+1)*jacobi_p(n+1,a,b,x), 193 makelist(ratsimp(jacobi_rec(n)),n,0, 7)); 194''(makelist (0, 8)); 195 196/* See A&S 22.7.17 page 782 197*/ 198 199(print (test_name : "A&S 22.7.17"), 200 jacobi_rec(n) := (1-x)*jacobi_p(n,a+1,b,x) 201 + (1+x)*jacobi_p(n,a,b+1,x) - 2*jacobi_p(n,a,b,x), 202 makelist(ratsimp(jacobi_rec(n)),n,0, 7)); 203''(makelist (0, 8)); 204 205/* See A&S 22.7.18 page 782 206*/ 207 208test_name : "A&S 22.7.18"$ 209 210jacobi_rec(n) := (2*n+a+b)*jacobi_p(n,a-1,b,x) 211 - (n+a+b)*jacobi_p(n,a,b,x) + (n+b)*jacobi_p(n-1,a,b,x)$ 212 213check_zero_list(makelist(ratsimp(jacobi_rec(n)),n,1, 7))$ 214 215/* See A&S 22.7.19 page 782 216*/ 217 218test_name : "A&S 22.7.19"$ 219 220jacobi_rec(n) := (2*n+a+b)*jacobi_p(n,a,b-1,x) 221 - (n+a+b)*jacobi_p(n,a,b,x) - (n+a)*jacobi_p(n-1,a,b,x)$ 222 223check_zero_list(makelist(ratsimp(jacobi_rec(n)),n, 1, 7))$ 224 225/* See A&S 22.7.20 page 782 226*/ 227 228test_name : "A&S 22.7.20"$ 229 230jacobi_rec(n) := jacobi_p(n,a,b-1,x) 231 - jacobi_p(n,a-1,b,x) - jacobi_p(n-1,a,b,x)$ 232 233check_zero_list(makelist(ratsimp(jacobi_rec(n)),n,1, 2))$ 234 235/* See A&S 22.7.21 page 782 236*/ 237 238test_name : "A&S 22.7.21"$ 239 240ultraspherical_rec(n) := 2*a*(1-x^2)*ultraspherical(n-1,a+1,x) 241- (2*a+n-1) * ultraspherical(n-1,a,x) + n*x*ultraspherical(n,a,x)$ 242 243foo : makelist(ratsimp(ultraspherical_rec(n)),n,1, 7)$ 244foo : append(ev(foo,a=4/5), ev(foo, a= -2/3), ev(foo, a=8))$ 245check_zero_list(foo)$ 246 247 248/* See A&S 22.7.23 page 782; funny 22.7.22 is missing a lhs. 249 Maxima lacks simplification rules to simplify linear 250 combinations of gamma functions. 251*/ 252 253test_name : "A&S 22.7.23"$ 254 255ultraspherical_rec(n) := (n+a)*ultraspherical(n+1,a-1,x) 256- (a-1) * (ultraspherical(n+1,a,x) - ultraspherical(n-1,a,x))$ 257 258foo : makelist(ratsimp(ultraspherical_rec(n)),n,1, 7)$ 259foo : append(ev(foo,a=4/5,rat), ev(foo, a=8,rat))$ 260check_zero_list(foo)$ 261 262 263/* See A&S 22.7.24 page 782 264*/ 265 266test_name : "A&S 22.7.24"$ 267 268cheb_rec(n,m) := 2 * chebyshev_t(m,x) * chebyshev_t(n,x) 269- chebyshev_t(n+m,x) - chebyshev_t(n-m,x)$ 270 271foo : makelist(makelist(cheb_rec(n, m), m,1,n),n,1, 9)$ 272foo : rat(foo)$ 273foo : apply(append,foo)$ 274check_zero_list(foo)$ 275 276 277/* See A&S 22.7.25 page 782 278*/ 279 280test_name : "A&S 22.7.25"$ 281 282cheb_rec(n,m) := 2*(x^2-1)*chebyshev_u(m-1,x) * chebyshev_u(n-1,x) 283 - chebyshev_t(n+m,x) + chebyshev_t(n-m,x)$ 284 285foo : makelist(makelist(cheb_rec(n, m), m,1,n),n,1, 9)$ 286foo : rat(foo)$ 287foo : apply(append,foo)$ 288check_zero_list(foo)$ 289 290 291/* See A&S 22.7.26 page 782 292*/ 293 294test_name : "A&S 22.7.26"$ 295 296cheb_rec(n,m) := 2*chebyshev_t(m,x) * chebyshev_u(n-1,x) 297 -chebyshev_u(n+m-1,x) - chebyshev_u(n-m-1,x)$ 298 299foo : makelist(makelist(cheb_rec(n, m), m, 0, n-1),n,1, 9)$ 300foo : rat(foo)$ 301foo : apply(append,foo)$ 302check_zero_list(foo)$ 303 304/* See A&S 22.7.27 page 782 305*/ 306 307test_name : "A&S 22.7.27"$ 308 309cheb_rec(n,m) := 2*chebyshev_t(n,x) * chebyshev_u(m-1,x) 310 -chebyshev_u(n+m-1,x) + chebyshev_u(n-m-1,x)$ 311 312foo : makelist(makelist(cheb_rec(n, m), m,1,n-1),n, 2, 10)$ 313foo : rat(foo)$ 314foo : apply(append,foo)$ 315check_zero_list(foo)$ 316 317/* See A&S 22.7.28 page 782 318*/ 319 320test_name : "A&S 22.7.28"$ 321 322cheb_rec(n) := 2*chebyshev_t(n,x) * chebyshev_u(n-1,x) 323 -chebyshev_u(2*n-1,x)$ 324 325foo : makelist(cheb_rec(n), n,2,10)$ 326foo : rat(foo)$ 327check_zero_list(foo)$ 328 329/* See A&S 22.7.29 page 783 330*/ 331 332test_name : "A&S 22.7.29"$ 333 334lag_rec(n) := gen_laguerre(n,1/3+1,x) -((x-n) * gen_laguerre(n,1/3,x) + (1/3+n) * gen_laguerre(n-1,1/3,x))/x$ 335check_zero_list(ratsimp(makelist(lag_rec(n),n,1, 10)))$ 336 337lag_rec(n) := gen_laguerre(n,5+1,x) -((x-n) * gen_laguerre(n,5,x) + (5+n) * gen_laguerre(n-1,5,x))/x$ 338check_zero_list(ratsimp(makelist(lag_rec(n),n,1,10)))$ 339 340lag_rec(n) := gen_laguerre(n,1,x) -((x-n) * gen_laguerre(n,0,x) + (n) * gen_laguerre(n-1,0,x))/x$ 341check_zero_list(ratsimp(makelist(lag_rec(n),n,1,10)))$ 342 343 344/* See A&S 22.7.30 page 783 345*/ 346 347test_name : "A&S 22.7.30"$ 348 349declare(a,integer)$ 350 351lag_rec(n) := gen_laguerre(n,a-1,x) - gen_laguerre(n,a,x) + gen_laguerre(n-1,a,x)$ 352 353check_zero_list(ratsimp(makelist(lag_rec(n),n,1, 10)))$ 354remove(a,integer)$ 355 356 /* See A&S 22.7.31 page 783 357 */ 358 359test_name : "A&S 22.7.31"$ 360 361lag_rec(n) := gen_laguerre(n,a+1,x) - (n+a+1)*gen_laguerre(n,a,x) / x + (n+1)*gen_laguerre(n+1,a,x)/x$ 362check_zero_list(ratsimp(makelist(lag_rec(n),n,0, 7)))$ 363 364 /* See A&S 22.7.32 page 783 365 */ 366 367test_name : "A&S 22.7.32"$ 368 369lag_rec(n) := gen_laguerre(n,a-1,x) - (n+1)*gen_laguerre(n+1,a,x) / (n+a) + (n+1-x)*gen_laguerre(n,a,x)/(n+a)$ 370check_zero_list(ratsimp(makelist(lag_rec(n),n,0,7)))$ 371 372/* See A&S 22.2.1 page 774 373*/ 374 375test_name : "A&S 22.2.1"$ 376 377f(a,b,n,m) := integrate((1-x)^a*(1+x)^b * jacobi_p(n,a,b,x) 378 * jacobi_p(m,a,b,x),x,-1,1)$ 379foo : makelist(makelist(f(1/2,-1/2,n,m), n, 0, m - 1), m, 1, 5)$ 380foo : apply(append,foo)$ 381check_zero_list(foo)$ 382 383foo : makelist(makelist(f(1/3,2/3,n,m), n, 0, m - 1), m, 1, 5)$ 384foo : apply(append,foo)$ 385check_zero_list(foo)$ 386 387 388test_name : "jacobi normalization"$ 389 390stand(n) := expand(jacobi_p(n,a,b,1) - binomial(n+a,n))$ 391check_zero_list(makelist(stand(n),n,0, 7))$ 392 393test_name : "A&S 22.2.3"$ 394 395assume(a > 1/2)$ 396baz(n,m) := 'integrate((1-x^2)^(a-1/2) * ultraspherical(n,a,x)* ultraspherical(m,a,x),x,-1,1)$ 397foo : makelist(makelist(ev(baz(n,m),integrate),n,0,m-1),m,1, 2)$ 398foo : ev(foo, rat)$ 399forget(a > 1/2)$ 400foo : apply(append,foo)$ 401check_zero_list(foo)$ 402 403/* See A&S 22.2.3 page 774; Maxima doesn't know enough about 404the gamma functions to simplify these expressions to zero. 405*/ 406 407test_name : "A&S 22.2.3"$ 408 409stand(n) := ultraspherical(n,a,1) - binomial(n+2*a-1,n)$ 410foo : makelist(stand(n),n,0, 7)$ 411foo : append(ev(foo,a=2/3), ev(foo, a=7), ev(foo,a=1/3))$ 412foo : rat(foo)$ 413check_zero_list(foo)$ 414 415/* See A&S 22.2.10 page 774 416*/ 417 418test_name : "A&S 22.2.10"$ 419 420f(n,m) := integrate(legendre_p(n,(rat(x))) * legendre_p(m,rat(x)),x,-1,1) - kron_delta(n,m) * 2 /(2*n+1)$ 421foo : makelist(makelist(f(n,m),n,0,3),m,0,3)$ 422foo : apply(append, foo)$ 423check_zero_list(foo)$ 424 425test_name : "legendre poly normalization"$ 426stand(n) := legendre_p(n,1) -1$ 427foo : makelist(stand(n),n,0,35)$ 428check_zero_list(foo)$ 429 430/* See A&S 22.2.12 page 774 431*/ 432 433test_name : "A&S 22.2.12"$ 434baz(n,m) := 'integrate(gen_laguerre(n,a,x) * gen_laguerre(m,a,x) * exp(-x)*x^a,x,0,inf) - kron_delta(n,m) * gamma(a+n+1)/n!$ 435foo : makelist(makelist(baz(n,m),n,0,5),m,0,5)$ 436foo : apply(append, ev(foo,a=2/3, integrate))$ 437check_zero_list(foo)$ 438 439 440/* See A&S 22.2.13 page 775 441*/ 442 443test_name : "A&S 22.2.13"$ 444baz(n,m) := 'integrate(laguerre(n,x) * laguerre(m,x) * exp(-x),x,0,inf) - kron_delta(n,m)$ 445foo : makelist(makelist(ev(baz(n,m), integrate), n,0,4),m,0,4)$ 446foo : rat(foo)$ 447foo : apply(append, foo)$ 448check_zero_list(foo)$ 449 450/* See A&S 22.2.14 page 775 451*/ 452 453test_name : "A&S 22.2.14"$ 454baz(n,m) := 'integrate(hermite(n,x) * hermite(m,x) * exp(-x^2),x,-inf,inf) - kron_delta(n,m) * sqrt(%pi) * 2^n * n!$ 455foo : makelist(makelist(ev(baz(n,m), integrate), n,0,4),m,0,4)$ 456foo : apply(append, foo)$ 457check_zero_list(foo)$ 458 459 460/* Some Christoffel-Darboux sum formulae 461*/ 462 463/* See A&S 22.12.2 page 785 464*/ 465 466test_name : "A&S 22.12.2"$ 467baz(n) := sum(chebyshev_t(2*k,x) ,k,0,n) - (1 + chebyshev_u(2*n,x))/2$ 468foo : makelist(rat(baz(n)),n,0,7)$ 469check_zero_list(foo)$ 470 471/* See A&S 22.12.3 page 785 472*/ 473 474test_name : "A&S 22.12.3"$ 475baz(n) := sum(chebyshev_t(2*k+1,x) ,k,0,n-1) - chebyshev_u(2*n-1,x)/2$ 476foo : makelist(rat(baz(n)),n,1, 7)$ 477check_zero_list(foo)$ 478 479/* See A&S 22.12.4 page 785 480*/ 481 482test_name : "A&S 22.12.4"$ 483baz(n) := sum(chebyshev_u(2*k,x) ,k,0,n) - (1-chebyshev_t(2*n+2,x))/(2 * (1-x^2))$ 484foo : makelist(rat(baz(n)),n,1, 11)$ 485check_zero_list(foo)$ 486 487/* See A&S 22.2.12.5 page 785 488*/ 489 490test_name : "A&S 22.12.5"$ 491baz(n) := sum(chebyshev_u(2*k+1,x) ,k,0,n-1) - (x-chebyshev_t(2*n+1,x))/(2 * (1-x^2))$ 492foo : makelist(rat(baz(n)),n,1,7)$ 493check_zero_list(foo)$ 494 495/* See A&S 22.12.6 page 785 496*/ 497 498test_name : "A&S 22.12.6"$ 499baz(n) := gen_laguerre(n,a+b+1,x+y) - sum(gen_laguerre(k,a,x) * gen_laguerre(n-k,b,y),k,0,n)$ 500foo : makelist(rat(baz(n)),n,1, 7)$ 501check_zero_list(foo)$ 502 503/* See A&S 22.12.7 page 785 504*/ 505 506test_name : "A&S 22.12.7"$ 507baz(n) := gen_laguerre(n,a,x*y) - 508sum(binomial(n+a,m) * y^(n-m) * (1-y)^m * gen_laguerre(n-m,a,x),m,0,n)$ 509foo : makelist(rat(baz(n)),n,1,3)$ 510check_zero_list(foo)$ 511 512/* See A&S 22.12.8 page 785 513*/ 514 515test_name : "A&S 22.12.8"$ 516baz(n) := hermite(n,x+y) - sum(binomial(n, k) * hermite(k,sqrt(2)*x) * hermite(n-k,sqrt(2)*y),k,0,n) / 2^(n/2)$ 517foo : makelist(rat(baz(n)),n,0, 7)$ 518check_zero_list(foo)$ 519 520/* See A&S 22.5.17 page 778 521*/ 522 523test_name : "A&S 22.5.17"$ 524baz(n,m) := gen_laguerre(n,m,x) - (-1)^m * diff(laguerre(n+m,x),x,m)$ 525foo : makelist(makelist(baz(i,j),j,0,i),i,0,7)$ 526foo : rat(foo)$ 527foo : apply(append,foo)$ 528check_zero_list(foo)$ 529 530/* See A&S 22.7.29 page 783 531*/ 532 533test_name : "A&S 22.7.29"$ 534baz(n) := gen_laguerre(n,a+1,x) - ((x-n)*gen_laguerre(n,a,x) + (a+n)*gen_laguerre(n-1,a,x))/x$ 535foo : makelist(baz(i),i,1,8)$ 536foo : rat(foo)$ 537check_zero_list(foo)$ 538 539/* float tests 540*/ 541 542test_name : "A&S 22.12.2 - float"$ 543orthopoly_returns_intervals : true$ 544listarith : true$ 545 546xargs(e) := if intervalp(e) then args(e) else [e,0]$ 547 548f(n,x) := block([p,q], 549 p : sum(xargs(chebyshev_t(2*i,x)),i,0,n), 550 q : xargs(chebyshev_u(2*n,x)) / 2, 551 [part(p,1) - part(q,1), part(p,2) + part(q,2)])$ 552 553gomer : []$ 554for i : -50 thru 50 do ( 555 foo : f(10, 1.25 * i), 556 e : part(foo,2), 557 foo : part(foo,1), 558 if (abs(foo - 0.5) <= e) then foo : 0 else foo : 1, 559 gomer : cons(foo, gomer))$ 560check_zero_list(gomer)$ 561 562gomer : []$ 563for i : -100 thru 100 do ( 564 foo : f(35, 0.01 * i), 565 e : part(foo,2), 566 foo : part(foo,1), 567 if (abs(foo - 0.5) <= e) then foo : 0 else foo : 1, 568 gomer : cons(foo, gomer))$ 569check_zero_list(gomer)$ 570 571gomer : []$ 572for i : 1 thru 100 do ( 573 foo : f(35, -0.01 * i), 574 e : part(foo,2), 575 foo : part(foo,1), 576 if (abs(foo - 0.5) <= e) then foo : 0 else foo : 1, 577 gomer : cons(foo, gomer))$ 578check_zero_list(gomer)$ 579 580gomer : []$ 581for i : -100 thru 100 do ( 582 foo : f(35, 0.01 * i + %i * 0.2), 583 e : part(foo,2), 584 foo : part(foo,1), 585 if (abs(foo - 0.5) <= e) then foo : 0 else foo : 1, 586 gomer : cons(foo, gomer))$ 587check_zero_list(gomer)$ 588 589gomer : []$ 590for i : -100 thru 100 do ( 591 foo : f(36, -0.01 * i + %i * 0.2), 592 e : part(foo,2), 593 foo : part(foo,1), 594 if (abs(foo - 0.5) <= e) then foo : 0 else foo : 1, 595 gomer : cons(foo, gomer))$ 596check_zero_list(gomer)$ 597 598test_name : "A&S 22.12.3 - float"$ 599 600f(n,x) := block([p,q], 601 p : sum(xargs(chebyshev_t(2*i+1,x)),i,0,n-1), 602 q : xargs(chebyshev_u(2*n-1,x)) / 2, 603 [part(p,1) - part(q,1), part(p,2) + part(q,2)])$ 604 605 606gomer : []$ 607for i : -50 thru 50 do ( 608 foo : f(10, 1.25 * i), 609 e : part(foo,2), 610 foo : part(foo,1), 611 if (abs(foo) <= e) then foo : 0 else foo : 1, 612 gomer : cons(foo, gomer))$ 613check_zero_list(gomer)$ 614 615gomer : []$ 616for i : -100 thru 100 do ( 617 foo : f(35, 0.01 * i), 618 e : part(foo,2), 619 foo : part(foo,1), 620 if (abs(foo) <= e) then foo : 0 else foo : 1, 621 gomer : cons(foo, gomer))$ 622check_zero_list(gomer)$ 623 624gomer : []$ 625for i : 1 thru 100 do ( 626 foo : f(35, -0.01 * i), 627 e : part(foo,2), 628 foo : part(foo,1), 629 if (abs(foo) <= e) then foo : 0 else foo : 1, 630 gomer : cons(foo, gomer))$ 631check_zero_list(gomer)$ 632 633gomer : []$ 634for i : -100 thru 100 do ( 635 foo : f(35, 0.01 * i + %i * 0.2), 636 e : part(foo,2), 637 foo : part(foo,1), 638 if (abs(foo) <= e) then foo : 0 else foo : 1, 639 gomer : cons(foo, gomer))$ 640check_zero_list(gomer)$ 641 642gomer : []$ 643for i : -100 thru 100 do ( 644 foo : f(36, -0.01 * i + %i * 0.2), 645 e : part(foo,2), 646 foo : part(foo,1), 647 if (abs(foo) <= e) then foo : 0 else foo : 1, 648 gomer : cons(foo, gomer))$ 649check_zero_list(gomer)$ 650 651 652test_name : "G&R 8.974-3-float"$ 653 654f(n,a,x) := block([p,q], 655 p : sum(xargs(gen_laguerre(i,a,x)),i,0,n), 656 q : xargs(gen_laguerre(n,a+1,x)), 657 [part(p,1) - part(q,1), part(p,2) + part(q,2)])$ 658 659gomer : []$ 660for i : -100 thru 100 do ( 661 foo : f(11,0.4, float(0.01 * i)), 662 e : part(foo,2), 663 foo : part(foo,1), 664 if (abs(foo) <= e) then foo : 0 else foo : 1, 665 gomer : cons(foo,gomer))$ 666check_zero_list(gomer)$ 667 668gomer : []$ 669for i : -100 thru 100 do ( 670 foo : f(12,0.4, float(0.01 * i)), 671 e : part(foo,2), 672 foo : part(foo,1), 673 if (abs(foo) <= e) then foo : 0 else foo : 1, 674 gomer : cons(foo,gomer))$ 675check_zero_list(gomer)$ 676 677gomer : []$ 678for i : -100 thru 100 do ( 679 foo : f(12,-0.4, float(0.01 * i)), 680 e : part(foo,2), 681 foo : part(foo,1), 682 if (abs(foo) <= e) then foo : 0 else foo : 1, 683 gomer : cons(foo,gomer))$ 684check_zero_list(gomer)$ 685 686gomer : []$ 687for i : -100 thru 100 do ( 688 foo : f(12,-41.0, float(0.01 * i)), 689 e : part(foo,2), 690 foo : part(foo,1), 691 if (abs(foo) <= e) then foo : 0 else foo : 1, 692 gomer : cons(foo,gomer))$ 693check_zero_list(gomer)$ 694 695gomer : []$ 696for i : -100 thru 100 do ( 697 foo : f(12,-4.1, float(0.01 * i * %i)), 698 e : part(foo,2), 699 foo : part(foo,1), 700 if (abs(foo) <= e) then foo : 0 else foo : 1, 701 gomer : cons(foo,gomer))$ 702check_zero_list(gomer)$ 703 704test_name : "random jacobi float"$ 705 706gomer : [ ]$ 707for i : 1 thru 500 do ( 708 n : random(10), 709 a : random(11) / (1 + random(9)), 710 b : random(11) / (1 + random(9)), 711 x : (random(11) - 5) / (1 + random(10)), 712 exact : float(jacobi_p(n,a,b,x)), 713 approx : jacobi_p(n,float(a),float(b),float(x)), 714 if (abs(exact-part(approx,1)) <= part(approx,2)) then foo : 0 else foo : 1, 715 if (foo = 1) then print ("n = ",n, " a = ",a," b = ",b, "x = ",x), 716 gomer : cons(foo,gomer))$ 717check_zero_list(gomer)$ 718 719test_name : "random ultraspherical float"$ 720 721gomer : [ ]$ 722for i : 1 thru 500 do ( 723 n : random(10), 724 a : random(11) / (1 + random(9)), 725 x : (random(11) - 5) / (1 + random(10)), 726 exact : float(ultraspherical(n,a,x)), 727 approx : ultraspherical(n,float(a),float(x)), 728 if (abs(exact-part(approx,1)) <= part(approx,2)) then foo : 0 else foo : 1, 729 if (foo = 1) then print ("n = ",n, " a = ",a,"x = ",x), 730 gomer : cons(foo,gomer))$ 731check_zero_list(gomer)$ 732 733test_name : "random chebyshev_t float"$ 734gomer : [ ]$ 735for i : 1 thru 500 do ( 736 n : random(10), 737 x : (random(11) - 5) / (1 + random(10)), 738 exact : float(chebyshev_t(n,x)), 739 approx : chebyshev_t(n,float(x)), 740 if (abs(exact-part(approx,1)) <= part(approx,2)) then foo : 0 else foo : 1, 741 if (foo = 1) then print("chebyshev_t float error; n = ",n, "x = ",x), 742 gomer : cons(foo,gomer))$ 743check_zero_list(gomer)$ 744 745test_name : "random chebyshev_u float"$ 746gomer : [ ]$ 747for i : 1 thru 500 do ( 748 n : random(10), 749 x : (random(11) - 5) / (1 + random(10)), 750 exact : float(chebyshev_u(n,x)), 751 approx : chebyshev_u(n,float(x)), 752 if (abs(exact-part(approx,1)) <= part(approx,2)) then foo : 0 else foo : 1, 753 if (foo = 1) then print("chebyshev_t float error; n = ",n, "x = ",x), 754 gomer : cons(foo,gomer))$ 755check_zero_list(gomer)$ 756 757 758test_name : "random legendre float"$ 759gomer : [ ]$ 760for i : 1 thru 500 do ( 761 n : random(10), 762 x : (random(11) - 5) / (1 + random(10)), 763 exact : float(legendre_p(n,x)), 764 approx : legendre_p(n,float(x)), 765 if (abs(exact-part(approx,1)) <= part(approx,2)) then foo : 0 else foo : 1, 766 if (foo = 1) then print("legendre float error; n = ",n, "x = ",x), 767 gomer : cons(foo,gomer))$ 768check_zero_list(gomer)$ 769 770test_name : "random hermite float"$ 771gomer : [ ]$ 772for i : 1 thru 500 do ( 773 n : random(10), 774 x : (random(11) - 5) / (1 + random(10)), 775 exact : float(hermite(n,x)), 776 approx : hermite(n,float(x)), 777 if (abs(exact-part(approx,1)) <= part(approx,2)) then foo : 0 else foo : 1, 778 if (foo = 1) then print("hermite float error; n = ",n, "x = ",x), 779 gomer : cons(foo,gomer))$ 780check_zero_list(gomer)$ 781 782 783test_name : "random gen_laguerre float"$ 784 785gomer : [ ]$ 786for i : 1 thru 500 do ( 787 n : random(10), 788 a : random(11) / (1 + random(9)), 789 x : (random(11) - 5) / (1 + random(10)), 790 exact : float(gen_laguerre(n,a,x)), 791 approx : gen_laguerre(n,float(a),float(x)), 792 if (abs(exact-part(approx,1)) <= part(approx,2)) then foo : 0 else foo : 1, 793 if (foo = 1) then print ("n = ",n, " a = ",a," b = ",b, "x = ",x), 794 gomer : cons(foo,gomer))$ 795check_zero_list(gomer)$ 796 797test_name : "random assoc_legendre_p float"$ 798gomer : [ ]$ 799for i : 1 thru 500 do ( 800 n : random(10), 801 m : random(10), 802 x : (random(11) - 5)/ (1 + random(10)), 803 exact : float(assoc_legendre_p(n,m,x)), 804 approx : assoc_legendre_p(n,m,float(x)), 805 if not(intervalp(approx)) then approx : interval(approx,0), 806 if (abs(exact-part(approx,1)) <= part(approx,2)) then foo : 0 else foo : 1, 807 if (foo = 1) then print ("n =", n, " m = ", m, " x = ",x), 808 gomer : cons(foo,gomer))$ 809check_zero_list(gomer)$ 810 811remvalue(n)$ 812remvalue(x)$ 813remvalue(k)$ 814q : []$ 815 816test_name : "simple spherical_bessel_j test"$ 817 818fpprec : 75$ 819foo : ratsimp(spherical_bessel_j(5,x))$ 820q : makelist(ev(foo / bfloat(spherical_bessel_j(5,0.1 * k)) - 1.0b0, 821 x = 0.1b0 * k),k,1,10)$ 822q : map(lambda([s], if abs(s) < 5.0b-11 then 0 else 1), q)$ 823check_zero_list(q)$ 824 825foo : ratsimp(spherical_bessel_j(5,x))$ 826q : makelist(ev(foo / bfloat(spherical_bessel_j(5,0.1 * k)) - 1.0b0, 827 x = 0.1b0 * k),k,-10,-1)$ 828q : map(lambda([s], if abs(s) < 5.0b-11 then 0 else 1), q)$ 829check_zero_list(q)$ 830 831test_name : "simple spherical_bessel_y test"$ 832 833q : []$ 834foo : expand(spherical_bessel_y(5,x))$ 835q : makelist(ev(foo / bfloat(spherical_bessel_y(5,0.1 * k)) - 1.0b0, 836 x = 0.1b0*k),k,1,10)$ 837q : map(lambda([s], if abs(s) < 5.0b-11 then 0 else 1), q)$ 838check_zero_list(q)$ 839 840q : []$ 841foo : spherical_bessel_y(5,x)$ 842q : makelist(ev(foo / bfloat(spherical_bessel_y(5,0.1 * k)) - 1.0b0, 843 x = 0.1b0*k),k,-10,-1)$ 844q : map(lambda([s], if abs(s) < 5.0b-11 then 0 else 1), q)$ 845check_zero_list(q)$ 846 847 848remvalue(n)$ 849remvalue(i)$ 850remvalue(x)$ 851remvalue(a)$ 852remvalue(b)$ 853remvalue(gomer)$ 854remvalue(exact)$ 855remvalue(approx)$ 856remfunction(f)$ 857 858/* See A & S 8.5.4 859*/ 860 861test_name : "A&S 8.5.3"$ 862 863foo(n) := (n + 1)*legendre_q(n+1,x) - (2*n+1)*x*legendre_q(n,x) 864 + n*legendre_q(n-1,x)$ 865gomer : makelist(rat(foo(i)),i,1,15)$ 866check_zero_list(gomer)$ 867 868test_name : "A&S 8.5.3"$ 869foo(n) := (n + 1)*legendre_p(n+1,x) - (2*n+1)*x*legendre_p(n,x) 870 + n*legendre_p(n-1,x)$ 871gomer : makelist(rat(foo(i)),i,1,15)$ 872check_zero_list(gomer)$ 873 874test_name : "A&S 8.5.3"$ 875 876foo(n,m) := (n - m + 1) * assoc_legendre_q(n+1,m,x) - 877 (2*n + 1) * x * assoc_legendre_q(n,m,x) + (n + m) * assoc_legendre_q(n-1,m,x)$ 878 879gomer : makelist(makelist(foo(n,m),m,-n+1,n-1),n,1,5)$ 880gomer : map(ratsimp,apply(append, gomer))$ 881check_zero_list(gomer)$ 882 883 884 885/* See A&S 8.6.7 page 334 886*/ 887 888test_name : "A&S 8.6.7"$ 889foo : makelist(makelist(assoc_legendre_q(n,m,x) - (-1)^m * (1-x^2)^(m/2) * diff(legendre_q(n,x),x,m),m,0,n),n,0,4)$ 890foo : apply(append,foo)$ 891foo : ratsimp(foo)$ 892check_zero_list(foo)$ 893 894 895/* See G&R 8.810 page 1014 896*/ 897 898test_name : "G&R 8.810"$ 899foo : makelist(makelist(assoc_legendre_p(n,m,x) - (-1)^m * (1-x^2)^(m/2) * diff(legendre_p(n,x),x,m),n,0,5),m,0,5)$ 900foo : apply(append,foo)$ 901foo : rat(foo)$ 902check_zero_list(foo)$ 903 904 905/* See G&R 8.813 page 1015 906*/ 907 908test_name : "G&R 8.813 (1-6)"$ 909 910foo : [assoc_legendre_p(1,1,x) + sqrt(1-x^2), 911 assoc_legendre_p(2,1,x) + 3 * x * sqrt(1-x^2), 912 assoc_legendre_p(2,2,x) - 3 *(1-x^2), 913 assoc_legendre_p(3,1,x) + 3 * sqrt(1-x^2) *(5*x^2-1) / 2, 914 assoc_legendre_p(3,2,x) - 15*x*(1-x^2), 915 assoc_legendre_p(3,3,x) + 15 * (1-x^2)^(3/2)]$ 916 917foo : rat(foo)$ 918 919check_zero_list(foo)$ 920 921/* See G&R 8.950 (1) page 1033 922*/ 923 924test_name : "G&R 8.950 (1)"$ 925foo : makelist(hermite(n,x) - (-1)^n * exp(x^2) * diff(exp(-x^2),x,n),n,0,9)$ 926foo : rat(foo)$ 927check_zero_list(foo)$ 928 929 930/* See G&R 8.952 (1) page 1033 931*/ 932 933test_name : "G&R 8.952 (1)"$ 934foo : makelist(diff(hermite(n,x),x) - 2 * n * hermite(n-1,x),n,1,7)$ 935foo : rat(foo)$ 936check_zero_list(foo)$ 937 938/* See G&R 8.952 (2) page 1033 939*/ 940test_name : "G&R 8.952 (2)"$ 941foo : makelist(hermite(n+1,x) - 2 * x * hermite(n,x) + 2 * n * hermite(n-1,x),n,1,7)$ 942foo : rat(foo)$ 943check_zero_list(foo)$ 944 945/* See G&R 8.956 (1-3) page 1034 946*/ 947 948test_name : "G&R 8.956 (1-5)"$ 949 950foo : [hermite(0,x) - 1, 951 hermite(1,x) - 2*x, 952 hermite(2,x) - (4*x^2 -2), 953 hermite(3,x) - (8*x^3-12*x), 954 hermite(4,x) - (16*x^4 - 48 * x^2 + 12)]$ 955foo : rat(foo)$ 956check_zero_list(foo)$ 957 958/* See A&S 10.1.19 page 439 959*/ 960 961test_name : "A&S 10.1.19 spherical_hankel2"$ 962baz(n) := spherical_hankel2(n-1,x) + spherical_hankel2(n+1,x) - (2*n+1) * spherical_hankel2(n,x) /x$ 963foo : makelist(baz(k),k,-7,7)$ 964foo : rat(expand(foo))$ 965check_zero_list(foo)$ 966 967test_name : "A&S 10.1.19 spherical_hankel1"$ 968baz(n) := spherical_hankel1(n-1,x) + spherical_hankel1(n+1,x) - (2*n+1) * spherical_hankel1(n,x) /x$ 969foo : makelist(baz(k),k,-7,7)$ 970foo : rat(expand(foo))$ 971check_zero_list(foo)$ 972 973test_name : "A&S 10.1.19 spherical_bessel_j"$ 974baz(n) := spherical_bessel_j(n-1,x) + spherical_bessel_j(n+1,x) - (2*n+1) * spherical_bessel_j(n,x) /x$ 975foo : makelist(baz(k),k,-7,7)$ 976foo : rat(foo)$ 977check_zero_list(foo)$ 978 979test_name : "A&S 10.1.19 spherical_bessel_y"$ 980baz(n) := spherical_bessel_y(n-1,x) + spherical_bessel_y(n+1,x) - (2*n+1) * spherical_bessel_y(n,x) /x$ 981foo : makelist(baz(k),k,-7,7)$ 982foo : rat(foo)$ 983check_zero_list(foo)$ 984 985/*--------------*/ 986 987/* See A&S 10.1.20 page 439 988*/ 989 990test_name : "A&S 10.1.20 spherical_hankel1"$ 991remvalue(q)$ 992baz(n) := n * spherical_hankel1(n-1,q) - (n+1)*spherical_hankel1(n+1,q) 993 -(2*n+1) * diff(spherical_hankel1(n,q),q)$ 994foo : makelist(baz(k),k,-7,7)$ 995foo : ratsimp(foo)$ 996check_zero_list(foo)$ 997 998test_name : "A&S 10.1.20 spherical_hankel2"$ 999baz(n) := n * spherical_hankel2(n-1,q) - (n+1)*spherical_hankel2(n+1,q) 1000 -(2*n+1) * diff(spherical_hankel2(n,q),q)$ 1001foo : makelist(baz(k),k,-7,7)$ 1002foo : ratsimp(foo)$ 1003check_zero_list(foo)$ 1004 1005test_name : "A&S 10.1.20 spherical_bessel_j"$ 1006baz(n) := n * spherical_bessel_j(n-1,q) - (n+1)*spherical_bessel_j(n+1,q) 1007 -(2*n+1) * diff(spherical_bessel_j(n,q),q)$ 1008foo : makelist(baz(k),k,-7,7)$ 1009foo : ratsimp(foo)$ 1010check_zero_list(foo)$ 1011 1012test_name : "A&S 10.1.20 spherical_bessel_y"$ 1013baz(n) := n * spherical_bessel_y(n-1,q) - (n+1)*spherical_bessel_y(n+1,q) 1014 -(2*n+1) * diff(spherical_bessel_y(n,q),q)$ 1015foo : makelist(baz(k),k,-7,7)$ 1016foo : ratsimp(foo)$ 1017check_zero_list(foo)$ 1018 1019/*--------------*/ 1020 1021kill(labels)$ 1022 1023test_name : "A&S 10.1.31"$ 1024f(n) := spherical_bessel_j(n,t) *spherical_bessel_y(n-1,t) - 1025 spherical_bessel_j(n-1,t) *spherical_bessel_y(n,t) - 1/t^2$ 1026 1027foo : makelist(f(k),k,-3, 3)$ 1028foo : trigreduce(ratsimp(foo))$ 1029check_zero_list(foo)$ 1030 1031/*--------------*/ 1032 1033remove(n,integer)$ 1034remove(k,integer)$ 1035remove(i,integer)$ 1036remvalue(a,b,x,n,m,i,j)$ 1037test_name : "jacobi_p gradef test"$ 1038q : [ ]$ 1039for i : 0 thru 15 do ( 1040 foo : diff(jacobi_p(n,a,b,x^2) - jacobi_p(i,a,b,x^2),x), 1041 foo : ev(foo, n=i), 1042 foo : rat(foo), 1043 q : cons(foo,q))$ 1044check_zero_list(q)$ 1045 1046/*--------------*/ 1047 1048test_name : "ultraspherical gradef test"$ 1049q : [ ]$ 1050for i : 0 thru 15 do ( 1051 foo : diff(ultraspherical(n,a,x^2) - ultraspherical(i,a,x^2),x), 1052 foo : ev(foo, n = i), 1053 foo : rat(foo), 1054 q : cons(foo, q))$ 1055check_zero_list(q)$ 1056 1057/*--------------*/ 1058 1059test_name : "assoc_legendre_p gradef test"$ 1060remvalue(q,i,j,n,m,x,foo)$ 1061q : []$ 1062for i : 0 thru 15 do ( 1063 for j : -15 thru 15 do ( 1064 foo : diff(assoc_legendre_p(n,m,x) - assoc_legendre_p(i,j,x),x), 1065 foo : ev(foo, n=i,m=j), 1066 foo : radcan(foo), 1067 q : cons(foo,q)))$ 1068check_zero_list(q)$ 1069 1070/*--------------*/ 1071 1072test_name : "assoc_legendre_q gradef test"$ 1073remvalue(q,i,j,n,m,x,foo,w)$ 1074q : []$ 1075for i : 0 thru 5 do ( 1076 for j : -i thru i do ( 1077 w : assoc_legendre_q(i,j,x), 1078 w : diff(w,x), 1079 foo : diff(assoc_legendre_q(n, m, x),x) - w, 1080 foo : ev(foo, n=i,m=j), 1081 foo : ratsimp(expand(foo)), 1082 q : cons(foo,q)))$ 1083check_zero_list(q)$ 1084 1085/*--------------*/ 1086 1087test_name : "chebyshev_t gradef test"$ 1088q : []$ 1089for i : 0 thru 15 do ( 1090 foo : diff(chebyshev_t(n,x^2) - chebyshev_t(i,x^2),x), 1091 foo : ev(foo, n = i), 1092 foo : rat(foo), 1093 q : cons(foo,q))$ 1094check_zero_list(q)$ 1095 1096/*--------------*/ 1097 1098test_name : "chebyshev_u gradef test"$ 1099q : []$ 1100for i : 0 thru 15 do ( 1101 foo : diff(chebyshev_u(n,x^2) - chebyshev_u(i,x^2),x), 1102 foo : ev(foo, n=i), 1103 foo : rat(foo), 1104 q : cons(foo,q))$ 1105check_zero_list(q)$ 1106 1107/*--------------*/ 1108 1109test_name : "laguerre gradef test"$ 1110q : [ ]$ 1111for i : 0 thru 15 do ( 1112 foo : diff(laguerre(n,x^2) - laguerre(i,x^2),x), 1113 foo : ev(foo, n=i), 1114 foo : rat(foo), 1115 q : cons(foo,q))$ 1116check_zero_list(q)$ 1117 1118/*--------------*/ 1119 1120test_name : "gen_laguerre gradef test"$ 1121q : [ ]$ 1122for i : 0 thru 15 do ( 1123 foo : diff(gen_laguerre(n,a,x^2) - gen_laguerre(i,a,x^2),x), 1124 foo : ev(foo, n=i), 1125 foo : rat(foo), 1126 q : cons(foo,q))$ 1127check_zero_list(q)$ 1128 1129/*--------------*/ 1130 1131test_name : "legendre_p gradef test"$ 1132q : [ ]$ 1133for i : 0 thru 15 do ( 1134 foo : diff(legendre_p(n,x^2) - legendre_p(i,x^2),x), 1135 foo : ev(foo, n=i), 1136 foo : rat(foo), 1137 q : cons(foo,q))$ 1138check_zero_list(q)$ 1139 1140/*--------------*/ 1141 1142test_name : "legendre_q gradef test"$ 1143q : [ ]$ 1144for i : 0 thru 15 do ( 1145 foo : diff(legendre_q(n,x) - legendre_q(i,x),x), 1146 foo : ev(foo, n=i), 1147 foo : rat(foo), 1148 q : cons(foo,q))$ 1149check_zero_list(q)$ 1150 1151 1152test_name : "hermite gradef test"$ 1153q : []$ 1154for i : 0 thru 15 do ( 1155 foo : diff(hermite(n,x^2) - hermite(i,x^2),x), 1156 foo : ev(foo, n=i), 1157 foo : rat(foo), 1158 q : cons(foo,q))$ 1159check_zero_list(q)$ 1160 1161/*--------------*/ 1162 1163test_name : "spherical_hankel2 gradef test"$ 1164q : []$ 1165for i : 0 thru 15 do ( 1166 foo : diff(spherical_hankel2(n,x^2) - spherical_hankel2(i,x^2),x), 1167 foo : ev(foo, n=i), 1168 foo : rat(foo), 1169 q : cons(foo,q))$ 1170check_zero_list(q)$ 1171 1172/*--------------*/ 1173 1174test_name : "spherical_hankel1 gradef test"$ 1175q : [ ]$ 1176for i : 0 thru 15 do ( 1177 foo : diff(spherical_hankel1(n,x^2) - spherical_hankel1(i,x^2),x), 1178 foo : ev(foo, n=i), 1179 foo : ratsimp(foo), 1180 q : cons(foo,q))$ 1181check_zero_list(q)$ 1182 1183/*--------------*/ 1184 1185test_name : "spherical_bessel_j gradef test"$ 1186q : []$ 1187for i : 0 thru 15 do ( 1188 foo : diff(spherical_bessel_j(n,x^2) - spherical_bessel_j(i,x^2),x), 1189 foo : ev(foo, n=i), 1190 foo : ratsimp(foo), 1191 q : cons(foo,q))$ 1192check_zero_list(q)$ 1193 1194/*--------------*/ 1195 1196test_name : "spherical_bessel_y gradef test"$ 1197q : []$ 1198for i : 0 thru 15 do ( 1199 foo : diff(spherical_bessel_y(n,x^2) - spherical_bessel_y(i,x^2),x), 1200 foo : ev(foo, n=i), 1201 foo : rat(foo), 1202 q : cons(foo,q))$ 1203check_zero_list(q)$ 1204 1205 1206/*--------------*/ 1207 1208test_name : "spherical_harmonic gradef test"$ 1209remvalue(q,i,j,n,m,x,foo)$ 1210q : []$ 1211for i : 0 thru 5 do ( 1212 for j : -i thru i do ( 1213 foo : [diff(spherical_harmonic(n, m, x, y) - spherical_harmonic(i, j, x, y),x), 1214 diff(spherical_harmonic(n, m, x, y) - spherical_harmonic(i,j, x, y),y)], 1215 foo : ev(foo, n=i, m=j), 1216 foo : radcan(foo), 1217 foo : trigreduce(rat(foo)), 1218 q : append(foo,q)))$ 1219check_zero_list(q)$ 1220remvalue(q)$ 1221 1222/*--------------*/ 1223 1224test_name : "jacobi sum representation"$ 1225declare(n,integer)$ 1226foo : jacobi_p(n,p,q,2/3)$ 1227foo : ev(foo, sum, n=7)$ 1228foo : rat(foo - jacobi_p(7,p,q,2/3))$ 1229check_zero_list([foo])$ 1230 1231/*--------------*/ 1232 1233test_name : "ultraspherical sum representation"$ 1234foo : ultraspherical(n,p,-2/3)$ 1235foo : ev(foo, sum, n=2)$ 1236foo : rat(foo - ultraspherical(2,p,-2/3))$ 1237check_zero_list([foo])$ 1238 1239/*--------------*/ 1240test_name : "legendre_p sum representation"$ 1241foo : legendre_p(n,1/8)$ 1242foo : ev(foo, sum, n=8)$ 1243foo : rat(foo - legendre_p(8,1/8))$ 1244check_zero_list([foo])$ 1245 1246/*----------------*/ 1247 1248test_name : "chebyshev_t sum representation"$ 1249foo : chebyshev_t(n,2)$ 1250foo : ev(foo, sum, n=5)$ 1251foo : rat(foo - chebyshev_t(5,2))$ 1252check_zero_list([foo])$ 1253 1254/*----------------*/ 1255 1256test_name : "chebyshev_u sum representation"$ 1257foo : chebyshev_u(n,-1/4)$ 1258foo : ev(foo, sum, n=15)$ 1259foo : rat(foo - chebyshev_u(15,-1/4))$ 1260check_zero_list([foo])$ 1261 1262/*---------------*/ 1263 1264test_name : "laguerre sum representation"$ 1265foo : laguerre(n,2/3)$ 1266foo : ev(foo,sum,n=4)$ 1267foo : rat(foo - laguerre(4,2/3))$ 1268check_zero_list([foo])$ 1269 1270/*---------------*/ 1271 1272test_name : "generalized laguerre sum representation"$ 1273foo : gen_laguerre(n,a,-2/3)$ 1274foo : ev(foo,sum,n=4)$ 1275foo : rat(foo - gen_laguerre(4,a,-2/3))$ 1276check_zero_list([foo])$ 1277 1278 1279remvalue([x,n,k])$ 1280 1281test_name : "pochhammer-1"$ 1282check_zero_list(ratsimp([pochhammer(x,0) - 1, pochhammer(x,1) - x, 1283pochhammer(x,2) - x*(x+1), pochhammer(x,5)/pochhammer(x,4) - (x+4)]))$ 1284 1285test_name : "pochhammer-2"$ 1286check_zero_list(makelist(ratsimp(pochhammer(x,-k) * pochhammer(1-x,k) - (-1)^k),k,-5,5))$ 1287 1288test_name : "pochhammer-grad"$ 1289foo : pochhammer(x,n)$ 1290dfoo : diff(foo,x)$ 1291goober : makelist(diff(ev(foo,n = k),x) - ev(dfoo,n = k),k,-5,5)$ 1292check_zero_list(expand(subst([x = 7/2], goober))); 1293check_zero_list(expand(subst([x = -7/2], goober))); 1294 1295/*-----------------*/ 1296 1297test_name : "unit_step"$ 1298check_zero_list([unit_step(-2), unit_step(-1/9), unit_step(-1.2), 1299unit_step(-1.5b-2), unit_step(0), unit_step(0.0), unit_step(0.0b0), 1300unit_step(2)-1, unit_step(1/9)-1, unit_step(4.5)-1, unit_step(8.23b3)-1, 1301unit_step(x^2+1)-1,unit_step(exp(x))-1])$ 1302 1303/*----------------*/ 1304 1305test_name : "kron_delta"$ 1306 1307check_zero_list([kron_delta(1,2), kron_delta(1,1)-1, kron_delta(x,rat(x)) - 1, 1308kron_delta(x,y)-kron_delta(y,x), kron_delta(1.0,1.0)-1])$ 1309 1310 1311 1312/*----------------*/ 1313 1314test_name : "jacobi_p recursion"$ 1315 1316foo : orthopoly_recur(jacobi_p,[m,a,b,x])$ 1317foo : rhs(foo)-lhs(foo)$ 1318foo : makelist(ev(foo,m=k),k,1,6)$ 1319foo : rat(foo)$ 1320check_zero_list(foo)$ 1321 1322/*-----------------*/ 1323test_name : "ultraspherical recursion"$ 1324 1325foo : orthopoly_recur(ultraspherical,[m,a,x])$ 1326foo : rhs(foo)-lhs(foo)$ 1327foo : makelist(ev(foo,m=k),k,1,6)$ 1328foo : rat(foo)$ 1329check_zero_list(foo)$ 1330 1331/*-----------------*/ 1332test_name : "chebyshev_t_recursion"$ 1333 1334foo : orthopoly_recur(chebyshev_t,[m,x])$ 1335foo : rhs(foo)-lhs(foo)$ 1336foo : makelist(ev(foo,m=k),k,1,6)$ 1337foo : rat(foo)$ 1338check_zero_list(foo)$ 1339 1340/*-----------------*/ 1341test_name : "chebyshev_u recursion"$ 1342 1343foo : orthopoly_recur(chebyshev_u,[m,x])$ 1344foo : rhs(foo)-lhs(foo)$ 1345foo : makelist(ev(foo,m=k),k,1,6)$ 1346foo : rat(foo)$ 1347check_zero_list(foo)$ 1348 1349/*-----------------*/ 1350test_name : "legendre_p recursion"$ 1351 1352foo : orthopoly_recur(legendre_p,[m,x])$ 1353foo : rhs(foo)-lhs(foo)$ 1354foo : makelist(ev(foo,m=k),k,1,6)$ 1355foo : rat(foo)$ 1356check_zero_list(foo)$ 1357 1358/*-----------------*/ 1359test_name : "legendre_q recursion"$ 1360foo : orthopoly_recur(legendre_q,[m,x])$ 1361foo : rhs(foo)-lhs(foo)$ 1362foo : makelist(ev(foo,m=k),k,1,6)$ 1363foo : rat(foo)$ 1364check_zero_list(foo)$ 1365 1366/*-----------------*/ 1367test_name : "assoc_legendre_p recursion"$ 1368foo : orthopoly_recur(assoc_legendre_p,[m,1,x])$ 1369foo : rhs(foo)-lhs(foo)$ 1370foo : makelist(ev(foo,m=k),k,2,6)$ 1371foo : rat(foo)$ 1372check_zero_list(foo)$ 1373 1374/*-----------------*/ 1375test_name : "assoc_legendre_p recursion"$ 1376foo : orthopoly_recur(assoc_legendre_p,[m,2,x])$ 1377foo : rhs(foo)-lhs(foo)$ 1378foo : makelist(ev(foo,m=k),k,3,6)$ 1379foo : rat(foo)$ 1380check_zero_list(foo)$ 1381 1382/*-----------------*/ 1383test_name : "assoc_legendre_q recursion"$ 1384foo : orthopoly_recur(assoc_legendre_q,[m,1,x])$ 1385foo : rhs(foo)-lhs(foo)$ 1386foo : makelist(ev(foo,m=k),k,2,6)$ 1387foo : rat(foo)$ 1388check_zero_list(foo)$ 1389 1390/*----------------*/ 1391 1392test_name : "laguerre recursion"$ 1393 1394foo : orthopoly_recur(laguerre,[m,x])$ 1395foo : rhs(foo)-lhs(foo)$ 1396foo : makelist(ev(foo,m=k),k,1,6)$ 1397foo : rat(foo)$ 1398check_zero_list(foo)$ 1399 1400/*----------------*/ 1401 1402test_name : "gen_laguerre recursion"$ 1403 1404foo : orthopoly_recur(gen_laguerre,[m,a,x])$ 1405foo : rhs(foo)-lhs(foo)$ 1406foo : makelist(ev(foo,m=k),k,1,6)$ 1407foo : rat(foo)$ 1408check_zero_list(foo)$ 1409 1410/*-----------------*/ 1411test_name : "hermite recursion"$ 1412 1413foo : orthopoly_recur(hermite,[m,x])$ 1414foo : rhs(foo)-lhs(foo)$ 1415foo : makelist(ev(foo,m=k),k,1,6)$ 1416foo : rat(foo)$ 1417check_zero_list(foo)$ 1418 1419/*----------------*/ 1420test_name : "spherical_bessel_j recursion"$ 1421foo : orthopoly_recur(spherical_bessel_j,[m,x])$ 1422foo : rhs(foo)-lhs(foo)$ 1423foo : makelist(ev(foo,m=k),k,1,6)$ 1424foo : rat(foo)$ 1425check_zero_list(foo)$ 1426 1427 1428/*---------------------------------------*/ 1429test_name : "spherical_bessel_j recursion"$ 1430foo : orthopoly_recur(spherical_bessel_y,[m,x])$ 1431foo : rhs(foo)-lhs(foo)$ 1432foo : makelist(ev(foo,m=k),k,1,6)$ 1433foo : rat(foo)$ 1434check_zero_list(foo)$ 1435 1436/*---------------------------------------*/ 1437test_name : "spherical_hankel1 recursion"$ 1438foo : orthopoly_recur(spherical_hankel1,[m,x])$ 1439foo : rhs(foo)-lhs(foo)$ 1440foo : makelist(ev(foo,m=k),k,1,6)$ 1441foo : rat(foo)$ 1442check_zero_list(foo)$ 1443 1444/*---------------------------------------*/ 1445test_name : "spherical_hankel2 recursion"$ 1446foo : orthopoly_recur(spherical_hankel2,[m,x])$ 1447foo : rhs(foo)-lhs(foo)$ 1448foo : makelist(ev(foo,m=k),k,1,6)$ 1449foo : rat(foo)$ 1450check_zero_list(foo)$ 1451 1452/*---------------------------------------*/ 1453 1454test_name : "jacobi_weight"$ 1455foo : orthopoly_weight(jacobi_p,[n,2,3,x])$ 1456foo : integrate(jacobi_p(5,2,3,x) * jacobi_p(4,2,3,x) * foo[1],x,foo[2],foo[3])$ 1457check_zero_list([foo])$ 1458 1459/*---------------------------------------*/ 1460 1461test_name : "ultraspherical_weight"$ 1462foo : orthopoly_weight(ultraspherical,[n,2,x])$ 1463foo : integrate(ultraspherical(5,2,x) * ultraspherical(4,2,x) 1464 * foo[1],x,foo[2],foo[3])$ 1465check_zero_list([foo])$ 1466 1467/*---------------------------------------*/ 1468 1469test_name : "chebyshev_t_weight"$ 1470foo : orthopoly_weight(chebyshev_t,[n,x])$ 1471foo : integrate(chebyshev_t(5,x) * chebyshev_t(4,x) 1472 * foo[1],x,foo[2],foo[3])$ 1473check_zero_list([foo])$ 1474 1475/*---------------------------------------*/ 1476 1477test_name : "chebyshev_u_weight"$ 1478foo : orthopoly_weight(chebyshev_u,[n,x])$ 1479foo : integrate(chebyshev_u(5,x) * chebyshev_u(4,x) 1480 * foo[1],x,foo[2],foo[3])$ 1481check_zero_list([foo])$ 1482 1483/*---------------------------------------*/ 1484 1485test_name : "legendre_p_weight"$ 1486foo : orthopoly_weight(legendre_p,[n,x])$ 1487foo : integrate(legendre_p(5,x) * legendre_p(4,x) 1488 * foo[1],x,foo[2],foo[3])$ 1489check_zero_list([foo])$ 1490 1491/*---------------------------------------*/ 1492 1493test_name : "laguerre_weight"$ 1494foo : orthopoly_weight(laguerre,[n,x])$ 1495foo : integrate(laguerre(5,x) * laguerre(4,x) 1496 * foo[1],x,foo[2],foo[3])$ 1497check_zero_list([foo])$ 1498 1499/*---------------------------------------*/ 1500 1501test_name : "gen_laguerre_weight"$ 1502foo : orthopoly_weight(gen_laguerre,[n,1/2,x])$ 1503foo : integrate(gen_laguerre(5,1/2,x) * gen_laguerre(4,1/2,x) 1504 * foo[1],x,foo[2],foo[3])$ 1505check_zero_list([foo])$ 1506 1507/*---------------------------------------*/ 1508 1509test_name : "hermite_weight"$ 1510foo : orthopoly_weight(hermite,[n,x])$ 1511foo : integrate(hermite(5,x) * hermite(4,x) 1512 * foo[1],x,foo[2],foo[3])$ 1513check_zero_list([foo])$ 1514 1515 1516/*---------------------------------------*/ 1517 1518orthopoly_returns_intervals : false$ 1519 1520test_name : "legendre_p negative degree--symbolic argument"$ 1521foo1 : makelist (legendre_p (-k, u), k, 1, 8); 1522foo2 : makelist (legendre_p (k - 1, u), k, 1, 8); 1523check_zero_list(foo1 - foo2)$ 1524 1525test_name : "legendre_p negative degree--rational argument"$ 1526foo1 : makelist (legendre_p (-k, 11/7), k, 1, 8); 1527foo2 : makelist (legendre_p (k - 1, 11/7), k, 1, 8); 1528check_zero_list(foo1 - foo2)$ 1529 1530test_name : "legendre_p negative degree--float argument"$ 1531foo1 : makelist (legendre_p (-k, float (17/16)), k, 1, 8); 1532foo2 : makelist (legendre_p (k - 1, float (17/16)), k, 1, 8); 1533foo : map (lambda ([a, b], is (a = b)), foo1, foo2); 1534check_true_list (foo); 1535 1536/*---------------------------------------*/ 1537 1538test_name : "assoc_legendre_p negative degree--symbolic argument"$ 1539foo1 : makelist (assoc_legendre_p (-k, 1, u), k, 1, 8); 1540foo2 : makelist (assoc_legendre_p (k - 1, 1, u), k, 1, 8); 1541check_zero_list(foo1 - foo2)$ 1542 1543test_name : "assoc_legendre_p negative degree--rational argument"$ 1544foo1 : makelist (assoc_legendre_p (-k, 1, 11/7), k, 1, 8); 1545foo2 : makelist (assoc_legendre_p (k - 1, 1, 11/7), k, 1, 8); 1546check_zero_list(foo1 - foo2)$ 1547 1548test_name : "assoc_legendre_p negative degree--float argument"$ 1549foo1 : makelist (assoc_legendre_p (-k, 1, float (17/16)), k, 1, 8); 1550foo2 : makelist (assoc_legendre_p (k - 1, 1, float (17/16)), k, 1, 8); 1551foo : map (lambda ([a, b], is (a = b)), foo1, foo2); 1552check_true_list (foo); 1553 1554reset (orthopoly_returns_intervals); 1555 1556 1557print("orthopoly version = ", get('orthopoly,'version))$ 1558print("errors found = ", errors_found)$ 1559print("number of tests passed =", length(tests_pass))$ 1560 1561 1562/* Generate A&S Figure 22.4 page 776 1563*/ 1564 1565 1566orthopoly_returns_intervals : false$ 1567 1568foo : ['(ultraspherical(2,0.5,x)), 1569'(ultraspherical(3,0.5,x)), 1570'(ultraspherical(4,0.5,x)), 1571'(ultraspherical(5,0.5,x))]$ 1572 1573plot2d(foo,[x,-1,1])$ 1574 1575/* Generate A&S Figure 22.5 page 777 1576*/ 1577 1578 1579foo : ['(ultraspherical(5, 0.2, x)), 1580'(ultraspherical(5, 0.4, x)), 1581'(ultraspherical(5, 0.6, x)), 1582'(ultraspherical(5, 0.8, x)), 1583'(ultraspherical(5, 1.0, x))]$ 1584 1585plot2d(foo, [x,-0.8,0.8])$ 1586 1587/* Generate A&S Figure 22.6 page 778 1588*/ 1589 1590foo : ['(chebyshev_t(1,x)), 1591'(chebyshev_t(2,x)), 1592'(chebyshev_t(3,x)), 1593'(chebyshev_t(4,x)), 1594'(chebyshev_t(5,x))]$ 1595 1596 1597plot2d(foo,[x,-1.0,1.0])$ 1598 1599/* Generate A&S Figure 22.7 page 779 1600*/ 1601 1602 1603foo : ['(chebyshev_u(1,x)), 1604'(chebyshev_u(2,x)), 1605'(chebyshev_u(3,x)), 1606'(chebyshev_u(4,x)), 1607'(chebyshev_u(5,x))]$ 1608 1609 1610plot2d(foo,[x,-1.0,1.0])$ 1611 1612reset (float2bf); 1613[float2bf]; 1614