1COMMENT 2 3 THE REDUCE INTEGRATION TEST PACKAGE 4 5 Edited By 6 7 Anthony C. Hearn 8 The RAND Corporation 9 10 11This file is designed to provide a set of representative tests of the 12Reduce integration package. Not all examples go through, even when an 13integral exists, since some of the arguments are outside the domain of 14applicability of the current package. However, future improvements to 15the package will result in more closed-form evaluations in later 16releases. We would appreciate any additional contributions to this test 17file either because they illustrate some feature (good or bad) of the 18current package, or suggest domains which future versions should handle. 19Any suggestions for improved organization of this test file (e.g., in a 20way which corresponds more directly to the organization of a standard 21integration table book such as Gradshteyn and Ryznik) are welcome. 22 23Acknowledgments: 24 25The examples in this file have been contributed by the following. 26Any omissions to this list should be reported to the Editor. 27 28David M. Dahm 29James H. Davenport 30John P. Fitch 31Steven Harrington 32Anthony C. Hearn 33K. Siegfried Koelbig 34Ernst Krupnikov 35Arthur C. Norman 36Herbert Stoyan 37; 38 39COMMENT we first set up a suitable testing functions; 40 41fluid '(gcknt!*); 42 43global '(faillist!* gcnumber!* inittime number!-of!-integrals 44 unintlist!*); 45 46symbolic operator time; 47 48symbolic procedure initialize!-integral!-test; 49 begin 50 faillist!* := unintlist!* := nil; 51 number!-of!-integrals := 0; 52 gcnumber!* := gcknt!*; 53 inittime := time() 54 end; 55 56symbolic procedure summarize!-integral!-test; 57 begin scalar totaltime; 58 totaltime := time()-inittime; 59 prin2t 60 " ***** SUMMARY OF INTEGRAL TESTS *****"; 61 terpri(); 62 prin2 "Number of integrals tested: "; 63 prin2t number!-of!-integrals; 64 terpri(); 65 prin2 "Total time taken: "; 66 prin2 totaltime; 67 prin2t " ms"; 68 terpri(); 69 if gcnumber!* 70 then <<prin2 "Number of garbage collections: "; 71 prin2t (gcknt!* - gcnumber!*); 72 terpri()>>; 73 prin2 "Number of incorrect integrals: "; 74 prin2t length faillist!*; 75 terpri(); 76 prin2 "Number of unevaluated integrals: "; 77 prin2t length unintlist!*; 78 terpri(); 79 if faillist!* 80 then <<prin2t "Integrands of incorrect integrals are:"; 81 for each x in reverse faillist!* do mathprint car x>>; 82 if unintlist!* 83 then <<prin2t "Integrands of unevaluated integrals are:"; 84 terpri(); 85 for each x in reverse unintlist!* do mathprint car x>> 86 end; 87 88procedure testint(a,b); 89 begin scalar der,diffce,res,tt; 90 tt:=time(); 91 symbolic (number!-of!-integrals := number!-of!-integrals + 1); 92 res:=int(a,b); 93% write "time for integral: ",time()-tt," ms"; 94 off precise; 95 der := df(res,b); 96 diffce := der-a; 97 if diffce neq 0 98 then begin for all x let cot x=cos x/sin x, 99 sec x=1/cos x, 100 sin x**2=1-cos x**2, 101 tan(x/2)=sin x/(1+cos x), 102 tan x=sin x/cos x, 103 tanh x= 104 (e**(x)-e**(-x))/(e**x+e**(-x)), 105 coth x= 1/tanh x; 106 diffce := diffce; 107 for all x clear cot x,sec x,sin x**2,tan x,tan(x/2), 108 tanh x,coth x 109 end; 110 %hopefully, difference appeared non-zero due to absence of 111 %above transformations; 112 if diffce neq 0 113 then <<on combineexpt; diffce := diffce; off combineexpt>>; 114 if diffce neq 0 115 then begin scalar !*reduced; 116 symbolic(!*reduced := t); 117 for all x let cos(2x)= 1-2sin x**2, sin x**2=1-cos x**2; 118 diffce := diffce; 119 for all x clear cos(2x),sin x**2 120 end; 121 if diffce neq 0 122 then <<write 123 " ***** DERIVATIVE OF INTEGRAL NOT EQUAL TO INTEGRAND *****"; 124 symbolic(faillist!* := list(a,b,res,der) . faillist!*)>>; 125 symbolic if smemq('int,res) 126 then unintlist!* := list(a,b,res) . unintlist!*; 127 on precise; 128 return res 129 end; 130 131symbolic initialize!-integral!-test(); 132 133% References are to Gradshteyn and Ryznik. 134 135testint(1+x+x**2,x); 136testint(x**2*(2*x**2+x)**2,x); 137testint(x*(x**2+2*x+1),x); 138testint(1/x,x); % 2.01 #2 139testint((x+1)**3/(x-1)**4,x); 140testint(1/(x*(x-1)*(x+1)**2),x); 141testint((a*x+b)/((x-p)*(x-q)),x); 142testint(1/(a*x**2+b*x+c),x); 143testint((a*x+b)/(1+x**2),x); 144testint(1/(x**2-2*x+3),x); 145 146% Rational function examples from Hardy, Pure Mathematics, p 253 et seq. 147 148testint(1/((x-1)*(x**2+1))**2,x); 149testint(x/((x-a)*(x-b)*(x-c)),x); 150testint(x/((x**2+a**2)*(x**2+b**2)),x); 151testint(x**2/((x**2+a**2)*(x**2+b**2)),x); 152testint(x/((x-1)*(x**2+1)),x); 153testint(x/(1+x**3),x); 154testint(x**3/((x-1)**2*(x**3+1)),x); 155testint(1/(1+x**4),x); 156testint(x**2/(1+x**4),x); 157testint(1/(1+x**2+x**4),x); 158 159% Examples involving a+b*x. 160 161z := a+b*x; 162 163testint(z**p,x); 164testint(x*z**p,x); 165testint(x**2*z**p,x); 166testint(1/z,x); 167testint(1/z**2,x); 168testint(x/z,x); 169testint(x**2/z,x); 170testint(1/(x*z),x); 171testint(1/(x**2*z),x); 172testint(1/(x*z)**2,x); 173testint(1/(c**2+x**2),x); 174testint(1/(c**2-x**2),x); 175 176% More complicated rational function examples, mostly contributed 177% by David M. Dahm, who also developed the code to integrate them. 178 179testint(1/(2*x**3-1),x); 180testint(1/(x**3-2),x); 181testint(1/(a*x**3-b),x); 182testint(1/(x**4-2),x); 183testint(1/(5*x**4-1),x); 184testint(1/(3*x**4+7),x); 185testint(1/(x**4+3*x**2-1),x); 186testint(1/(x**4-3*x**2-1),x); 187testint(1/(x**4-3*x**2+1),x); 188testint(1/(x**4-4*x**2+1),x); 189testint(1/(x**4+4*x**2+1),x); 190testint(1/(x**4+x**2+2),x); 191testint(1/(x**4-x**2+2),x); 192testint(1/(x**6-1),x); 193testint(1/(x**6-2),x); 194testint(1/(x**6+2),x); 195testint(1/(x**8+1),x); 196testint(1/(x**8-1),x); 197testint(1/(x**8-x**4+1),x); 198testint(x**7/(x**12+1),x); 199 200% Examples involving logarithms. 201 202testint(log x,x); 203testint(x*log x,x); 204testint(x**2*log x,x); 205testint(x**p*log x,x); 206testint((log x)**2,x); 207testint(x**9*log x**11,x); 208testint(log x**2/x,x); 209testint(1/log x,x); 210testint(1/log(x+1),x); 211testint(1/(x*log x),x); 212testint(1/(x*log x)**2,x); 213testint((log x)**p/x,x); 214testint(log x *(a*x+b),x); 215testint((a*x+b)**2*log x,x); 216testint(log x/(a*x+b)**2,x); 217testint(x*log (a*x+b),x); 218testint(x**2*log(a*x+b),x); 219testint(log(x**2+a**2),x); 220testint(x*log(x**2+a**2),x); 221testint(x**2*log(x**2+a**2),x); 222testint(x**4*log(x**2+a**2),x); 223testint(log(x**2-a**2),x); 224testint(log(log(log(log(x)))),x); 225 226% Examples involving circular functions. 227 228testint(sin x,x); % 2.01 #5 229testint(cos x,x); % #6 230testint(tan x,x); % #11 231testint(1/tan(x),x); % 2.01 #12 232testint(1/(1+tan(x))**2,x); 233testint(1/cos x,x); 234testint(1/sin x,x); 235testint(sin x**2,x); 236testint(x**3*sin(x**2),x); 237testint(sin x**3,x); 238testint(sin x**p,x); 239testint((sin x**2+1)**2*cos x,x); 240testint(cos x**2,x); 241testint(cos x**3,x); 242testint(sin(a*x+b),x); 243testint(1/cos x**2,x); 244testint(sin x*sin(2*x),x); 245testint(x*sin x,x); 246testint(x**2*sin x,x); 247testint(x*sin x**2,x); 248testint(x**2*sin x**2,x); 249testint(x*sin x**3,x); 250testint(x*cos x,x); 251testint(x**2*cos x,x); 252testint(x*cos x**2,x); 253testint(x**2*cos x**2,x); 254testint(x*cos x**3,x); 255testint(sin x/x,x); 256testint(cos x/x,x); 257testint(sin x/x**2,x); 258testint(sin x**2/x,x); 259testint(tan x**3,x); 260% z := a+b*x; 261testint(sin z,x); 262testint(cos z,x); 263testint(tan z,x); 264testint(1/tan z,x); 265testint(1/sin z,x); 266testint(1/cos z,x); 267testint(sin z**2,x); 268testint(sin z**3,x); 269testint(cos z**2,x); 270testint(cos z**3,x); 271testint(1/cos z**2,x); 272testint(1/(1+cos x),x); 273testint(1/(1-cos x),x); 274testint(1/(1+sin x),x); 275testint(1/(1-sin x),x); 276testint(1/(a+b*sin x),x); 277testint(1/(a+b*sin x+cos x),x); 278testint(x**2*sin z**2,x); 279testint(cos x*cos(2*x),x); 280testint(x**2*cos z**2,x); 281testint(1/tan x**3,x); 282testint(x**3*tan(x)**4,x); 283testint(x**3*tan(x)**6,x); 284testint(x*tan(x)**2,x); 285testint(sin(2*x)*cos(3*x),x); 286testint(sin x**2*cos x**2,x); 287testint(1/(sin x**2*cos x**2),x); 288testint(d**x*sin x,x); 289testint(d**x*cos x,x); 290testint(x*d**x*sin x,x); 291testint(x*d**x*cos x,x); 292testint(x**2*d**x*sin x,x); 293testint(x**2*d**x*cos x,x); 294testint(x**3*d**x*sin x,x); 295testint(x**3*d**x*cos x,x); 296testint(sin x*sin(2*x)*sin(3*x),x); 297testint(cos x*cos(2*x)*cos(3*x),x); 298testint(sin(x*kx)**3*x**2,x); 299testint(x*cos(xi/sin(x))*cos(x)/sin(x)**2,x); 300 301% Mixed angles and half angles. 302 303int(cos(x)/(sin(x)*tan(x/2)),x); 304 305% This integral produces a messy result because the code for 306% converting half angle tans to sin and cos is not effective enough. 307 308testint(sin(a*x)/(b+c*sin(a*x))**2,x); 309 310% Examples involving logarithms and circular functions. 311 312testint(sin log x,x); 313testint(cos log x,x); 314 315% Examples involving exponentials. 316 317testint(e**x,x); % 2.01 #3 318testint(a**x,x); % 2.01 #4 319testint(e**(a*x),x); 320testint(e**(a*x)/x,x); 321testint(1/(a+b*e**(m*x)),x); 322testint(e**(2*x)/(1+e**x),x); 323testint(e**(2*x)*e**(a*x),x); 324testint(1/(a*e**(m*x)+b*e**(-m*x)),x); 325testint(x*e**(a*x),x); 326testint(x**20*e**x,x); 327testint(a**x/b**x,x); 328testint(a**x*b**x,x); 329testint(a**x/x**2,x); 330testint(x*a**x/(1+b*x)**2,x); 331testint(x*e**(a*x)/(1+a*x)**2,x); 332testint(x*k**(x**2),x); 333testint(e**(x**2),x); 334testint(2**(x**2),x); 335testint(2**(2*x**2),x); 336testint(e**(a*x**2),x); 337testint(e**(a**2*x**2),x); 338testint(x*e**(x**2),x); 339testint((x+1)*e**(1/x)/x**4,x); 340testint((2*x**3+x)*(e**(x**2))**2*e**(1-x*e**(x**2))/(1-x*e**(x**2))**2, 341 x); 342testint(e**(e**(e**(e**x))),x); 343 344% Examples involving exponentials and logarithms. 345 346testint(e**x*log x,x); 347testint(x*e**x*log x,x); 348testint(e**(2*x)*log(e**x),x); 349 350% Examples involving square roots. 351 352testint(sqrt(2)*x**2 + 2*x,x); 353testint(log x/sqrt(a*x+b),x); 354u:=sqrt(a+b*x); 355v:=sqrt(c+d*x); 356testint(u*v,x); 357testint(u,x); 358testint(x*u,x); 359testint(x**2*u,x); 360testint(u/x,x); 361testint(u/x**2,x); 362testint(1/u,x); 363testint(x/u,x); 364testint(x**2/u,x); 365testint(1/(x*u),x); 366testint(1/(x**2*u),x); 367testint(u**p,x); 368testint(x*u**p,x); 369testint(atan((-sqrt(2)+2*x)/sqrt(2)),x); 370testint(1/sqrt(x**2-1),x); 371testint(sqrt(x+1)*sqrt x,x); 372 373testint(sin(sqrt x),x); 374testint(x*(1-x^2)^(-9/4),x); 375testint(x/sqrt(1-x^4),x); 376testint(1/(x*sqrt(1+x^4)),x); 377testint(x/sqrt(1+x^2+x^4),x); 378testint(1/(x*sqrt(x^2-1-x^4)),x); 379 380% Examples from James Davenport's thesis: 381 382testint(1/sqrt(x**2-1)+10/sqrt(x**2-4),x); % p. 173 383testint(sqrt(x+sqrt(x**2+a**2))/x,x); 384 385% Examples generated by differentiating various functions. 386 387testint(df(sqrt(1+x**2)/(1-x),x),x); 388testint(df(log(x+sqrt(1+x**2)),x),x); 389testint(df(sqrt(x)+sqrt(x+1)+sqrt(x+2),x),x); 390testint(df(sqrt(x**5-2*x+1)-sqrt(x**3+1),x),x); 391 392% Another such example from James Davenport's thesis (p. 146). 393% It contains a point of order 3, which is found by use of Mazur's 394% bound on the torsion of elliptic curves over the rationals; 395 396testint(df(log(1+sqrt(x**3+1)),x),x); 397 398% Examples quoted by Joel Moses: 399 400testint(1/sqrt(2*h*r**2-alpha**2),r); 401testint(1/(r*sqrt(2*h*r**2-alpha**2-epsilon**2)),r); 402testint(1/(r*sqrt(2*h*r**2-alpha**2-2*k*r)),r); 403testint(1/(r*sqrt(2*h*r**2-alpha**2-epsilon**2-2*k*r)),r); 404testint(r/sqrt(2*e*r**2-alpha**2),r); 405testint(r/sqrt(2*e*r**2-alpha**2-epsilon**2),r); 406testint(r/sqrt(2*e*r**2-alpha**2-2*k*r**4),r); 407testint(r/sqrt(2*e*r**2-alpha**2-2*k*r),r); 408 409% These two integrals will evaluate, but they take a very long time 410% and the results are messy (compared with the algint results). 411 412% testint(1/(r*sqrt(2*h*r**2-alpha**2-2*k*r**4)),r); 413% testint(1/(r*sqrt(2*h*r**2-alpha**2-epsilon**2-2*k*r**4)),r); 414 415 416COMMENT many of these integrals used to require Steve Harrington's 417 code to evaluate. They originated in Novosibirsk as examples 418 of using Analytik. There are still a few examples that could 419 be evaluated using better heuristics; 420 421testint(a*sin(3*x+5)**2*cos(3*x+5),x); 422testint(log(x**2)/x**3,x); 423testint(x*sin(x+a),x); 424testint((log(x)*(1-x)-1)/(e**x*log(x)**2),x); 425testint(x**3*(a*x**2+b)**(-1),x); 426testint(x**(1/2)*(x+1)**(-7/2),x); 427testint(x**(-1)*(x+1)**(-1),x); 428testint(x**(-1/2)*(2*x-1)**(-1),x); 429testint((x**2+1)*x**(1/2),x); 430testint(x**(-1)*(x-a)**(1/3),x); 431testint(x*sinh(x),x); 432testint(x*cosh(x),x); 433testint(sinh(2*x)/cosh(2*x),x); 434testint((i*eps*sinh x-1)/(eps*i*cosh x+i*a-x),x); 435testint(sin(2*x+3)*cos(x)**2,x); 436testint(x*atan(x),x); 437testint(x*acot(x),x); 438testint(x*log(x**2+a),x); 439testint(sin(x+a)*cos(x),x); 440testint(cos(x+a)*sin(x),x); 441testint((1+sin(x))**(1/2),x); 442testint((1-sin(x))**(1/2),x); 443testint((1+cos(x))**(1/2),x); 444testint((1-cos(x))**(1/2),x); 445testint(1/(x**(1/2)-(x-1)**(1/2)),x); 446testint(1/(1-(x+1)**(1/2)),x); 447testint(x/(x**4+36)**(1/2),x); 448testint(1/(x**(1/3)+x**(1/2)),x); 449testint(log(2+3*x**2),x); 450testint(cot(x),x); 451testint(cot x**4,x); 452testint(tanh(x),x); 453testint(coth(x),x); 454testint(b**x,x); 455testint((x**4+x**(-4)+2)**(1/2),x); 456testint((2*x+1)/(3*x+2),x); 457testint(x*log(x+(x**2+1)**(1/2)),x); 458testint(x*(e**x*sin(x)+1)**2,x); 459testint(x*e**x*cos(x),x); 460 461COMMENT the following set came from Herbert Stoyan; 462 463testint(1/(x-3)**4,x); 464testint(x/(x**3-1),x); 465testint(x/(x**4-1),x); 466testint(log(x)*(x**3+1)/(x**4+2),x); 467testint(log(x)+log(x+1)+log(x+2),x); 468testint(1/(x**3+5),x); 469testint(1/sqrt(1+x**2),x); 470testint(sqrt(x**2+3),x); 471testint(x/(x+1)**2,x); 472 473COMMENT The following integrals were used among others as a test of 474 Moses' SIN program; 475 476testint(asin x,x); 477testint(x**2*asin x,x); 478testint(sec x**2/(1+sec x**2-3*tan x),x); 479testint(1/sec x**2,x); 480testint((5*x**2-3*x-2)/(x**2*(x-2)),x); 481testint(1/(4*x**2+9)**(1/2),x); 482testint((x**2+4)**(-1/2),x); 483testint(1/(9*x**2-12*x+10),x); 484testint(1/(x**8-2*x**7+2*x**6-2*x**5+x**4),x); 485testint((a*x**3+b*x**2+c*x+d)/((x+1)*x*(x-3)),x); 486testint(1/(2-log(x**2+1))**5,x); 487 488% The next integral appeared in Risch's 1968 paper. 489 490testint(2*x*e**(x**2)*log(x)+e**(x**2)/x+(log(x)-2)/(log(x)**2+x)**2+ 491 ((2/x)*log(x)+(1/x)+1)/(log(x)**2+x),x); 492 493% The following integral would not evaluate in REDUCE 3.3. 494 495testint(exp(x*ze+x/2)*sin(pi*ze)**4*x**4,ze); 496 497% This one evaluates: 498 499testint(erf(x),x); 500 501COMMENT So why not this one? 502 503Solved by adding a rule for the pattern matcher; 504 505testint(erf(x+a),x); 506 507COMMENT Three integrals added Sep. 2018, after improvements to the integrator; 508 509testint(atan(1/(1-x^2)),x); 510 511COMMENT The following integral can be done with extra substitution heuristics for exp; 512 513testint((1-e^x)^(1/2)/((2*e^x)-(e^(-x))-1),x); 514 515COMMENT The third one can be done easily using the substitution z=sin(x), but it currently fails; 516 517%testint((sin(x)-sin(2*x))/(2*cos(x)+((cos(x))^(1/2))),x); 518 519COMMENT Results are better if sin(2*x) is rewriten first; 520 521testint((sin(x)-2*sin(x)*cos(x))/(2*cos(x)+((cos(x))^(1/2))),x); 522 523 524COMMENT here is an example of using the integrator with pattern 525 matching; 526 527for all m,n let int(k1**m*log(k1)**n/(p**2-k1**2),k1)=foo(m,n), 528 int(k1*log(k1)**n/(p**2-k1**2),k1)=foo(1,n), 529 int(k1**m*log(k1)/(p**2-k1**2),k1)=foo(m,1), 530 int(k1*log(k1)/(p**2-k1**2),k1)=foo(1,1), 531 int(log(k1)**n/(k1*(p**2-k1**2)),k1)=foo(-1,n); 532 533int(k1**2*log(k1)/(p**2-k1**2),k1); 534 535COMMENT It is interesting to see how much of this one can be done; 536 537let f1s= (12*log(s/mc**2)*s**2*pi**2*mc**3*(-8*s-12*mc**2+3*mc) 538 + pi**2*(12*s**4*mc+3*s**4+176*s**3*mc**3-24*s**3*mc**2 539 -144*s**2*mc**5-48*s*mc**7+24*s*mc**6+4*mc**9-3*mc**8)) 540 /(384*e**(s/y)*s**2); 541 542int(f1s,s); 543 544factor Ei,log; 545 546ws; 547 548COMMENT the following is an example of integrals that used to loop 549 forever. They were first revealed by problems with Bessel 550 function integration when specfn was loaded, 551 e.g., int(x*besseli(2,x),x) or int(besselj(n,x),x); 552 553operator f; let {df(f(~x),x) => x*f(x-1)}; 554 555int(f x,x); 556 557 558COMMENT the following integrals reveal deficiencies in the current 559integrator; 560 561%high degree denominator; 562%testint(1/(2-log(x**2+1))**5,x); 563 564%this example should evaluate; 565testint(sin(2*x)/cos(x),x); 566 567%this example, which appeared in Tobey's thesis, needs factorization 568%over algebraic fields. It currently gives an ugly answer and so has 569%been suppressed; 570 571% testint((7*x**13+10*x**8+4*x**7-7*x**6-4*x**3-4*x**2+3*x+3)/ 572% (x**14-2*x**8-2*x**7-2*x**4-4*x**3-x**2+2*x+1),x); 573 574symbolic summarize!-integral!-test(); 575 576end; 577