share/orthopoly/test_orthopoly.mac

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