tests/rtestint.mac

   1 (kill(all),0);
   2 0$
   3
   4 /*
   5  * This is a set of test integrals for testing various parts of the
   6  * integration routines.
   7  */
   8
   9 /*
  10  * We're testing irinte now, which handles things like x^m*sqrt(R^(2*p+1))
  11  * for integral values of m and p.
  12  */
  13 R:c*x^2+b*x+a;
  14 c*x^2+b*x+a;
  15
  16 /* matchsum - derivative divides */
  17 integrate((4*x^3-x^2+4*x)^(-2/3)*(12*x^2-2*x+4),x);
  18 3*(4*x^3-x^2+4*x)^(1/3);
  19
  20 /*
  21  * Tests of den1den1, 1/x/sqrt(R)
  22  */
  23
  24 /* G&R a < 0, del < 0 */
  25 integrate(1/x/sqrt(x^2-1),x);
  26 -asin(1/abs(x));
  27
  28 /* G&R a > 0, del > 0 */
  29 integrate(1/x/sqrt(x^2+1),x);
  30 -asinh(1/abs(x));
  31
  32 /* G&R a > 0, del = 0 */
  33 /*integrate(1/x/sqrt(x^2-2*x+1),x);*/
  34
  35 /* G&R a < 0, del < 0 */
  36 integrate(1/x/sqrt(-6+5*x-x^2),x);
  37 asin(5*x/abs(x)-12/abs(x))/sqrt(6);
  38
  39 /* Last case */
  40 integrate(1/x/sqrt(-x^2+x-1),x);
  41 %i*asinh(x/(sqrt(3)*abs(x))-2/(sqrt(3)*abs(x)));
  42
  43
  44 /*
  45  * Test of denn, 1/sqrt(R^(2*p+1))
  46  */
  47
  48 (assume(4*a*c-b^2>0), assume(c > 0), assume(a > 0), assume(b > 0),
  49  assume(R > 0), 0);
  50 0;
  51
  52 /* 1/sqrt(R), repeated roots */
  53 integrate(1/sqrt(c*x^2+b*x+b^2/4/c),x);
  54 log(x+b/(2*c))/sqrt(c);
  55
  56 /* 1/sqrt(R^(2*p+1), p > 0, repeated roots */
  57 integrate(1/sqrt(c*x^2+b*x+b^2/4/c)^3,x);
  58 -1/(2*c^(3/2)*(x+b/(2*c))^2);
  59
  60
  61 /*
  62  * 1/sqrt(R), distinct roots,
  63  * G&R 2.261 case c > 0, del > 0
  64  */
  65 integrate(1/sqrt(R),x);
  66 1/sqrt(c)*asinh((2*c*x+b)/sqrt(4*a*c-b^2));
  67
  68 /*
  69  * 1/sqrt(R^3), distinct roots,
  70  * G&R 2.264 eq. 5
  71  */
  72 factor(integrate(1/R^(3/2),x));
  73 2*(2*c*x+b)/((4*a*c-b^2)*sqrt(c*x^2+b*x+a));
  74
  75 /*
  76  * General case.
  77  * 1/sqrt(R^(2*p+1),
  78  * G&R 2.263 eq. 3, for reduction; eq 4 for explicit solution.
  79  */
  80 factor(integrate(1/R^(5/2),x));
  81 2*(2*c*x+b)*(8*c^2*x^2+8*b*c*x+12*a*c-b^2)
  82          /(3*(4*a*c-b^2)^2*(c*x^2+b*x+a)^(3/2));
  83 /*
  84  * Tests of den1denn, 1/x/sqrt(R^(2*p+1))
  85  */
  86
  87 /*
  88  * G&R 2.268 gives the general solution for this:
  89  *
  90  * integrate(1/x/sqrt(R^(2*n+1)),x) =
  91  *   1/(2*n-1)/a/sqrt(R^(2*n-1))
  92  *   - b/2/a*integrate(1/sqrt(R^(2*n+1)),x)
  93  *   + 1/a*integrate(1/x/sqrt(R^(2*n-1)),x)
  94  *
  95  * For this example, n = 1
  96  *
  97  * integrate(1/x/sqrt(R^3),x) =
  98  *   1/a/sqrt(R)
  99  *   - b/2/a*integrate(1/sqrt(R^3),x)
 100  *   + 1/a*integrate(1/x/sqrt(R),x)
 101  *
 102  * and we already know the answers to the other two integrals.
 103  */
 104 integrate(1/x/sqrt(R^3),x);
 105 -asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))+2*a/(sqrt(4*a*c-b^2)*abs(x)))
 106         /a^(3/2)
 107         -2*b*c*x/(a*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))
 108         -b^2/(a*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))+1/(a*sqrt(c*x^2+b*x+a));
 109
 110 /*
 111  * Tests of nummdenn x^m/sqrt(R^(2*p+1))
 112  */
 113
 114 /*
 115  * See G&R 2.264, eq 2:
 116  *
 117  * integrate(x/sqrt(R),x) =
 118  *   sqrt(R)/c - b/(2*c)*integrate(1/sqrt(R),x)
 119  */
 120 integrate(x/sqrt(R),x);
 121 sqrt(''R)/c-b*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(2*c^(3/2));
 122
 123 /*
 124  * See G&R 2.264, eq 3:
 125  *
 126  * integrate(x^2/sqrt(R),x) =
 127  *   (x/(2*c)-3*b/(4*c^2))*sqrt(R)
 128  *     + (3*b^2/(8*c^2) - a/(2*c))*integrate(1/sqrt(R),x)
 129  */
 130 integrate(x^2/sqrt(R),x);
 131 -a*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(2*c^(3/2))
 132         +3*b^2*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(8*c^(5/2))
 133         +x*sqrt(c*x^2+b*x+a)/(2*c)-3*b*sqrt(c*x^2+b*x+a)/(4*c^2);
 134
 135
 136 /*
 137  * Test of repeated roots
 138  */
 139 integrate(x^2/sqrt(c*x^2+b*x+b^2/4/c),x);
 140 b^2*log(x+b/(2*c))/(4*c^(5/2))+x^2/(2*sqrt(c))-b*x/(2*c^(3/2));
 141
 142 integrate(x^2/sqrt(c*x^2+b*x+b^2/4/c)^3,x);
 143 log(x+b/(2*c))/c^(3/2)+b*x/(c^(5/2)*(x+b/(2*c))^2)
 144                               +3*b^2/(8*c^(7/2)*(x+b/(2*c))^2);
 145
 146 integrate(x^2/sqrt(c*x^2+b*x+b^2/4/c)^5,x);
 147 -1/(2*c^(5/2)*(x+b/(2*c))^2)+b/(3*c^(7/2)*(x+b/(2*c))^3)
 148                                     -b^2/(16*c^(9/2)*(x+b/(2*c))^4);
 149
 150 /*
 151  * Test denmnumn sqrt(R^(2*p+1))/x^m.
 152  */
 153
 154 integrate(sqrt(R)/x,x);
 155 -sqrt(a)*asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))
 156                         +2*a/(sqrt(4*a*c-b^2)*abs(x)))
 157          +b*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(2*sqrt(c))+sqrt(c*x^2+b*x+a);
 158
 159 integrate(sqrt(R^3)/x,x),radexpand:all;
 160 -a^(3/2)*asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))
 161                         +2*a/(sqrt(4*a*c-b^2)*abs(x)))
 162          +3*a*b*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(4*sqrt(c))
 163          -b^3*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(16*c^(3/2))
 164          +(c*x^2+b*x+a)^(3/2)/3+b*x*sqrt(c*x^2+b*x+a)/4
 165          +b^2*sqrt(c*x^2+b*x+a)/(8*c)+a*sqrt(c*x^2+b*x+a);
 166
 167 integrate(sqrt(R)/x^2,x);
 168 -b*asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))+2*a/(sqrt(4*a*c-b^2)*abs(x)))
 169          /(2*sqrt(a))
 170          +sqrt(c)*asinh((2*c*x+b)/sqrt(4*a*c-b^2))-sqrt(c*x^2+b*x+a)/x;
 171
 172 integrate(sqrt(R^3)/x^2,x),radexpand:all;
 173 -3*sqrt(a)*b
 174           *asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))+2*a/(sqrt(4*a*c-b^2)*abs(x)))
 175          /2
 176          +3*a*sqrt(c)*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/2
 177          +3*b^2*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(8*sqrt(c))
 178          -(c*x^2+b*x+a)^(3/2)/x+3*c*x*sqrt(c*x^2+b*x+a)/2
 179          +9*b*sqrt(c*x^2+b*x+a)/4;
 180
 181 /* Same as above, but test a = 0 (nonconstquadenum) */
 182 integrate(sqrt(R-a)/x,x);
 183 b*log(2*sqrt(c)*sqrt(c*x^2+b*x)+2*c*x+b)/(2*sqrt(c))+sqrt(c*x^2+b*x);
 184
 185 integrate(sqrt((R-a)^3)/x,x),radexpand:all;
 186 -b^3*log(2*sqrt(c)*sqrt(c*x^2+b*x)+2*c*x+b)/(16*c^(3/2))
 187          +(c*x^2+b*x)^(3/2)/3+b*x*sqrt(c*x^2+b*x)/4+b^2*sqrt(c*x^2+b*x)/(8*c);
 188
 189 integrate(sqrt(R-a)/x^2,x);
 190 sqrt(c)*log(2*sqrt(c)*sqrt(c*x^2+b*x)+2*c*x+b)-2*sqrt(c*x^2+b*x)/x;
 191
 192 integrate(sqrt((R-a)^3)/x^2,x),radexpand:all;
 193 3*b^2*log(2*sqrt(c)*sqrt(c*x^2+b*x)+2*c*x+b)/(8*sqrt(c))
 194          +(c*x^2+b*x)^(3/2)/(2*x)+3*b*sqrt(c*x^2+b*x)/4;
 195
 196 integrate(sqrt((R-a)^3)/x^3,x),radexpand:all;
 197 3*b*sqrt(c)*log(2*sqrt(c)*sqrt(c*x^2+b*x)+2*c*x+b)/2
 198          +(c*x^2+b*x)^(3/2)/x^2-3*b*sqrt(c*x^2+b*x)/x;
 199
 200 /*
 201  * Test denmdenn 1/sqrt(R^(2*p+1))/x^m
 202  */
 203 integrate(1/sqrt(R)/x^2,x);
 204 b*asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))+2*a/(sqrt(4*a*c-b^2)*abs(x)))
 205          /(2*a^(3/2))
 206          -sqrt(c*x^2+b*x+a)/(a*x);
 207
 208 integrate(1/sqrt(R)/x^3,x);
 209 c*asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))+2*a/(sqrt(4*a*c-b^2)*abs(x)))
 210          /(2*a^(3/2))
 211          -3*b^2
 212            *asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))+2*a/(sqrt(4*a*c-b^2)*abs(x)))
 213           /(8*a^(5/2))+3*b*sqrt(c*x^2+b*x+a)/(4*a^2*x)
 214          -sqrt(c*x^2+b*x+a)/(2*a*x^2);
 215
 216 integrate(1/sqrt(R^3)/x^2,x);
 217 3*b*asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))+2*a/(sqrt(4*a*c-b^2)*abs(x)))
 218          /(2*a^(5/2))
 219          -8*c^2*x/(a*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))
 220          +3*b^2*c*x/(a^2*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))
 221          -1/(a*x*sqrt(c*x^2+b*x+a))-4*b*c/(a*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))
 222          +3*b^3/(2*a^2*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))
 223          -3*b/(2*a^2*sqrt(c*x^2+b*x+a));
 224
 225 integrate(1/sqrt(R^3)/x^3,x);
 226 3*c*asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))+2*a/(sqrt(4*a*c-b^2)*abs(x)))
 227          /(2*a^(5/2))
 228          -15*b^2
 229             *asinh(b*x/(sqrt(4*a*c-b^2)*abs(x))+2*a/(sqrt(4*a*c-b^2)*abs(x)))
 230           /(8*a^(7/2))+13*b*c^2*x/(a^2*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))
 231          -15*b^3*c*x/(4*a^3*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))
 232          +5*b/(4*a^2*x*sqrt(c*x^2+b*x+a))-1/(2*a*x^2*sqrt(c*x^2+b*x+a))
 233          +13*b^2*c/(2*a^2*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))
 234          -15*b^4/(8*a^3*(4*a*c-b^2)*sqrt(c*x^2+b*x+a))
 235          -3*c/(2*a^2*sqrt(c*x^2+b*x+a))+15*b^2/(8*a^3*sqrt(c*x^2+b*x+a));
 236
 237 /* Same as above, but test the case of a = 0 (test noconstquad) */
 238 integrate(1/sqrt(R-a)/x^2,x);
 239 4*c*sqrt(c*x^2+b*x)/(3*b^2*x)-2*sqrt(c*x^2+b*x)/(3*b*x^2);
 240
 241 integrate(1/sqrt(R-a)/x^3,x);
 242 -16*c^2*sqrt(c*x^2+b*x)/(15*b^3*x)+8*c*sqrt(c*x^2+b*x)/(15*b^2*x^2)
 243                                           -2*sqrt(c*x^2+b*x)/(5*b*x^3);
 244
 245 integrate(1/sqrt((R-a)^3)/x^2,x),radexpand:all;
 246 -32*c^3*x/(5*b^4*sqrt(c*x^2+b*x))+4*c/(5*b^2*x*sqrt(c*x^2+b*x))
 247                                          -2/(5*b*x^2*sqrt(c*x^2+b*x))
 248                                          -16*c^2/(5*b^3*sqrt(c*x^2+b*x));
 249
 250 integrate(1/sqrt((R-a)^3)/x^3,x),radexpand:all;
 251 256*c^4*x/(35*b^5*sqrt(c*x^2+b*x))-32*c^2/(35*b^3*x*sqrt(c*x^2+b*x))
 252                                           +16*c/(35*b^2*x^2*sqrt(c*x^2+b*x))
 253                                           -2/(7*b*x^3*sqrt(c*x^2+b*x))
 254                                           +128*c^3/(35*b^4*sqrt(c*x^2+b*x));
 255
 256 /*
 257  * Test nummnumn x^m*sqrt(R^(2*p+1))
 258  */
 259
 260 integrate(sqrt(R),x);
 261 a*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(2*sqrt(c))
 262         -b^2*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(8*c^(3/2))
 263         +x*sqrt(c*x^2+b*x+a)/2+b*sqrt(c*x^2+b*x+a)/(4*c);
 264
 265 integrate(x*sqrt(R),x);
 266 -a*b*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(4*c^(3/2))
 267         +b^3*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(16*c^(5/2))
 268         +(c*x^2+b*x+a)^(3/2)/(3*c)-b*x*sqrt(c*x^2+b*x+a)/(4*c)
 269         -b^2*sqrt(c*x^2+b*x+a)/(8*c^2);
 270
 271 integrate(x^2*sqrt(R),x);
 272 a*(5*b^2-4*a*c)*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(32*c^(5/2))
 273         -b^2*(5*b^2-4*a*c)*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(128*c^(7/2))
 274         +x*(c*x^2+b*x+a)^(3/2)/(4*c)-5*b*(c*x^2+b*x+a)^(3/2)/(24*c^2)
 275         +(5*b^2-4*a*c)*x*sqrt(c*x^2+b*x+a)/(32*c^2)
 276         +b*(5*b^2-4*a*c)*sqrt(c*x^2+b*x+a)/(64*c^3);
 277
 278 integrate(x^3*sqrt(R),x);
 279 -7*a*b*(5*b^2-4*a*c)*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(320*c^(7/2))
 280         +7*b^3*(5*b^2-4*a*c)*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(1280*c^(9/2))
 281         +a^2*b*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(10*c^(5/2))
 282         -a*b^3*asinh((2*c*x+b)/sqrt(4*a*c-b^2))/(40*c^(7/2))
 283         +x^2*(c*x^2+b*x+a)^(3/2)/(5*c)-7*b*x*(c*x^2+b*x+a)^(3/2)/(40*c^2)
 284         -2*a*(c*x^2+b*x+a)^(3/2)/(15*c^2)+7*b^2*(c*x^2+b*x+a)^(3/2)/(48*c^3)
 285         -7*b*(5*b^2-4*a*c)*x*sqrt(c*x^2+b*x+a)/(320*c^3)
 286         +a*b*x*sqrt(c*x^2+b*x+a)/(10*c^2)
 287         -7*b^2*(5*b^2-4*a*c)*sqrt(c*x^2+b*x+a)/(640*c^4)
 288         +a*b^2*sqrt(c*x^2+b*x+a)/(20*c^3);
 289
 290 /* Integration is not correct: invalid 'false' term in results - ID: 3470668 */
 291 (assume(4*a*c > 1), 0);
 292 0;
 293
 294 integrate(x^4 / (a*x^2 + x + c)^(5/2), x);
 295 'integrate(x^4 / (a*x^2 + x + c)^(5/2), x);
 296
 297 (forget(4*a*c > 1), 0);
 298 0;
 299
 300 /*------------------------------------------------------------------------------*/
 301 /*
 302  * I (rtoy) have not verified any of these, except for the cases so marked.
 303  */
 304
 305 /*
 306  * Sample integrals from MIT LCS TR-092, Wang's thesis on definite integration
 307  */
 308
 309 /* p. 59 */
 310 /* Verified by doing the indefinite integral and taking limits */
 311 /* Tests mtoinf and zmtorat */
 312 factor(integrate((x^2+a1*x+b1)/(x^4+10*x^2+9),x,minf,inf));
 313 %pi*(b1+3)/12;
 314
 315 /* p. 62 */
 316 /*
 317  * Evaluating the indefinite integral and taking appropriate limits
 318  * gives the equivalent answer
 319  */
 320 /* Tests ztoinf and zmtorat */
 321 factor(integrate((x^2+a1*x+b1)/(x^4+10*x^2+9),x,0,inf));
 322 (%pi*b1+3*log(3)*a1+3*%pi)/24;
 323
 324 /*
 325  * p. 62
 326  *
 327  * Wang gives 2*%pi/3/sqrt(3), which is equivalent.  Why doesn't maxima
 328  * know that atan(sqrt(3)) is %pi/3?
 329  */
 330 /* Verified by doing evaluating indefinite integral */
 331 /* Tests ztoinf and zmtorat */
 332 integrate(1/(x^2+x+1),x,0,inf);
 333 2*sqrt(3)*atan(sqrt(3))/3;
 334
 335 /* p. 64 */
 336 /*
 337  * Using contour integration, we see this has poles at %i and 2*%i.  The
 338  * residues are 0 and -%i/3, respectively, so the integral is 2*%pi/3.
 339  */
 340 /* Tests ztoinf and zmtorat */
 341 integrate((x-%i)/((x-2*%i)*(x^2+1)),x,minf,inf);
 342 2*%pi/3;
 343
 344 /* p. 64 */
 345 /* Tests ztoinf and zmtorat */
 346 integrate((x-%i)/((x-2*%i)*(x^2+1)),x,0,inf);
 347 (2*%i*log(2)+2*%pi)/6;
 348
 349 /*
 350  * p. 67
 351  *
 352  * Wang gives 4*%pi/9/sqrt(3)-1/2, which is equivalent.
 353  */
 354
 355 /* Verified via indefinite integral */
 356 /* Tests ztoinf and zmtorat */
 357 integrate(1/(x^2+x+1)^3,x,0,inf);
 358 (8*sqrt(3)*%pi-27)/54;
 359
 360
 361 /* p. 70 */
 362 /* Verified via indefinite integral */
 363 /* Tests ztoinf and batapp */
 364 integrate(1/((1+x)*sqrt(x)),x,0,inf);
 365 %pi;
 366
 367 /* p. 71 */
 368 /*
 369  * First, we need assume(k < 0, k+1>0).  Then we also need to tell
 370  * maxima k is not an integer.  Then maxima can evaluate this integral.
 371  */
 372 /* Tests ztoinf and batapp */
 373 (assume(kj<0,kj+1>0,exp(kj)<1),0);
 374 0;
 375 (declare(kj,noninteger),0);
 376 0;
 377 integrate(x^kj/(x+3),x,0,inf);
 378 3^kj*beta(kj+1,-kj);
 379
 380 /* p. 74 */
 381 /* I think Wang meant the integral from a to inf, not 0 to inf. */
 382 /*
 383  * Tests parse-integrand, method-by-limits, dintegrate,
 384  * method-radical-poly, intsubs.
 385  */
 386 integrate(1/(x^2*sqrt(x^2-a^2)),x,a,inf);
 387 1/a^2;
 388
 389 /* p. 76 */
 390 /*
 391  * Wang gives 1024/6567561
 392  *
 393  * This particular integral matches form III with M = 6, N = 6, A = B
 394  * = C = 1.
 395  *
 396  * Based on his procedure, the case N>=M holds, and the integral is
 397  *
 398  * H*diff(U(1,1+z2,1+z3),z2,6)
 399  *
 400  * evaluated at z2=z3=0, where U(a,b,c) = 1/(sqrt(c)*(b/2+sqrt(a*c)))
 401  * and H = (-1)^N*sqrt(%pi)/2/gamma(3/2+N).  If we substitute the
 402  * values in, we find the answer is 2048/6567561.
 403  *
 404  * Thus, Wang's answer is wrong.
 405  */
 406 /*
 407  * Tests method-by-limits, dintegrate, method-radical-poly, intsubs.
 408  */
 409 integrate(x^6/(x^2+x+1)^(15/2),x,0,inf);
 410 2048/6567561;
 411
 412 /* p. 79 */
 413 /*
 414  * The single pole in the upper half plane is %i, and the residue of
 415  * exp(%i*x)/(x^2+1) is -%i/2*exp(-1), so the integral is as given.
 416  */
 417 /* Tests mtoinf, mtosc */
 418 integrate(cos(x)/(x^2+1),x,minf,inf);
 419 %e^(-1)*%pi;
 420
 421 /* p. 79 */
 422 /*
 423  * The single pole in the upper half plane is %i, and the residue of
 424  * x*exp(%i*x)/(x^2+1) is exp(-1)/2, so the integral is as given.
 425  */
 426 /* Tests mtoinf, mtosc */
 427 integrate(x*sin(x)/(x^2+1),x,minf,inf);
 428 %e^(-1)*%pi;
 429
 430 /* p. 79 */
 431 /* This follows from the above example.  Besides the function is odd. */
 432 integrate(x*cos(x)/(x^2+1),x,minf,inf);
 433 0;
 434
 435 /* p. 79 */
 436 /*
 437  * Wang gives
 438  * %pi/2*((-1/sqrt(3)-%i)*%e^((sqrt(3)-%i)/2) + (%i-1/sqrt(3))*%e^((-sqrt(3)-%i)/2))
 439  *
 440  * That can't be right because the result has real and imaginary parts.
 441  *
 442  * Maxima calculates this correctly.
 443  */
 444 /* Tests mtoinf, mtosc */
 445 integrate(x*cos(x)/(x^2+x+1),x,minf,inf);
 446 %e^-(sqrt(3)/2)*(3*%pi*sin(1/2)-sqrt(3)*%pi*cos(1/2))/3;
 447
 448 /* p. 80 */
 449 /* Tests mtoinf, mtosc */
 450 integrate(1/(%e^(%i*x)*(x^2+1)),x,minf,inf);
 451 %e^(-1)*%pi;
 452
 453 /* p. 80 */
 454 /*
 455  * Wang gives %pi/%e
 456  *
 457  * I think Wang's answer is wrong.
 458  *
 459  * If we exponentialize the integrand, we get
 460  *
 461  * %i*exp(-2*%i*x)/2/(x^2+1) - %i/2/(x^2+1).
 462  *
 463  * The integral of the second term is obviously %i*%pi/2.  The
 464  * integral of the first term can be down via residues, using the
 465  * residue at x = -%i, which is -exp(-2)/4.  Thus, the integral should
 466  * be %i*%pi*exp(-2)/2-%i*%pi/2.  (We need to change the sign of the
 467  * residue because we traverse the real axis in the opposite
 468  * direction.)
 469  */
 470 /* Tests mtoinf, mtosc */
 471 integrate(sin(x)/(%e^(%i*x)*(x^2+1)),x,minf,inf);
 472 -%e^(-2)*(%e^2-1)*%i*%pi/2;
 473
 474 /* p. 80 */
 475 /*
 476  * Wang gives %pi/%e^2
 477  *
 478  * I think Wang's answer is wrong.
 479  *
 480  * Using the same derivation as above, we get
 481  *
 482  * %i/2/(x^2+1)*(exp(-3*%i*x)-exp(-%i*x).
 483  *
 484  * We use the lower half plane and compute the residue at x = -%i,
 485  * which is exp(-3)*(exp(2)-1)/4.  Thus, the integral is
 486  * -%i*%pi*exp(-3)*(exp(2)-1)/2.
 487  */
 488 /* Tests mtoinf, mtosc */
 489 integrate(sin(x)/(%e^(2*%i*x)*(x^2+1)),x,minf,inf);
 490 -%e^(-3)*(%e^2-1)*%i*%pi/2;
 491
 492 /* Integration of products of sin(nz)/(nz) sometimes fails - ID: 3497046 */
 493 integrate(sin(z)/z*sin(3*z)/(3*z)*sin(5*z)/(5*z)*sin(7*z)/(7*z),z,minf,inf);
 494 44*%pi/315;
 495
 496 /* p. 83 */
 497 /*
 498  * Wang gives 3*gamma(3/7)*cos(3*%pi/14)/(7*9^(3/7))
 499  *
 500  * I've verified the derivation and the value, so I think this is right.
 501  */
 502 /* Tests ztoinf, scaxn */
 503 integrate(cos(9*x^(7/3)),x,0,inf);
 504 3*gamma(3/7)*cos(3*%pi/14)/(7*9^(3/7));
 505
 506
 507 /* p. 83 */
 508 /*
 509  * Wang gives 3*gamma(3/7)*sin(3*%pi/14)/(7*9^(3/7))
 510  *
 511  * Verified.
 512  */
 513 /* Tests ztoinf, scaxn */
 514 integrate(sin(9*x^(7/3)),x,0,inf);
 515 3*gamma(3/7)*sin(3*%pi/14)/(7*9^(3/7));
 516
 517
 518 /* p. 85 */
 519 /*
 520  * Wang gives gamma((2*%i+4)/3)/3/(%i/2)^((2*%i+4)/3)
 521  *
 522  * I think this is right, in a sense.  The final formula he gives is
 523  * wrong, but the derivation leading up to it is correct, except that
 524  * he has the wrong value for integrate(x^m*exp(k*x^n),x,0,inf).
 525  *
 526  * However, Wang also says these integrals only converge if
 527  * n-realpart(m)>1, and we have n=3 and realpart(m)=2, so we don't
 528  * satisfy the convergence criterion.
 529  */
 530 /*
 531 integrate(x^(3/2*%i+3)/exp(%i*x^3/2),x,0,inf);
 532 (-sqrt(3)*%i/2-1/2)*2^((2*%i+4)/3)*gamma((2*%i+4)/3)/3;
 533 */
 534
 535 /* p. 86 */
 536 /*
 537  * We can verify this by integrate(exp(%i*s*x)/sqrt(x),x,0,inf) and
 538  * taking the imaginary part.  This satisfies the convergence
 539  * criteria, and maxima returns the value
 540  * sqrt(%pi)/sqrt(s)*exp(%i*%pi/4).  Thus, we get the answer below.
 541  */
 542 /* Tests ztoinf, sconvert, ggr */
 543 (assume(s > 0), 0);
 544 0;
 545 ratsimp(integrate(sin(s*x)/sqrt(x),x,0,inf));
 546 sqrt(%pi)/(sqrt(2)*sqrt(s));
 547
 548 /* p. 87 */
 549 /* Tests ztoinf, ssp */
 550 (assume(r>0),0);
 551 0;
 552
 553 integrate(sin(r*x)/x,x,0,inf);
 554 %pi/2;
 555
 556 /* p. 87 */
 557 /* Tests ztoinf, ssp */
 558 (assume(q>0),0);
 559 0;
 560 integrate(sin(q*x)^2/x^2,x,0,inf);
 561 %pi*q/2;
 562
 563 /* p. 92 */
 564 /*
 565  * We can verify this integral by using
 566  * integrate(x^k*log(x)^2/(x^2+1),x,0,inf), and the evaluate the
 567  * result at k = 0.  The integral is a bit messy in terms of a beta
 568  * function, the digamma, and trigamma functions.  But maxima can
 569  * evaluate these at k = 0 and the result is %pi^3/8.
 570  *
 571  */
 572 /* Tests method-by-limits, zto1, batap-inf */
 573 integrate(log(x)^2/(x^2+1),x,0,inf);
 574 %pi^3/8;
 575
 576 /* p. 93 */
 577 /*
 578  * Wang gives
 579  *
 580  * log(3)*3^k*beta(k+1,-k) + 3^k*diff(beta(r,-k),r) - diff(beta(k+1,r),r)
 581  *
 582  * where the first derivative is evaluated at r=1+k and the second at r=-k.
 583  *
 584  * Actually, we can have maxima compute it for us, because we want the
 585  * derivative of x^k*beta(k+1,-k) which is
 586  *
 587  * log(3)*3^k*beta(k+1,-k)+3^k*(psi[0](k+1)-psi[0](-k))*beta(k+1,-k)
 588  *
 589  * Some simple numerical evaluations for k=-1/2 and k=-1/3 indicates
 590  * that this result is probably correct.
 591  */
 592 /* Tests method-by-limits, zto1, batap-inf */
 593 integrate(x^kj*log(x)/(x+3),x,0,inf);
 594 3^kj*(psi[0](kj+1)-psi[0](-kj))*beta(kj+1,-kj)+log(3)*3^kj*beta(kj+1,-kj);
 595
 596 /*
 597  * This is the same integral as above, using the substitution
 598  * y=log(x).  See also Bug 1451351.
 599  */
 600 /* Tests method-by-limits, zto1, batap-inf */
 601 integrate(x*exp((kj+1)*x)/(exp(x)+3),x,minf,inf);
 602 3^kj*(psi[0](kj+1)-psi[0](-kj))*beta(kj+1,-kj)+log(3)*3^kj*beta(kj+1,-kj);
 603
 604 /* p. 94
 605  *
 606  * subres/red gcd gets a "quotient is not exact" error.  spmod works better.
 607  *
 608  * The substitution y=1/x produces the same integral with a negative
 609  * sign, so the answer must be zero.
 610  */
 611 /* Tests method-by-limits, dintlog */
 612
 613 integrate((atan(x^(1/3)) + atan(x^(-1/3)))*log(x)/(x^2+1),x,0,inf);
 614 0;
 615
 616 /* p. 95 */
 617 /*
 618  * Wang gives %pi/sin(5*%pi/7)
 619  *
 620  * Based on his thesis, integration by parts for
 621  * integrate(log(1+x^(7/2))/x^2,0,inf) gives
 622  *
 623  *                  |inf
 624  * -log(1+x^(7/2))/x|    - integrate(-1/x*7/2*x^(5/2)/(1+x^(7/2)),x,0,inf)
 625  *                  |0
 626  *
 627  * = 7/2*integrate(x^(3/2)/(x^(7/2)+1),x,0,inf)
 628  *
 629  * Use the substitution (maxima doesn't do this) x = y^2 to get
 630  *
 631  * = 7*integrate(y^4/(y^7+1),y,0,inf)
 632  *
 633  * = %pi/sin(5*%pi/7)
 634  *
 635  * gcd:red gives a "quotient is not exact" error, so use
 636  * spmod.  We need intanalysis false to give the routine a chance.
 637  * (If true, we return the integral unevaluated.  I think this is
 638  * because solve fails to find the roots of x^(7/2)+1=0 correctly.)
 639  */
 640 /* Tests method-by-limits, dintlog, dintbypart, ztoinf, batapp */
 641 (intanalysis:false, 0);
 642 0;
 643
 644 integrate(log(x^(7/2)+1)/x^2,x,0,inf);
 645 /*%pi/sin(5*%pi/7);*/
 646 /* The result of the integral is beta(5/7,2/7). This simplifies to
 647    %pi/sin(5*%pi/z). Because we have introduced a symmetry rule for
 648    the beta function, Maxima first simplifies to beta(2/7,5/7).
 649    The equivalent result is %pi/sin(2*%pi/7). 03/2009 (DK) */
 650 %pi/sin(2*%pi/7);
 651
 652 (intanalysis:true, 0);
 653 0;
 654
 655 /* p. 97 */
 656 /* Verified via indefinite integral */
 657 /* Tests dintegrate, method-by-limits, method-radical-poly */
 658 factor(integrate(1/(exp(x/3)*(5/exp(x/3)+7)),x,0,inf));
 659 3*(log(12)-log(7))/5;
 660
 661 /* p. 97 */
 662 /* Verified via indefinite integral */
 663 integrate(exp(x/4)/(9*exp(x/2)+4),x,minf,inf);
 664 %pi/3;
 665
 666 /* p. 99 */
 667 /* Wang says %pi
 668  *
 669  * Not sure if this is right or not.  Maxima is using rectzto%pi2 to
 670  * evaluate this integral.
 671  *
 672  * If we follow Wang's derivation, we need to find q(z) such that
 673  * q(z)-q(z+2*%i*%pi) = z.  We find that q(z) = %i/(4*%pi)*z^2+z/2.
 674  *
 675  * R(z) = 1/((z-1/z)/2-%i.  We need to find the poles of R(z) such
 676  * that 0 <= Im(z) < 2*%pi.  In this case, the (only) pole is z = %i.
 677  * We need to then find the residue of q(z)/(sinh(z)-%i) at z =
 678  * glog(%i) = %i*%pi/2.  Maxima says the residue is -%i/2.  Hence the
 679  * value of the integral is 2*%i*%pi*(-%i/2) = %pi.
 680  */
 681 /* Tests mtoinf, rectzto%pi2 */
 682 integrate(x/(sinh(x)-%i),x,minf,inf);
 683 %pi;
 684
 685
 686 /* p. 101 */
 687 /* Tests ztoinf, ggr */
 688 integrate(exp(-x^2),x,0,inf);
 689 sqrt(%pi)/2;
 690
 691 /* p. 102 */
 692 /*
 693  * The substitution y=s*t produces 1/sqrt(s)*integrate(y^(-1/2)*exp(-y),y,0,inf),
 694  * which is gamma(1/2)/sqrt(s).
 695  */
 696 /* Tests ztoinf, ggr */
 697 integrate(1/sqrt(t)/exp(s*t),t,0,inf);
 698 sqrt(%pi)/sqrt(s);
 699
 700 /*
 701  * Write this as imagpart(integrate(exp(%i*s*x)/exp(x),x)).  Evaluate
 702  * the result at x=inf and at x=0.  We get the expected answer.
 703  */
 704
 705 integrate(sin(s*x)/exp(x),x,0,inf);
 706 s/(s^2+1);
 707
 708 /* p. 103 */
 709 /*
 710  * There's a typo in Wang's thesis.  He has the integrand as
 711  * exp(-s*t)*x^(1/3)*log(x).  I changed the x's to t's.
 712  *
 713  * Also, Wang says the answer is
 714  *
 715  * -gamma(1/3)*(log(s)-psi(4/3))/(6*s^(4/3))
 716  *
 717  * Assuming psi(4/3) is maxima's digamma function, psi[0](4/3), this
 718  * result is 1/2 of the result maxima produces.
 719  *
 720  * Consider integrate(t^(z+1/3)*exp(-s*t),t,0,inf).  Differentiating
 721  * this wrt to z gives us the integral we're looking for, when we
 722  * evaluate the result at z = 0.
 723  *
 724  * But this integral is equal to
 725  * 1/s^(z+4/3)*integrate(y^(z+1/3)*exp(-y),y,0,inf), which is a Gamma
 726  * function.  Thus, we have
 727  *
 728  * gamma(z+4/3)/s^(z+4/3)
 729  *
 730  * Differentiate wrt to z:
 731  *
 732  * s^(-z-4/3)*(psi[0](z+4/3)-log(s))*gamma(z+4/3)
 733  *
 734  * Evaluate at z = 0:
 735  *
 736  * gamma(1/3)*(psi[0](4/3)-log(s))/(3*s^(4/3))
 737  *
 738  * This is the answer we get below, once we evaluate psi[0](4/3).
 739  * Wang's thesis is off by a factor of 2.
 740  */
 741 integrate(exp(-s*t)*t^(1/3)*log(t),t,0,inf);
 742 gamma(1/3)*(-3*log(3)/2-%pi/(2*sqrt(3))-%gamma+3)/(3*s^(4/3))
 743         -gamma(1/3)*log(s)/(3*s^(4/3));
 744
 745 /* p. 106 */
 746 /* Verified via indefinite integral */
 747 ratsimp(integrate(1/(x^2-3),x,0,1) - sqrt(3)*log(2-sqrt(3))/6);
 748 0$
 749
 750 /* p. 109 */
 751 /* Verified via indefinite integral */
 752 integrate(cos(x)^2-sin(x),x,0,2*%pi);
 753 %pi;
 754
 755 /* p. 109 */
 756 /*
 757  * Wang says 2*%pi
 758  *
 759  * Evaluation of the indefinite integral says 0.
 760  */
 761 integrate(exp(2*%i*x)/(exp(%i*x)+3),x,0,2*%pi);
 762 0;
 763
 764 /* p. 110 */
 765 /*
 766  * Wang gives the equivalent 6/5*gamma(2/3)*gamma(3/4)/gamma(5/12)
 767  */
 768 integrate(cos(x)^(1/3)*sin(x)^(1/2),x,0,%pi/2);
 769 beta(3/4,2/3)/2;
 770
 771 /* p. 112 */
 772 /*
 773  * Wang gives the equivalent 18/7*sqrt(%pi)*gamma(2/3)/gamma(1/6)
 774  */
 775 integrate(cos(x)^(1/3)*sin(x)^2,x,-%pi/2,%pi/2);
 776 beta(2/3,3/2);
 777
 778 /* p. 112 */
 779 integrate(cos(x)^3*sin(x)^2,x,3*%pi/2,3*%pi);
 780 2/15;
 781
 782 /* p. 114 */
 783 /*
 784  * Wang gives -%i/3*plog((3+%i*sqrt(7))/4)
 785  *
 786  * The rectform is -atan(sqrt(7)/3)/3, which is -acos(3/4)/3,
 787  * which is -1/3*(%pi/2-asin(3/4))=-%pi/6+asin(3/4)/3.
 788  * Maxima gives an answer with the opposite sign.
 789  *
 790  * But we can evaluate the indefinite integral to be
 791  *
 792  *   -asin(3/x)/3.
 793  *
 794  * Hence the answer should be %pi/6-asin(3/4)/3.
 795  *
 796  * Wang's thesis is wrong.
 797  */
 798 integrate(1/(x*sqrt(x^2-9)),x,3,4);
 799 %pi/6-asin(3/4)/3;
 800
 801
 802 /* p. 114 */
 803 /* Verified via indefinite integral */
 804 logcontract(integrate(1/((x+1)*sqrt(4-x^2)),x,0,2) - sqrt(3)*log(sqrt(3)+2)/3);
 805 0$
 806
 807 /* p. 120 */
 808 /*
 809  * Maxima asks if k is an integer, but it doesn't matter,
 810  * and I don't know how to tell maxima that k isn't an integer.
 811  *
 812  */
 813 (assume(k1>0),declare(k1,noninteger),0);
 814 0;
 815 /*
 816  * This is not such a good integral to use.  log(x) is negative over
 817  * the integration limits, so integrand is complex.  Let's change the
 818  * integrand to be (-log(x))^k1, instead, which is real-valued.
 819  *
 820  * In this case, it's easy to see that the substitution y=-log(x)
 821  * gives as a gamma function.
 822  */
 823 /*
 824 integrate(log(x)^k1,x,0,1);
 825 (-1)^k1*gamma(k1+1);
 826 */
 827 integrate((-log(x))^k1,x,0,1);
 828 k1*gamma(k1);
 829
 830 /* p. 120 */
 831 integrate(sqrt(log(x))/x^2,x,1,inf);
 832 sqrt(%pi)/2;
 833
 834 /* p. 120 */
 835 /*
 836  * Wang says -sqrt(%pi)/sqrt(4/3) = -sqrt(3)*sqrt(%pi)/2.
 837  *
 838  * This is wrong.  Clearly, the integrand is positive, so the integral
 839  * must be positive.
 840  *
 841  * The issue is that radexpand is true so the Risch integrator is
 842  * doing something not quite right and getting the sign wrong.  If we
 843  * set radexpand to false, we get the correct answer.
 844  */
 845 radexpand:false;
 846 false;
 847
 848 integrate(x^(1/3)/sqrt(-log(x)),x,0,1);
 849 sqrt(3)*sqrt(%pi)/2;
 850
 851 domain:complex;
 852 complex;
 853
 854 /* [ 1731624 ] asked about sign of yx in integral containing only z */
 855 integrate(exp(sqrt(x^3)),x,0,1);
 856 2*((-1)^(1/3)*gamma_incomplete(2/3,-1)-(-1)^(1/3)*gamma(2/3))/3;
 857 /* with radexpand:true and domain:real we get
 858    integrate(exp(sqrt(x^3)),x)  ->  -2*gamma_incomplete(2/3,-x^(3/2))/3
 859    (which is not correct) due to invalid simplification in
 860    integrate-exp-special
 861 */
 862
 863 domain:real;
 864 real;
 865
 866 radexpand:true;
 867 true;
 868
 869 /* p. 121 */
 870 /*
 871  * Wang gives 4*diff(beta(r,3/2),r)
 872  * evaluated at r = 1.
 873  *
 874  * Maxima says 4*factor(ratsimp(subst(1,r,diff(makegamma(beta(r,3/2)),r))))
 875  * = 16*(3*log(2)-4)/9.  Good.
 876  */
 877 factor(integrate(sqrt(1-sqrt(x))*log(x)/sqrt(x),x,0,1));
 878 16*(3*log(2)-4)/9;
 879
 880 /*------------------------------------------------------------------------------*/
 881
 882 /* Problems from Moses' thesis, Appendix D, problems proposed by McIntosh */
 883
 884 /* Problem 1 */
 885 /*
 886  * Answer = 1/a*asin(a/r/sqrt(2*h))
 887  *
 888  * (This answer seems wrong.  The derivative has the wrong sign?)
 889  */
 890 (assume(h>0),0);
 891 0;
 892
 893 ratsimp(integrate(1/r/sqrt(2*h*r^2-a^2),r) + asin(a/sqrt(2)/sqrt(h)/abs(r))/a);
 894 0$
 895
 896 /* Problem 2 */
 897 /*
 898  * Answer = -1/sqrt(a^2-eps^2)*asin(sqrt(a^2+eps^2)/sqrt(2*h)/r)
 899  */
 900 /*
 901 integrate(1/sqrt(2*h*r^2-a^2-eps^2),r);
 902 */
 903
 904 /* Problem 3 */
 905 /*
 906  * Answer = 1/2/a*asin((h*r^2-a^2)/r^2/sqrt(h^2-2*k*a^2))
 907  * h^2 > 2*a^2*k
 908  *
 909  */
 910 (assume(h^2>2*a^2*k),0);
 911 0;
 912 factor(integrate(1/r/sqrt(2*h*r^2-a^2-2*k*r^4),r));
 913 1/2/a*asin((h*r^2-a^2)/r^2/sqrt(h^2-2*k*a^2));
 914
 915 /* Problem 4 */
 916 /*
 917  * Answer = 1/2/sqrt(a^2+eps^2)*asin((h*r^2-(a^2+eps^2))/r^2/sqrt(h^2-2*k*(a^2+eps^2)))
 918  * h^2 > 2*(a^2+eps^2)*k
 919  *
 920  * Maxima doesn't get stuck anymore, but I can't figure out the magic
 921  * assume expression so that integrate doesn't ask questions anymore.
 922  */
 923 /*
 924 (assume(2*eps^2*k+2*a^2*k-h^2 < 0),0);
 925 0;
 926 integrate(1/r/sqrt(2*h*r^2-a^2-eps^2-2*k*r^4),r);
 927 1/2/sqrt(a^2+eps^2)*asin((h*r^2-(a^2+eps^2))/r^2/sqrt(h^2-2*k*(a^2+eps^2)));
 928 */
 929
 930 /* Problem 5 */
 931 /*
 932  * Answer = 1/a*asin((k*r-a^2)/r/sqrt(k^2+2*h*a^2))
 933  * k^2+2*h*a^2 > 0
 934  *
 935  * This answer seems wrong.  The derivative doesn't match.
 936  *
 937  * Maxima's answer has a derivative that matches.
 938  */
 939 factor(integrate(1/r/sqrt(2*h*r^2-a^2-2*k*r),r));
 940 -asin((k*r+a^2)/(sqrt(k^2+2*a^2*h)*r))/a;
 941
 942 /* Problem 6 */
 943 /*
 944  * Answer = 1/(a^2+eps^2)*asin((k*r-(a^2+eps^2))/r/sqrt(k^2+2*h*(a^2+eps^2)))
 945  * k^2+2*h*a^2 > 0
 946  *
 947  * Same issue as problem 5.
 948  *
 949  * Maxima's answer has a derivative that matches.
 950  */
 951
 952 factor(integrate(1/r/sqrt(2*h*r^2-a^2-eps^2-2*k*r),r));
 953 -asin((k*r+eps^2+a^2)/(sqrt(k^2+(2*eps^2+2*a^2)*h)*r))/sqrt(eps^2+a^2);
 954
 955 /* Problem 7 */
 956 integrate(x/sqrt(2*e*x^2-a^2),x);
 957 sqrt(2*e*x^2-a^2)/(2*e);
 958
 959 /* Problem 8 */
 960 integrate(x/sqrt(2*e*x^2-a^2-eps^2),x);
 961 sqrt(2*e*x^2-eps^2-a^2)/(2*e);
 962
 963 /* Problem 9 */
 964 /*
 965  * Answer is 1/2/sqrt(2*k)*asin((2*k*r^2-e)/sqrt(e^2-2*k*a^2))
 966  * e^2>2*k*a^2
 967  *
 968  */
 969 (assume(e^2>2*k*a^2,k>0,e>0),0);
 970 0;
 971 factor(integrate(r/sqrt(2*e*r^2-a^2-2*k*r^4),r));
 972 1/2/sqrt(2*k)*asin((2*k*r^2-e)/sqrt(e^2-2*k*a^2));
 973
 974 /* Problem 10 */
 975 /*
 976  * Answer is 1/2/sqrt(2*k)*asin((2*k*r^2-e)/sqrt(e^2-2*k*(a^2+eps^2)))
 977  * e^2>2*k*a^2
 978  *
 979  * For some reason, maxima returns
 980  * -%i*asinh((2*k*r^2-e)/sqrt((2*eps^2+2*a^2)*k-e^2))/(2*sqrt(2)*sqrt(k))
 981  *
 982  * Using the fact that asin(%i*x) = %i*asinh(x), we see that the
 983  * answer is identical.
 984  */
 985 (assume(2*eps^2*k+2*a^2*k-e^2>0),0);
 986 0;
 987 factor(integrate(r/sqrt(2*e*r^2-a^2-eps^2-2*k*r^4),r));
 988 -%i*asinh((2*k*r^2-e)/sqrt((2*eps^2+2*a^2)*k-e^2))/(2*sqrt(2)*sqrt(k));
 989
 990
 991 /* Problem 11 */
 992 /*
 993  * Answer is sqrt(2*E*r^2-a^2-2*k*r)/2*E
 994  *            + 1/(2*h*E*sqrt(-2*E))*asin((2*E*r+k)/sqrt(k^2-2*E-a^2)
 995  * E < 0
 996  */
 997 (assume(E<0,k^2+2*E*a^2>0,k>0),0);
 998 0;
 999 integrate(r/sqrt(2*E*r^2-a^2-2*k*r),r);
1000 sqrt(2*E*r^2-2*k*r-a^2)/(2*E)-k*asin((4*E*r-2*k)/sqrt(4*k^2+8*a^2*E))
1001                                       /(2*sqrt(2)*sqrt(-E)*E);
1002
1003 /* Bug 1741705 */
1004 factor(integrate(1/(sin(x)^2+1),x,0,8));
1005 (atan(sqrt(2)*sin(8)/cos(8))+3*%pi)/sqrt(2);
1006
1007 ratsimp(integrate(1/(sin(x)^2+1),x,-8,0));
1008 (atan(sqrt(2)*sin(8)/cos(8))+3*%pi)/sqrt(2);
1009
1010 ratsimp(integrate(1/(sin(x/3)^2+1),x,0,24));
1011 (3*atan(sqrt(2)*sin(8)/cos(8))+9*%pi)/sqrt(2);
1012
1013 ratsimp(integrate(1/(sin(x-3)^2+1),x,3,11));
1014 (atan(sqrt(2)*sin(8)/cos(8))+3*%pi)/sqrt(2);
1015
1016 /* defint log(sqrt(q^2-1)+1) asks about YX - ID: 924868 */
1017 integrate( log(sqrt(q^2-1)+1),q,0,1);
1018 'integrate( log(sqrt(q^2-1)+1),q,0,1);
1019 /*  should be
1020  *  ((4*asinh(1)+(2-sqrt(2))*%i*%pi+(-2^(3/2)))/2^(3/2));
1021  */
1022
1023 /* Bug 1748168 */
1024 integrate(1/(2+cos(x)),x,-%pi/2,%pi/2);
1025 (2/9)*sqrt(3)*%pi;
1026
1027 /* Bug 1781537 */
1028 integrate(cos(x-c),x,c,%pi/2+c);
1029 1;
1030
1031 integrate(cos(x)^3/sin(x)^4,x,1,2);
1032 1/sin(2)-1/(3*sin(2)^3)-1/sin(1)+1/(3*sin(1)^3);
1033
1034 /* [ 1828956 ] Integration yields wrong results */
1035 /* assume & integrate - ID: 3522750 */
1036 (assume(x>0,l>x),0);
1037 0;
1038 integrate(integrate((acos(x/l)+acos(y/l)-%pi/2)/(2*%pi),y,0,sqrt(l^2-x^2)),x,0,l);
1039 l^2/(4*%pi);
1040 is(l>x);
1041 true;
1042 (forget(x>0,l>x),0);
1043 0;
1044
1045 /* principal value integral */
1046 /* [ 657382 ] defint/limit infinite loop */
1047
1048 /* Because of changes to timesin this does no longer work
1049 integrate(1/(1-x^5), x, 0, inf);
1050 0-(2*sqrt(2*sqrt(5)+10)*atan((sqrt(5)-3)*sqrt(2*sqrt(5)+10)/(4*sqrt(5)))
1051         +2*sqrt(10-2*sqrt(5))*atan(sqrt(10-2*sqrt(5))*(sqrt(5)+3)/(4*sqrt(5)))
1052         -sqrt(2)*sqrt(sqrt(5)+5)*%pi-sqrt(2)*sqrt(5-sqrt(5))*%pi)
1053         /20;
1054 */
1055
1056 integrate(1/x, x, minf, inf);
1057 0;
1058
1059 /* #2513 wrong principle value integral */
1060 integrate(1/((x-1)*(x-2)),x,0,3);
1061 -2*log(2);
1062
1063 /* integrate(sin(x)^2/x,x,minf,inf) gives not zero - ID: 3480954 */
1064 integrate(sin(x)^2/x,x,minf,inf);
1065 0;
1066
1067 /* [ 1051437 ] Trig integral error */
1068 integrate(2*cot(x)^2*cos(2*x)/(csc(2*x)+cot(2*x)), x);
1069 log(sin(x)^2+cos(x)^2+2*cos(x)+1) + log(sin(x)^2+cos(x)^2-2*cos(x)+1) + cos(2*x);
1070
1071 /* [ 1960200 ] integrate function raises "too many contexts" error */
1072 float(rectform(integrate(integrate(x^2*y^2*(sin(2.*x)+%i*cos(2.*y))*exp(%i*20.*x)*exp(%i*30.*y),x,-1,-10),y,-5,15)));
1073 6.127592006717518 - 40.62195392101798 * %i;
1074
1075 /* [ 2210087 ] integrate((x+1)^2.0,x) loops endlessly */
1076 integrate((x+1)^2.0,x);
1077 .3333333333333333*(x+1)^3.0;
1078
1079 integrate(x^3/(1 + + 4*x^2 + 6*x^4 + 4*x^6 + x^8 ), x, 0, inf);
1080 1/12;
1081
1082 /* [ 1668087 ] integrate(1/cosh(x/R)^4,x,-inf,inf) = 0 */
1083 /* already assumed a>0 */
1084 integrate(1/cosh(x/a)^4,x,-inf,inf);
1085 4*a/3;
1086
1087 /* [ 1309432 ] integrate(1/cosh(a*x)^2,x,-inf,inf); */
1088 /* already assumed a>0 */
1089 integrate(1/cosh(x/a)^2,x,-inf,inf);
1090 2*a;
1091
1092 /* [ 1860487 ] MONSTERTRIG endless recursion */
1093 integrate (x*(cos(2*x) + sin(x)), x);
1094 (2*x*sin(2*x)+cos(2*x))/4+sin(x)-x*cos(x);
1095
1096 /* [ 1631094 ] integrate(sin(n*x)*x, x) => linear time when n is an integer */
1097 integrate(sin(1000000000*x)*x, x);
1098 (sin(1000000000*x)-1000000000*x*cos(1000000000*x))/1000000000000000000;
1099
1100 /* [ 1847543 ] Integration problem of a special periodic function */
1101 integrate( 1/(11/10+sin(2*%pi*x)), x);
1102 10*atan((22*sin(2*%pi*x)/(cos(2*%pi*x)+1)+20)/(2*sqrt(21)))
1103         /(sqrt(21)*%pi);
1104
1105 /* [ 1054472 ] defint(log(1+exp(A+B*cos(phi))),phi,0,%pi) wrong */
1106
1107 /* Commenting this example out
1108  * There is a problem: The integral does not return the noun form, but
1109  * gives an error: "sign: argument cannot be imaginary; found %i".
1110  * If we execute this integral within the testsuite,
1111  * this integral seems to loop endlessly
1112
1113 integrate(log(1+exp(cos(phi))),phi,0,%pi);
1114 'integrate(log(1+exp(cos(phi))),phi,0,%pi);
1115
1116 */
1117
1118 /* atan2 function is a special case in lisp function INTEGRATOR.
1119    There were no regression tests for it. */
1120 integrate(atan2(y,x),x);
1121 y*log(y^2+x^2)/2+x*atan(y/x);
1122
1123 /* Note: Fixing bug #3246 changed the answer to an equivalent one.
1124    Therefore, compare the ratsimp-ed difference of the expressions to zero. */
1125 ratsimp(integrate(atan2(y,x),y) - (y*atan(y/x)-(x*log(y^2/x^2+1)/2)));
1126 0;
1127
1128 /* Bug ID: 3023997 - integrate(x*atan2(x,y),x) -> noun form
1129  * We have generalized the special case for the atan2 function.
1130  */
1131 integrate(y*atan2(y,x),y);
1132 y^2*atan(y/x)/2-(x^2*y-x^3*atan(y/x))/(2*x);
1133
1134 /* [ 2501765 ] integrate((-14*x^2-32)/(x^4+3*x^2+1)^2,x,0,inf); */
1135 integrate(1/(x^4+3*x^2+1)^2,x,0,inf);
1136 (sqrt(3-sqrt(5))*(3*sqrt(2)*sqrt(5)+15*sqrt(2))*%pi
1137   +sqrt(sqrt(5)+3)*(15*sqrt(2)-3*sqrt(2)*sqrt(5))*%pi)
1138 /400;
1139
1140 /* integrate(exp(-x^n),x,0,1) bug for n >2 - ID: 3469184 */
1141 integrate(exp(-t^4),t,0,1);
1142 (gamma(1/4)-gamma_incomplete(1/4,1))/4;
1143
1144 integrate(exp(sqrt(x)),x,1,5);
1145 2*(sqrt(5)*%e^sqrt(5)-%e^sqrt(5));
1146
1147  /* [ 1899352 ] integrate asks about (y-1)(y+1) after assume(y^2>1) */
1148 (assume(y^2>1),0);
1149 0;
1150 ratsimp(integrate(log(x^2+y^2),x,0,1));
1151 log(y^2+1)+2*atan(1/y)*y-2;
1152 (forget(y^2>1),0);
1153 0;
1154
1155 /* integrate(sqrt(t)*log(t)^(1/2),t,0,1) wrong sign - ID: 2847436 */
1156 integrate(sqrt(t)*log(t)^(1/2),t,0,1);
1157 sqrt(2)*sqrt(%pi)*%i/(3^(3/2));
1158
1159 /* #2732 wrong answer for similar to gaussian integral */
1160 integrate(exp(-x^2-1/x^2),x,-inf,inf);
1161 %e^-2*sqrt(%pi);
1162
1163 /* Integrals of the type: c*z^w*log(z)^m*exp(-t^s) for 0 to inf */
1164
1165 integrate(exp(-t),t,0,inf);
1166 1;
1167 integrate(exp(-t^2),t,0,inf);
1168 sqrt(%pi)/2;
1169 integrate(exp(-t^3),t,0,inf);
1170 gamma(1/3)/3;
1171 integrate(exp(-t^4),t,0,inf);
1172 gamma(1/4)/4;
1173
1174 integrate(t*exp(-t),t,0,inf);
1175 1;
1176 integrate(t^2*exp(-t),t,0,inf);
1177 2;
1178 integrate(t^3*exp(-t),t,0,inf);
1179 6;
1180 integrate(t^4*exp(-t),t,0,inf);
1181 24;
1182
1183 integrate(t*exp(-t^2),t,0,inf);
1184 1/2;
1185 integrate(t^2*exp(-t^2),t,0,inf);
1186 sqrt(%pi)/4;
1187 integrate(t^3*exp(-t^2),t,0,inf);
1188 1/2;
1189 integrate(t^4*exp(-t^2),t,0,inf);
1190 3*sqrt(%pi)/8;
1191
1192 integrate(t^-2*exp(-t^-2),t,0,inf);
1193 sqrt(%pi)/2;
1194 integrate(t^-3*exp(-t^-2),t,0,inf);
1195 1/2;
1196 integrate(t^-4*exp(-t^-2),t,0,inf);
1197 sqrt(%pi)/4;
1198
1199 integrate(t^(1/2)*exp(-t^2),t,0,inf);
1200 gamma(3/4)/2;
1201 integrate(t^(1/3)*exp(-t^2),t,0,inf);
1202 gamma(2/3)/2;
1203 integrate(t^(2/3)*exp(-t^2),t,0,inf);
1204 gamma(5/6)/2;
1205
1206 integrate(exp(-t)*log(t),t,0,inf);
1207 -%gamma;
1208 integrate(t*exp(-t)*log(t),t,0,inf);
1209 1-%gamma;
1210 integrate(t^2*exp(-t)*log(t),t,0,inf);
1211 3-2*%gamma;
1212
1213 integrate(t*exp(-t)*log(t)^2,t,0,inf);
1214 (%pi^2+6*%gamma^2-12*%gamma)/6;
1215 integrate(t*exp(-t^(1/2))*log(t)^2,t,0,inf);
1216 8*%pi^2+48*%gamma^2-176*%gamma+96;
1217 integrate(t*exp(-t^2)*log(t),t,0,inf);
1218 -%gamma/4;
1219
1220
1221 /* Integrals of the type: c*z^w*log(z)^m*exp(-t^s) for x to inf */
1222
1223 (assume(x>0),done);
1224 done;
1225
1226 /* Expand Gamma and Incomplete Gamma functions */
1227 gamma_expand:true;
1228 true;
1229
1230 integrate(exp(-t),t,x,inf);
1231 %e^-x;
1232 integrate(exp(-t^2),t,x,inf);
1233 sqrt(%pi)*erfc(x)/2;
1234 integrate(exp(-t^3),t,x,inf);
1235 gamma_incomplete(1/3,x^3)/3;
1236 integrate(exp(-t^4),t,x,inf);
1237 gamma_incomplete(1/4,x^4)/4;
1238
1239 integrate(t*exp(-t),t,x,inf);
1240 (x+1)*%e^-x;
1241 integrate(t^2*exp(-t),t,x,inf);
1242 (x^2+2*x+2)*%e^-x;
1243 integrate(t^3*exp(-t),t,x,inf);
1244 (x^3+3*x^2+6*x+6)*%e^-x;
1245 integrate(t^4*exp(-t),t,x,inf);
1246 (x^4+4*x^3+12*x^2+24*x+24)*%e^-x;
1247
1248 integrate(t*exp(-t^2),t,x,inf);
1249 %e^-x^2/2;
1250 integrate(t^2*exp(-t^2),t,x,inf);
1251 %e^-x^2*(sqrt(%pi)*%e^x^2*erfc(x)+2*x)/4;
1252 integrate(t^3*exp(-t^2),t,x,inf);
1253 (x^2+1)*%e^-x^2/2;
1254 integrate(t^4*exp(-t^2),t,x,inf);
1255 %e^-x^2*(3*sqrt(%pi)*%e^x^2*erfc(x)+4*x^3+6*x)/8;
1256
1257 integrate(t^-2*exp(-t^-2),t,x,inf);
1258 -sqrt(%pi)*(erfc(1/x)-1)/2;
1259 integrate(t^-3*exp(-t^-2),t,x,inf);
1260 %e^-(1/x^2)*(%e^(1/x^2)-1)/2;
1261 integrate(t^-4*exp(-t^-2),t,x,inf);
1262 -%e^-(1/x^2)*(sqrt(%pi)*(erfc(1/x)-1)*x*%e^(1/x^2)+2)/(4*x);
1263
1264 integrate(t^(1/2)*exp(-t^2),t,x,inf);
1265 gamma_incomplete(3/4,x^2)/2;
1266 integrate(t^(1/3)*exp(-t^2),t,x,inf);
1267 gamma_incomplete(2/3,x^2)/2;
1268 integrate(t^(2/3)*exp(-t^2),t,x,inf);
1269 gamma_incomplete(5/6,x^2)/2;
1270
1271 integrate(exp(-t)*log(t),t,x,inf);
1272 -%e^-x*((%e^x-1)*log(x)+(%gamma-hypergeometric_regularized([1,1],[2,2],-x)*x)
1273                         *%e^x);
1274 integrate(t*exp(-t)*log(t),t,x,inf);
1275 -%e^-x*((%e^x-x-1)*log(x)+(-hypergeometric_regularized([2,2],[3,3],-x)*x^2
1276                           +%gamma-1)
1277                           *%e^x);
1278 integrate(t^2*exp(-t)*log(t),t,x,inf);
1279 -%e^-x*((2*%e^x-x^2-2*x-2)*log(x)+(-4*hypergeometric_regularized(
1280                                       [3,3],[4,4],-x)*x^3
1281                                   +2*%gamma-3)
1282                                   *%e^x);
1283
1284 integrate(t*exp(-t)*log(t),t,x,inf);
1285 -%e^-x*((%e^x-x-1)*log(x)+(-hypergeometric_regularized([2,2],[3,3],-x)*x^2
1286                           +%gamma-1)
1287                           *%e^x);
1288 integrate(t*exp(-t^(1/2))*log(t),t,x,inf);
1289 %e^-sqrt(x)*((-12*%e^sqrt(x)+6*x+12)*log(x)+sqrt(x)*(2*x+12)*log(x)
1290                                            +(144
1291                                             *hypergeometric_regularized(
1292                                              [4,4],[5,5],-sqrt(x))*x^2
1293                                             -24*%gamma+44)
1294                                             *%e^sqrt(x));
1295
1296 /* Integrals of the type: c*z^r*log(z)^n*(1-z)^s*log(1-z)^m */
1297
1298 integrate(t*log(t)*(1-t)*log(1-t),t,0,1);
1299 (3*%pi^2-37)/-108;
1300 integrate(t^2*log(t)*(1-t)*log(1-t),t,0,1);
1301 (3*%pi^2-37)/-216;
1302 integrate(t^2*log(t)^2*(1-t)*log(1-t),t,0,1);
1303 (576*zeta(3)+56*%pi^2-1317)/3456;
1304 integrate(t^2*log(t)^2*(1-t)^2*log(1-t),t,0,1);
1305 (72000*zeta(3)+9400*%pi^2-191243)/1080000;
1306 integrate(t^2*log(t)^2*(1-t)^2*log(1-t)^2,t,0,1);
1307 (3384000*zeta(3)+6000*%pi^4+747800*%pi^2-12135541)/-16200000;
1308
1309 /* integrate -> too many contexts - ID: 2838268 */
1310 integrate(t*cos(1.0*t),t);
1311 t*sin(t)+cos(t);
1312
1313 /* definite integral - bad answer - ID: 2880886 */
1314 integrate(cos(x)^2 * (1 + sin(x)^2)^-3,x,0,%pi/2);
1315 7*%pi/2^(11/2);
1316
1317 (forget(x>0),gamma_expand:false,done);
1318 done;
1319
1320 /* BUG ID: 2906049 - integration failure with option integrate_use_rootsof
1321  * This example is more simple than the example in the bug report.
1322  */
1323 %rnum:0;
1324 0;
1325 integrate((a*x^2+b*x+c)/(d*x^3+a*x^2+b),x),integrate_use_rootsof:true;
1326 'lsum((c+%r1*b+%r1^2*a)*log(x-%r1)/(3*%r1^2*d+2*%r1*a),%r1,
1327             rootsof(d*%r1^3+a*%r1^2+b,%r1));
1328
1329 /* fourier integral incorrect - ID: 2875784 */
1330 integrate(cos(x)*sin(x)/x, x, minf, inf);
1331 %pi/2;
1332
1333 integrate(exp(%i*x)*sin(x)/x, x, minf, inf);
1334 %pi/2;
1335
1336 integrate(exp(-%i*x)*sin(x)/x, x, minf, inf);
1337 %pi/2;
1338
1339 /* Bug ID: 655270 - Incomplete integration
1340  */
1341 integrate(sin(3*asin(x)),x);
1342 (4*x^4-6*x^2+1)/-4;
1343 integrate(sin(4*asin(x)),x);
1344 (15*((1-x^2)^(3/2)/5-x^2*(1-x^2)^(3/2)/5)
1345        +sqrt(1-x)*sqrt(x+1)*(93*x^4-106*x^2+13))/-60;
1346
1347 /* Bug ID: 3029610 - integrate and %e_to_numlog
1348 */
1349 integrate((exp(3/2*log(x))+1)/sqrt(x),x);
1350 2*(x^2/4+sqrt(x));
1351
1352 integrate(3*sqrt(x)*exp(3/2*log(x)),x);
1353 x^3;
1354
1355 /* Someday this should return 2*%pi. And it does now. */
1356 block([ans : integrate(exp(cos(x))*cos(sin(x)),x,0,2*%pi)],[ans, ratsimp(ans)]);
1357 [%i*(-expintegral_ei(1)-%i*%pi)-%i*(%i*%pi-expintegral_ei(1)),2*%pi]$
1358
1359 /* Bug ID: 3039452 - integrate(sqrt(t^c)/(t*(b*t^c+a)),t) hangs
1360  */
1361 /* These assumptions are already present. Consider to delete it more early.
1362 assume(a>0,b>0);
1363 [a>0,b>0]; */
1364 integrate(sqrt(t^c)/(t*(b*t^c+a)),t);
1365 2*atan(sqrt(b)*%e^(c*log(t)/2)/sqrt(a))/(sqrt(a)*sqrt(b)*c);
1366 forget(a>0,b>0);
1367 [a>0,b>0];
1368
1369 /* Bug ID: 3045559 integrate(exp(-u^2), u, minf, x) => incorrect gamma_incomple
1370  */
1371 integrate(exp(-u^2),u,minf,x);
1372 sqrt(%pi)*erf(x)/2+sqrt(%pi)/2;
1373
1374 /* integration error - ID: 3085498 */
1375 integrate(sqrt(1-cos(t)),t,0,2*%pi);
1376 2^(5/2);
1377
1378 /* incorrect integration - ID: 3153533 */
1379 assume(x>2);
1380 [x>2];
1381 integrate(1/(x-z),z,1,x-1);
1382 log(x-1);
1383 forget(x>2);
1384 [x>2];
1385
1386 /* integrate(x^2*exp(x)/(1+exp(x))^2, x, minf, inf); - ID: 3158526 */
1387 integrate(x^2*exp(x)/(1+exp(x))^2, x, minf, inf);
1388 %pi^2/3;
1389
1390 /* integrate(abs(sin(x)),x,0,2*%pi) wrong result - ID: 3165872 */
1391 integrate(abs(sin(x)),x,0,2*%pi);
1392 'integrate(abs(sin(x)),x,0,2*%pi);
1393
1394 /* integrate(sin(2x)atan(sin(x)),x) - ID: 3199708 */
1395 integrate(sin(2*x)*atan(sin(x)),x,0,%pi/2);
1396 %pi/2 - 1;
1397
1398 /* wrong integration answer - ID: 2989983 */
1399 integrate(cos(w+T)/(1+1/2*cos(T))^2,T,0,2*%pi);
1400 -8*%pi*cos(w)/3^(3/2);
1401
1402 /* wrong sign for integral of e^(-1/x^2) - ID: 3211937 */
1403 integrate(%e^(-1/x^2),x,0,1);
1404 sqrt(%pi)*erf(1)-sqrt(%pi)+%e^(-1)$
1405
1406 /* definite integration introduces xy variable - ID: 3292707 */
1407 integrate((log((2-x)/2)+log(2))/(1-x), x, 0, 1);
1408 li[2](2)-(6*log(2)^2-6*log(-2)*log(2)+%pi^2)/6;
1409
1410 /* Incorrect integral of log(1+a/(x*t)^2) - ID: 3291160 */
1411 integrate(log(1+7/(x^2)),x,1,inf);
1412 sqrt(7)*(2*atan(sqrt(7))-log(8)/sqrt(7))$
1413
1414 /* Wrong sign in exponential integral - ID: 3517785 */
1415 integrate(1/(%e^(2*t)*sqrt(1-1/%e^(2*t))), t, 0, inf);
1416 1;
1417
1418 /* integral from minf to inf of an absolute value combo fails - ID: 3313531 */
1419 integrate(abs(x - 1) + abs(x + 1) - 2*abs(x),x,minf,inf);
1420 2;
1421
1422 /* definite integration - ID: 3315476 */
1423 integrate(1/2^x, x, 0, 1);
1424 1/(2*log(2));
1425
1426 /* discontinuities in integral */
1427 /* error in integrating exp(-x)*sinh(sqrt(x)) - ID: 3292033 */
1428 integrate(exp(sqrt(x)-x), x, 0, inf);
1429 %e^(1/4)*gamma_incomplete(1,1/4)-%e^(1/4)*gamma_incomplete(1/2,1/4)/2
1430                                     +%e^(1/4)*sqrt(%pi);
1431
1432 /* integrate(log(t)*log(t+1),t,0,1) gives Lisp error - ID: 3381012 */
1433 integrate(log(t)*log(t+1),t,0,1);
1434 -2*log(2)-%pi^2/12+2;
1435
1436 /* integrate(log(2*sin(t/2)),t,0,%pi) contains false - ID: 3381037 */
1437 integrate(log(2*sin(t/2)),t,0,%pi);
1438 2*(%i*li[2](%i)+%i*li[2](-%i)+(%i*%pi^2)/24)$
1439 /* above result is really zero */
1440
1441 /* integrate error - ID: 3388801 */
1442 integrate( sqrt((x-474)^2/(107669-(x-474)^2)+1), x, 181, 474);
1443 sqrt(67)*sqrt(1607)*%i*log(2*sqrt(107669)*%i)-sqrt(67)*sqrt(1607)*%i*log(2*(2*sqrt(5455)*%i-293))$
1444
1445 /* integrate erf fails - ID: 3454370 */
1446 integrate( erf(x+a)-erf(x-a), x, minf, inf);
1447 4*a;
1448
1449 /* Wrong result for definite integral - ID: 3538167 */
1450 integrate(cos(3*x)/(5-4*cos(x)),x,0,2*%pi);
1451 %pi/12;
1452
1453 /* #2575 Integration error: integrate(sqrt(k-k*cos(2*x)), x) */
1454 integrate(sqrt(cos(x)+1), x, -%pi, %pi);
1455 2^(5/2);
1456
1457 /* #2776 Error when integrate sqrt */
1458 integrate ((1+exp(t))/sqrt(t+exp(t)), t, 0, 1);
1459 2*sqrt(%e+1)-2;
1460
1461 /* #2862 Incorrect result for integrate(u/(u+1)^2,u,0,inf); */
1462 errcatch(integrate(u/(u+1)^2,u,0,inf));
1463 [];
1464
1465 error;
1466 ["defint: integral is divergent."];
1467
1468
1469 /* #2903 Wrong integration over 1/(x+1) with bounds 0 to inf, since version 5.34.0 */
1470 errcatch(integrate(1/(x+1), x,0,inf));
1471 [];
1472
1473 error;
1474 ["defint: integral is divergent."];
1475
1476
1477 /* #2919 Definite integration broken: integrate(1/(x^2), x, -inf, inf) gives zero  */
1478 errcatch(integrate(1/(x^2), x, -inf, inf));
1479 [];
1480
1481 error;
1482 ["defint: integral is divergent."];
1483
1484 /* other tests mentioned in #2919 */
1485
1486 errcatch (integrate(1/x, x, 0, 1));
1487 [];
1488
1489 error;
1490 ["defint: integral is divergent."];
1491
1492 integrate(1/x, x, -1, 1);
1493 0;
1494 /* expected output: Principal Value -- dunno how to capture that */
1495
1496 errcatch (integrate(1/x, x, 0, inf));
1497 [];
1498
1499 error;
1500 ["defint: integral is divergent."];
1501
1502 errcatch (integrate(1/x, x, 1, inf));
1503 [];
1504
1505 error;
1506 ["defint: integral is divergent."];
1507
1508 integrate(1/x, x, -inf, inf);
1509 0;
1510 /* expected output: Principal Value -- dunno how to capture that */
1511
1512 errcatch (integrate(1/x^2, x, 0, 1));
1513 [];
1514
1515 error;
1516 ["defint: integral is divergent."];
1517
1518 errcatch (integrate(1/x^2, x, -1, 1));
1519 [];
1520
1521 error;
1522 ["defint: integral is divergent."];
1523
1524 errcatch (integrate(1/x^2, x, 0, inf));
1525 [];
1526
1527 error;
1528 ["defint: integral is divergent."];
1529
1530 /* #2987 Some divergent integrals give error, some don't */
1531 errcatch (integrate(exp(x)/x, x, 0, 1));
1532 [];
1533
1534 error;
1535 ["defint: integral is divergent."];
1536
1537 integrate(1/x^2, x, 1, inf);
1538 1;
1539
1540 /* SF # 2633: ev(integrate,numer) gives strange result */
1541 (foo : ev(integrate(sin(2*%pi*x),x,0,1),numer), 0);
1542 0;
1543
1544 freeof (false, foo);
1545 true;
1546
1547 float (foo);
1548 0.0;
1549
1550 /* Maxima used to call these integrals divergent even though it would yield 0
1551  * when both limits were inf.
1552  */
1553 integrate(exp(-x^2), x, minf, minf);
1554 0;
1555
1556 integrate(const, x, minf, minf);
1557 0;
1558
1559 /* Maxima used to give the bogus result "-false" for this integral even though
1560  * it would yield 0 when both limits were inf.
1561  */
1562 integrate(erf(x), x, minf, minf);
1563 0;
1564
1565 /* From the mailing list on 29 Aug 2015. Maxima was leaving integrate(0,x) in the result.
1566  */
1567 (reset(integration_constant_counter), integrate(x=0,x));
1568 x^2/2 = %c1;
1569
1570 /* SourceForge bug #3017 - "integrate_use_rootsof" leads to wrong result */
1571 %rnum:0;
1572 0;
1573 integrate(1/(x^(1/3)+x+1),x),integrate_use_rootsof:true;
1574 3*'lsum((%r1^2*log(x^(1/3)-%r1))/(3*%r1^2+1),%r1,
1575               rootsof(%r1^3+%r1+1,%r1));
1576
1577 /* Bug #2507: Strange non-evaluation of integral
1578  *
1579  * The problem was that the second integral was not evaluating
1580  * while the first one was.  The assume is also important because
1581  * the integral would evaluate if answered interactively.
1582  */
1583 (assume(n>0),0);
1584 0;
1585
1586 integrate(1/sqrt(1+x^2*n),x,1,2);
1587 asinh(2*sqrt(n))/sqrt(n)-asinh(sqrt(n))/sqrt(n);
1588
1589 integrate(1/sqrt(1+x^2/n),x,1,2);
1590 asinh(2/sqrt(n))*sqrt(n)-asinh(1/sqrt(n))*sqrt(n);
1591
1592 (forget(n>0),0);
1593 0;
1594
1595 /* mailing list 2015-11-23: "Change in behavior for taylor/integrate - a bug?" */
1596
1597 (kill(f, x, c, deltax), trunc(integrate(taylor(f(x), x, c, 6), x, c-deltax, c+deltax)));
1598 (deltax^7*('at('diff(f(x),x,6),x = c))+42*deltax^5*('at('diff(f(x),x,4),x = c))
1599                                       +840*deltax^3*('at('diff(f(x),x,2),x = c))+5040*f(c)*deltax)
1600                                        /2520$
1601
1602 integrate (taylor (f(x), x, 0, 1), x, a, b);
1603 ((b^2-a^2)*('at('diff(f(x),x,1),x = 0))+2*f(0)*b-2*f(0)*a)/2$
1604
1605 /* mailing list 2015-11-08: "un-nerving behavior of numer... yikes!" */
1606
1607 /* disable this test pending resolution of SF bug #3097 !!
1608 integrate(x*sin(%pi/2*x), x,0,1 ),    numer;
1609 0.4052847345693511;
1610  */
1611
1612 integrate(x*sin(1.5*x), x,0,1 ),    numer;
1613 0.3961729707122222;
1614
1615 integrate(x*cos(1.5*x), x,0,1 ),    numer;
1616 0.2519909695883491;
1617
1618 integrate(x*(cos(2.5*x) + 1.25), x, 0, 1), numer;
1619 0.5762058791540733;
1620
1621 integrate(x*(cos(2.5*x^2 - 1) + 1.25), x, 0, 1), numer;
1622 0.9927931942823904;
1623
1624 integrate (sin(4*asin(x)),x, 0.5, 0.75), numer;
1625 0.09667085773175951;
1626
1627 /* disable these next few tests, since these fail (with a floating point error)
1628  * in the absence of float-to-bigfloat promotion, which was recently reverted.
1629  *
1630 /* These next 4 examples yield big messy rational expressions.
1631  * It's maybe not so great to have to squash them (with bfloat and rectform)
1632  * but at least we get a result; these examples cause a floating-point error
1633  * in previous versions.
1634  */
1635 rectform (bfloat (ev (integrate ((x^2*(%e^(2*%i*x)+%e^-(2*%i*x))*%e^(30*%i*x-200*%i))/2,x, 1.25, 1.5), numer)));
1636 (- 6.358161510436547b-2*%i) - 2.859595163448958b-2;
1637
1638 rectform (bfloat (ev (integrate ((x^2*(%e^(2*%i*x)+%e^-(2*%i*x))*%e^(30*%i*x-20*%i))/2,x, 1.25, 1.5), numer)));
1639 6.096077978946073b-2*%i - 3.382504333513841b-2;
1640
1641 rectform (bfloat (ev (integrate (x^2*%e^(30*%i*x-200*%i)*cos(2*x),x, 1.25, 1.5), numer)));
1642 (- 6.35816151043802b-2*%i) - 2.859595163448362b-2;
1643
1644 rectform (bfloat (ev (integrate (x^2*%e^(30*%i*x-20*%i)*cos(2*x),x, 1.25, 1.5), numer)));
1645 6.096077978946119b-2*%i - 3.382504333513988b-2;
1646  *
1647  */
1648
1649
1650 /* bug #3114 */
1651 integrate(sec(x)/sqrt(tan(x)),x);
1652 'integrate(sec(x)/sqrt(tan(x)),x);
1653
1654 /* bug #3115 */
1655 integrate((x*sqrt(x^3-4*x))/(x^2-4),x);
1656 'integrate((x*sqrt(x^3-4*x))/(x^2-4),x);
1657
1658 /* SF bug #3144: "stackoverflow in integral" */
1659
1660 integrate(sin(k*x)/x*erf(x^2),x,0,inf);
1661 'integrate(sin(k*x)/x*erf(x^2),x,0,inf);
1662
1663 /* Borwein integrals; see: https://en.wikipedia.org/wiki/Borwein_integral */
1664
1665 integrate(sin(x)/x,x,0,inf)$
1666 %pi/2 $
1667
1668 integrate((sin(x)/x) * (sin(x/3)/(x/3)),x,0,inf)$
1669 %pi/2 $
1670
1671 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)),x,0,inf)$
1672 %pi/2 $
1673
1674 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)) * (sin(x/7)/(x/7)),x,0,inf)$
1675 %pi/2 $
1676
1677 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)) * (sin(x/7)/(x/7)) * (sin(x/9)/(x/9)),x,0,inf)$
1678 %pi/2 $
1679
1680 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)) * (sin(x/7)/(x/7)) * (sin(x/9)/(x/9))
1681                      * (sin(x/11)/(x/11)),x,0,inf)$
1682 %pi/2 $
1683
1684 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)) * (sin(x/7)/(x/7)) * (sin(x/9)/(x/9))
1685                      * (sin(x/11)/(x/11)) * (sin(x/13)/(x/13)),x,0,inf)$
1686 %pi/2 $
1687
1688 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)) * (sin(x/7)/(x/7)) * (sin(x/9)/(x/9))
1689                      * (sin(x/11)/(x/11)) * (sin(x/13)/(x/13)) * (sin(x/15)/(x/15)),x,0,inf)$
1690 467807924713440738696537864469*%pi/935615849440640907310521750000 $
1691
1692 /* Bug #2314: Two different results for an integral */
1693
1694 (assume (0 < x, x < 1, 0 < y, y < 1), 0)$
1695 0$
1696
1697 integrate (integrate (sqrt (4*x^2 - y^2), x, y, 1), y, 0, 1)$
1698 (2*sqrt(3)*%pi+9)/(2*3^(5/2))$
1699
1700 integrate (integrate (sqrt (4*x^2 - y^2), y, 0, x), x, 0, 1)$
1701 (2*sqrt(3)*%pi+9)/(2*3^(5/2))$
1702
1703 (forget (0 < x, x < 1, 0 < y, y < 1), 0)$
1704 0$
1705
1706 /* Bug #2295: Failure to evaluate definite integral */
1707
1708 integrate (1 / (sqrt (x) * (1 + sqrt (x))^2), x, 1, 4)$
1709 1/3$
1710
1711 /* SF bug #3274: "quotient by zero" error in integrate */
1712
1713 (expr_tau: (%e^-(t/tau3)*(%e^(t/tau3+t/tau1)-%e^t*tau3))
1714       /((tau1*%e^(t/tau1)+%e^(t/tau1))
1715        *%e^(t/tau2)*tau3+%e^(t/tau2+t/tau1)),
1716  I_tau : integrate(expr_tau, t));
1717 -(%e^((-t/tau2)-t/tau1)*(tau1*tau2*tau3^2*%e^(t-t/tau3)
1718                         +%e^(t/tau1)
1719                          *(tau1*(tau2^2*(tau3-1)-tau2*tau3)
1720                           -tau2^2*tau3)))
1721  /(tau1^2*(tau2*(tau3^2-tau3)-tau3^2)
1722   +tau1*((-tau3^2)-tau3-tau2)+tau2*((-tau3^2)-tau3))$
1723
1724 is(equal(diff(I_tau, t), expr_tau));
1725 true;
1726
1727 (expr_c : subst([tau1=c1, tau2=c2, tau3=c3], expr_tau),
1728  I_c : integrate (expr_c, t));
1729 (-(c1*c2*c3^2*%e^((-t/c3)-t/c2-t/c1+t))
1730  /(c1^2*(c2*(c3^2-c3)-c3^2)+c1*((-c3^2)-c3-c2)+c2*((-c3^2)-c3)))
1731  -(c2*c3*%e^-(t/c2))/(c1*c3^2+c3^2+c3)$
1732
1733 is(equal(diff(I_c, t), expr_c));
1734 true;
1735
1736 constantp (ratsimp (subst([c1=tau1, c2=tau2, c3=tau3], I_c) - I_tau));
1737 true;
1738
1739 /* following was incorrect (returned 0) at some point in recent past
1740  * ensure that it stays correct
1741  * from sage-support 2017-02-22
1742  */
1743
1744 errcatch (integrate(x/(1+x^2),x,0,inf));
1745 [];
1746
1747 /* Bug #3368: integrate('limit(...),...) internal error */
1748
1749 /* Original test cases from bug report: */
1750
1751 integrate('limit(x,x,0),x);
1752 x*'limit(x,x,0);
1753
1754 integrate('limit(x,x,y),x);
1755 x*'limit(x,x,y);
1756
1757 integrate('limit(x,x,x),x);
1758 'integrate('limit(x,x,x),x);
1759
1760 /* Some more: */
1761
1762 integrate('at(x,x=1), x);
1763 x*'at(x,x=1);
1764
1765 integrate('sum(x,x,1,1), x);
1766 x*'sum(x,x,1,1);
1767
1768 integrate('product(x,x,1,1), x);
1769 x*'product(x,x,1,1);
1770
1771 /* SF bug #3712: "Symbolic integration may fail when called with numerical constants in function and/or limits." */
1772
1773 /* original problem from bug report */
1774
1775 (E1:200,
1776  nu:0.33,
1777  rho1:7850,
1778  Pi:ev(%pi,numer),
1779  ah:5,
1780  e0:0.25,
1781  h:1/ah,
1782  em:1-sqrt(1-e0),
1783  aa:Pi*z/(2*h)+Pi/(4),
1784  Ez:E1*(1-e0*(cos(aa))),
1785  rhoz:rho1*(1-em*(cos(aa))),
1786  Q11:Ez,
1787  /* this call to integrate is okay in bug report */
1788  A11:ev(integrate(Q11,z,-h/2,h/2),numer));
1789 33.63380227632413;
1790
1791 /* this one triggers the bug in the bug report */
1792 B11:ev(integrate(Q11*z,z,-h/2,h/2),numer);
1793 0.1739496967711204;
1794
1795 /* double check numerical results */
1796
1797 is (abs (first (quad_qags (Q11,z,-h/2,h/2, 'epsabs = 1e-8)) - 33.63380227632413) < 1e-8);
1798 true;
1799
1800 /* for display_all = true, need nondefault value of linel. */
1801 block([linel : 107], is (abs (first (quad_qags (Q11*z,z,-h/2,h/2, 'epsabs = 1e-8)) - 0.1739496967711204) < 1e-8));
1802 true;
1803
1804 /* additional cases for #3712 */
1805
1806 /* try integral w/ symbolic constants;
1807  * some gyrations here to verify symbolic result
1808  */
1809 (kill(values),
1810  Q:a*(b-c*cos(d*z+e)),
1811  IQz_indefinite: integrate(Q*z, z),
1812  ratsimp (diff (IQz_indefinite, z) - Q*z));
1813 0;
1814
1815 (IQz_defint: integrate(Q*z,z,-f,f),
1816  IQz_defint_via_indefinite: ev (IQz_indefinite, z = f) - ev (IQz_indefinite, z = -f),
1817  ratsimp (IQz_defint - IQz_defint_via_indefinite));
1818 0;
1819
1820 (L1:[a=200, b=1, c=0.25, d=7.853981633974483, e=0.7853981633974483,f=1/10],
1821  subst (L1, IQz_defint));
1822 0.1739496967711208;
1823
1824 /* verify integrate doesn't take a long time with more digits */
1825
1826 /*
1827 integrate(z*(1-cos(7.85*z)),z,0,1) => 0 seconds
1828 integrate(z*(1-cos(7.8539*z)),z,0,1) => 3.2 seconds
1829 integrate(z*(1-cos(7.85398*z)),z,0,1) => 77 seconds
1830  */
1831
1832 (time_in_secs (f, e) ::= buildq ([f, e], (timer (f), e, block ([i: timer_info (f)], untimer (f), i[2, 4]/verbify(sec)))),
1833  is (time_in_secs (integrate, integrate(z*(1-cos(7.85*z)),z,0,1)) < 1));
1834 true;
1835
1836 is (time_in_secs (integrate, integrate(z*(1-cos(7.8539*z)),z,0,1)) < 1);
1837 true;
1838
1839 is (time_in_secs (integrate, integrate(z*(1-cos(7.85398*z)),z,0,1)) < 1);
1840 true;
1841
1842 /* factor_max_degree now controls loops in ROOTFAC -- verify result not changed */
1843
1844 block ([factor_max_degree: 1000], integrate(z*(1-cos(7.8*z)),z,0,1));
1845 25/1521-(390*sin(39/5)+50*cos(39/5)-1521)/3042$
1846
1847 block ([factor_max_degree: 10], integrate(z*(1-cos(7.8*z)),z,0,1));
1848 25/1521-(390*sin(39/5)+50*cos(39/5)-1521)/3042$
1849
1850 /*
1851 B11:'integrate(Q11*z,z,-h/2,h/2)$
1852 B11r: scanmap('ratsimp,B11)$
1853 ev(B11r,nouns) => result in 0.16 seconds
1854  */
1855
1856 (B11:'integrate(Q11*z,z,-h/2,h/2),
1857  B11r: scanmap('ratsimp,B11),
1858  is (time_in_secs (integrate, ev(B11r,nouns)) < 1));
1859 true;
1860
1861 (foo: integrate(z*cos(7.853981633974483*z + 0.7853981633974483),z),
1862  bar: subst (a = 7.853981633974483, integrate(z*cos(a*z + a/10), z)),
1863  is (lmax (makelist (ev (foo - bar, numer), z, [0, 1, 2, 3, 4, 5])) < 1e-12));
1864 true;
1865
1866 /* bug reported to mailing list 2021-02-18: "Maxima asking for the sign of 'und'" */
1867
1868 kill (a, b, c, d, u, v, w, x, y, z);
1869 done;
1870
1871 integrate(cos(x),x,0,inf);
1872 'integrate(cos(x),x,0,inf);
1873
1874 integrate(sin(u) + cos(u), u, 0, inf);
1875 'integrate(sin(u) + cos(u), u, 0, inf);
1876
1877 integrate (a*sin(y)+b*cos(y), y, minf, 1);
1878 'integrate (a*sin(y)+b*cos(y), y, minf, 1);
1879
1880 integrate ((a*sin(z)+b*cos(z))^3, z, minf, 1);
1881 'integrate ((a*sin(z)+b*cos(z))^3, z, minf, 1);
1882
1883 integrate (a*sin(5*w)+b*cos(7*w), w, minf, 1);
1884 'integrate (a*sin(5*w)+b*cos(7*w), w, minf, 1);
1885
1886 /* SF bug #3718: "incorrect trigonometric definite integral" */
1887
1888 (kill(e, e1, e2, e3, t, I),
1889  e:(tan(x)^2+1)^2/(16*tan(x)^4+12*tan(x)^2+1),
1890  [e1, e2, e3]: args (expand (e)),
1891  assume (t > acos(3/5)/2),
1892  0);
1893 0;
1894
1895 /* RISCHINT seems to be returning a discontinuous antiderivative
1896  * (a separate problem from #3718).
1897  * First just verify it's not zero (the previous, incorrect result).
1898  */
1899
1900 is ((I: integrate (e, x, 0, t)) # 0);
1901 true;
1902
1903 /* Now try some random intervals and see if the antiderivative
1904  * matches a numerical calculation on at least one such interval.
1905  *
1906  * Trying this repeated yields 14 failures (apparently hit a discontinuity)
1907  * and 39 successess (apparently missed).
1908  * If I am reading the tea leaves correctly, that suggests
1909  * the mean number of tries to the first success is 1/(39/(39 + 14)) = 53/39 = 1.36
1910  * and p(number of tries > 10) = (1 - 39/(39 + 14))^10 = (14/53)^10 = 1.65 x 10^-6 = 1/605000.
1911  * So at least for this problem, 10 tries seems enough.
1912  */
1913
1914 block ([t1, t2, I_defint_numerical, I_defint_FTC],
1915   for k thru 10
1916     do (t1: random(5.0),
1917         t2: t1 + random(0.5),
1918         I_defint_numerical: first (quad_qags (e, x, t1, t2)),
1919         I_defint_FTC: float (ev (I, t = t2) - ev (I, t = t1)),
1920         if abs (I_defint_numerical - I_defint_FTC) < 1e-4
1921           then (printf (true, "re: #3718, apparently valid antiderivative on (~f, ~f)~%", t1, t2),
1922                 return(true))
1923           else printf (true, "re: #3718, can't confirm antiderivative on (~f, ~f), keep trying.~%", t1, t2)));
1924 true;
1925
1926 /* Bug #4018: defint(foo,,0,inf) lisp error when denom(foo) contains %i */
1927
1928 defint(1/((x+%i/2)*(x+1)),x,0,inf);
1929 ((2*%i+4)*log(2)+(1-2*%i)*%pi)/5;
1930
1931 defint(1/((x+%i)*(x/2+1)),x,0,inf);
1932 ((2*%i+4)*log(2)+(1-2*%i)*%pi)/5;
1933
1934 defint(1/((x+%i/2)*(x-%i)),x,minf,inf);
1935 (4*%pi)/3;
1936
1937 block ([g_poles: subst(phi=%pi/3, exp(%i * [phi, -phi, %pi-phi, %pi+phi]))],
1938   defint(1/apply("*", makelist(w - pole, pole, g_poles)), w, 0, inf));
1939 %pi/(2*sqrt(3));
1940
1941 /* SF bug #3810: "integrate error \"not of type FIXNUM\" for integrand with floats in it"
1942  * The following 20 cases are derived from the example shown in the bug report;
1943  * I wasn't able to find a smaller set.
1944  * For the record, I've only observed the bug, in which there's
1945  * a Lisp error about some huge integer which doesn't fit into a fixnum,
1946  * with SBCL.
1947  */
1948 (if not ?fboundp ('read_matrix) then load (numericalio),
1949  construct_float (i) := ?construct\-float\-64\-from\-integer (i),
1950  ratprint:false,
1951  ratfac:true,
1952  algebraic:true,
1953  factor_max_degree_print_warning:false,
1954 0);
1955 0;
1956
1957 (integrate (x^4*((construct_float(13832102834422948280)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(4659109193942506099)*sqrt(x)*%e^(construct_float(13852053825502851604)*x))^2,x,0,inf),
1958 integrate (x^4*((construct_float(4603705080989522416)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(13876663004335033608)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(4680418835321724124)*x^(5/2)*%e^(construct_float(13850392950977358476)*x))^2,x,0,inf),
1959 integrate (x^4*((construct_float(13822985546564103824)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(4650040524666783766)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(13898864330314251671)*x^(5/2)*%e^(construct_float(13850392950977358476)*x)+construct_float(4689176719312114346)*x^(9/2)*%e^(construct_float(13849077382790281914)*x))^2,x,0,inf),
1960 integrate (x^4*((construct_float(13815380123310259942)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(4642544339172547822)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(13892017095959042546)*x^(5/2)*%e^(construct_float(13850392950977358476)*x)+construct_float(4682273596335467143)*x^(9/2)*%e^(construct_float(13849077382790281914)*x)+construct_float(13906083594305576361)*x^(13/2)*%e^(construct_float(13847501804967206881)*x)+construct_float(4676345781606752784)*x^(17/2)*%e^(construct_float(13845850997455375733)*x))^2,x,0,inf),
1961 integrate (x^4*((construct_float(4588287034702523055)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(13862081432535911133)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(4665543650886122250)*x^(5/2)*%e^(construct_float(13850392950977358476)*x)+construct_float(13901162363127812545)*x^(9/2)*%e^(construct_float(13849077382790281914)*x)+construct_float(4680543130755504491)*x^(13/2)*%e^(construct_float(13847501804967206881)*x)+construct_float(13894925586695210887)*x^(17/2)*%e^(construct_float(13845850997455375733)*x)+construct_float(4656829810530318229)*x^(21/2)*%e^(construct_float(13844543403286028560)*x))^2,x,0,inf),
1962 integrate (x^4*((construct_float(13807950446123059420)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(4634715679194593503)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(13884770780167979999)*x^(5/2)*%e^(construct_float(13850392950977358476)*x)+construct_float(4674631855308846949)*x^(9/2)*%e^(construct_float(13849077382790281914)*x)+construct_float(13899442112350166625)*x^(13/2)*%e^(construct_float(13847501804967206881)*x)+construct_float(4669380157924531903)*x^(17/2)*%e^(construct_float(13845850997455375733)*x)+construct_float(13875646849200254557)*x^(21/2)*%e^(construct_float(13844543403286028560)*x)+construct_float(4628334369197177667)*x^(25/2)*%e^(construct_float(13842950077923898657)*x))^2,x,0,inf),
1963 integrate (x^4*((construct_float(4580499919285849890)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(13853982548367093798)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(4657378199720877552)*x^(5/2)*%e^(construct_float(13850392950977358476)*x)+construct_float(13893812467321528715)*x^(9/2)*%e^(construct_float(13849077382790281914)*x)+construct_float(4672405349716730081)*x^(13/2)*%e^(construct_float(13847501804967206881)*x)+construct_float(13888488673177254699)*x^(17/2)*%e^(construct_float(13845850997455375733)*x)+construct_float(4649144676243844169)*x^(21/2)*%e^(construct_float(13844543403286028560)*x)+construct_float(13847692918362499544)*x^(25/2)*%e^(construct_float(13842950077923898657)*x)+construct_float(4591940524810724167)*x^(29/2)*%e^(construct_float(13841309276406812944)*x))^2,x,0,inf),
1964 integrate (x^4*((construct_float(13799730948534155654)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(4626522136374979242)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(13876608944864255440)*x^(5/2)*%e^(construct_float(13850392950977358476)*x)+construct_float(4666332164086245695)*x^(9/2)*%e^(construct_float(13849077382790281914)*x)+construct_float(13891692449046899557)*x^(13/2)*%e^(construct_float(13847501804967206881)*x)+construct_float(4661352201463965459)*x^(17/2)*%e^(construct_float(13845850997455375733)*x)+construct_float(13868358252270079858)*x^(21/2)*%e^(construct_float(13844543403286028560)*x)+construct_float(4621086115000026761)*x^(25/2)*%e^(construct_float(13842950077923898657)*x)+construct_float(13811978465200383126)*x^(29/2)*%e^(construct_float(13841309276406812944)*x)+construct_float(4547901422001664002)*x^(33/2)*%e^(construct_float(13840009607922498019)*x))^2,x,0,inf),
1965 integrate (x^4*((construct_float(4572184908458273386)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(13845758088780195921)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(4649087273240521486)*x^(5/2)*%e^(construct_float(13850392950977358476)*x)+construct_float(13885518154380124499)*x^(9/2)*%e^(construct_float(13849077382790281914)*x)+construct_float(4664285118797405184)*x^(13/2)*%e^(construct_float(13847501804967206881)*x)+construct_float(13880505870680953441)*x^(17/2)*%e^(construct_float(13845850997455375733)*x)+construct_float(4640946536814966454)*x^(21/2)*%e^(construct_float(13844543403286028560)*x)+construct_float(13840303268528715398)*x^(25/2)*%e^(construct_float(13842950077923898657)*x)+construct_float(4584986699114512377)*x^(29/2)*%e^(construct_float(13841309276406812944)*x)+construct_float(13768183034564966614)*x^(33/2)*%e^(construct_float(13840009607922498019)*x)+construct_float(4495887554581308620)*x^(37/2)*%e^(construct_float(13838398642593989805)*x))^2,x,0,inf),
1966 integrate (x^4*((construct_float(13791339657508427891)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(4618203091029411652)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(13868259498832129124)*x^(5/2)*%e^(construct_float(13850392950977358476)*x)+construct_float(4657923725573991758)*x^(9/2)*%e^(construct_float(13849077382790281914)*x)+construct_float(13883544990328347659)*x^(13/2)*%e^(construct_float(13847501804967206881)*x)+construct_float(4652887676647289995)*x^(17/2)*%e^(construct_float(13845850997455375733)*x)+construct_float(13860171559518208449)*x^(21/2)*%e^(construct_float(13844543403286028560)*x)+construct_float(4612726215435183695)*x^(25/2)*%e^(construct_float(13842950077923898657)*x)+construct_float(13804177960580528374)*x^(29/2)*%e^(construct_float(13841309276406812944)*x)+construct_float(4540904542485620193)*x^(33/2)*%e^(construct_float(13840009607922498019)*x)+construct_float(13716623792620785711)*x^(37/2)*%e^(construct_float(13838398642593989805)*x)+construct_float(4436689235760666063)*x^(41/2)*%e^(construct_float(13836767786422585277)*x))^2,x,0,inf),
1967 integrate (x^4*((construct_float(4563694232363278599)*%e^(construct_float(13853611547558191654)*x))/x^(3/2)+construct_float(13837329261576760585)*sqrt(x)*%e^(construct_float(13852053825502851604)*x)+construct_float(4640628279602305787)*x^(5/2)*%e^(construct_float(13850392950977358476)*x)+construct_float(13877017855321991145)*x^(9/2)*%e^(construct_float(13849077382790281914)*x)+construct_float(4655984262259314301)*x^(13/2)*%e^(construct_float(13847501804967206881)*x)+construct_float(13871961325783836989)*x^(17/2)*%e^(construct_float(13845850997455375733)*x)+construct_float(4632579512185702728)*x^(21/2)*%e^(construct_float(13844543403286028560)*x)+construct_float(13831820979204054045)*x^(25/2)*%e^(construct_float(13842950077923898657)*x)+construct_float(4576528486454867609)*x^(29/2)*%e^(construct_float(13841309276406812944)*x)+construct_float(13760093515846203489)*x^(33/2)*%e^(construct_float(13840009607922498019)*x)+construct_float(4489545647912134320)*x^(37/2)*%e^(construct_float(13838398642593989805)*x)+construct_float(13657416847333971758)*x^(41/2)*%e^(construct_float(13836767786422585277)*x)+construct_float(4369908098040789962)*x^(45/2)*%e^(construct_float(13835475995583563131)*x))^2,x,0,inf),
1968 integrate (x^4*(construct_float(13874334652917149491)*sqrt(x)*%e^(construct_float(13851818867533264930)*x)+construct_float(4676446327894293873)*x^(5/2)*%e^(construct_float(13850206842258521600)*x))^2,x,0,inf),
1969 integrate (x^4*(construct_float(4647407053860680595)*sqrt(x)*%e^(construct_float(13851818867533264930)*x)+construct_float(13897605418546479799)*x^(5/2)*%e^(construct_float(13850206842258521600)*x)+construct_float(4685350810584663793)*x^(9/2)*%e^(construct_float(13848929967253417878)*x))^2,x,0,inf),
1970 integrate (x^4*(construct_float(13867011523143065883)*sqrt(x)*%e^(construct_float(13851818867533264930)*x)+construct_float(4670847344307255816)*x^(5/2)*%e^(construct_float(13850206842258521600)*x)+construct_float(13906833880713031021)*x^(9/2)*%e^(construct_float(13848929967253417878)*x)+construct_float(4683948219867094378)*x^(13/2)*%e^(construct_float(13847268271141981213)*x))^2,x,0,inf),
1971 integrate (x^4*(construct_float(4639781555154818310)*sqrt(x)*%e^(construct_float(13851818867533264930)*x)+construct_float(13890485364949001659)*x^(5/2)*%e^(construct_float(13850206842258521600)*x)+construct_float(4680103856627696714)*x^(9/2)*%e^(construct_float(13848929967253417878)*x)+construct_float(13904723712727785251)*x^(13/2)*%e^(construct_float(13847268271141981213)*x)+construct_float(4672513506043682445)*x^(17/2)*%e^(construct_float(13845666016792262747)*x))^2,x,0,inf),
1972 integrate (x^4*(construct_float(13859309570165497800)*sqrt(x)*%e^(construct_float(13851818867533264930)*x)+construct_float(4663344075634158250)*x^(5/2)*%e^(construct_float(13850206842258521600)*x)+construct_float(13899734443221025783)*x^(9/2)*%e^(construct_float(13848929967253417878)*x)+construct_float(4678356998751580204)*x^(13/2)*%e^(construct_float(13847268271141981213)*x)+construct_float(13893686855614096759)*x^(17/2)*%e^(construct_float(13845666016792262747)*x)+construct_float(4652770339619897509)*x^(21/2)*%e^(construct_float(13844396881274955196)*x))^2,x,0,inf),
1973 integrate (x^4*(construct_float(4632067885786510673)*sqrt(x)*%e^(construct_float(13851818867533264930)*x)+construct_float(13882897799064936059)*x^(5/2)*%e^(construct_float(13850206842258521600)*x)+construct_float(4672535462110635372)*x^(9/2)*%e^(construct_float(13848929967253417878)*x)+construct_float(13898168598415961189)*x^(13/2)*%e^(construct_float(13847268271141981213)*x)+construct_float(4666779889548038673)*x^(17/2)*%e^(construct_float(13845666016792262747)*x)+construct_float(13873168701933081002)*x^(21/2)*%e^(construct_float(13844396881274955196)*x)+construct_float(4624476436427395444)*x^(25/2)*%e^(construct_float(13842717959610906452)*x))^2,x,0,inf),
1974 integrate (x^4*(construct_float(13851526180082158044)*sqrt(x)*%e^(construct_float(13851818867533264930)*x)+construct_float(4655653753038828547)*x^(5/2)*%e^(construct_float(13850206842258521600)*x)+construct_float(13892002975252393992)*x^(9/2)*%e^(construct_float(13848929967253417878)*x)+construct_float(4670725636947176122)*x^(13/2)*%e^(construct_float(13847268271141981213)*x)+construct_float(13886231423995319732)*x^(17/2)*%e^(construct_float(13845666016792262747)*x)+construct_float(4646523490691061153)*x^(21/2)*%e^(construct_float(13844396881274955196)*x)+construct_float(13844865077416599720)*x^(25/2)*%e^(construct_float(13842717959610906452)*x)+construct_float(4587941619969479187)*x^(29/2)*%e^(construct_float(13841125416961970319)*x))^2,x,0,inf),
1975 integrate (x^4*(construct_float(4624194709066097849)*sqrt(x)*%e^(construct_float(13851818867533264930)*x)+construct_float(13875102515253259746)*x^(5/2)*%e^(construct_float(13850206842258521600)*x)+construct_float(4664672840878753682)*x^(9/2)*%e^(construct_float(13848929967253417878)*x)+construct_float(13889968737336596784)*x^(13/2)*%e^(construct_float(13847268271141981213)*x)+construct_float(4658852307357948608)*x^(17/2)*%e^(construct_float(13845666016792262747)*x)+construct_float(13866104190951107225)*x^(21/2)*%e^(construct_float(13844396881274955196)*x)+construct_float(4617836750785529516)*x^(25/2)*%e^(construct_float(13842717959610906452)*x)+construct_float(13808553348881729110)*x^(29/2)*%e^(construct_float(13841125416961970319)*x)+construct_float(4543770050112687942)*x^(33/2)*%e^(construct_float(13839863974021311716)*x))^2,x,0,inf),
1976 integrate (x^4*(construct_float(13843560353427925598)*sqrt(x)*%e^(construct_float(13851818867533264930)*x)+construct_float(4647735622713936005)*x^(5/2)*%e^(construct_float(13850206842258521600)*x)+construct_float(13884036655858068600)*x^(9/2)*%e^(construct_float(13848929967253417878)*x)+construct_float(4662430308450520543)*x^(13/2)*%e^(construct_float(13847268271141981213)*x)+construct_float(13878145371141373487)*x^(17/2)*%e^(construct_float(13845666016792262747)*x)+construct_float(4638760849903044964)*x^(21/2)*%e^(construct_float(13844396881274955196)*x)+construct_float(13837200460116180214)*x^(25/2)*%e^(construct_float(13842717959610906452)*x)+construct_float(4581361216757031604)*x^(29/2)*%e^(construct_float(13841125416961970319)*x)+construct_float(13764215660170911913)*x^(33/2)*%e^(construct_float(13839863974021311716)*x)+construct_float(4491442966451646454)*x^(37/2)*%e^(construct_float(13838167931213425198)*x))^2,x,0,inf),
1977 0);
1978 0;
1979
1980 /*
1981  * SF bug #4265: integrate(sin(nn*x)*sin(mm*x),x,0,2*%pi) always zero.
1982  * We have similar issues with cos(nn*x)*cos(mm*x) and sin(nn*x)*cos(mm*x)
1983  */
1984 (assume(equal(nn, mm)),declare(nn, integer, mm, integer), 0);
1985 0;
1986 integrate(sin(nn*x)*sin(mm*x),x);
1987 x/2-sin(2*nn*x)/(4*nn);
1988
1989 integrate(cos(nn*x)*cos(mm*x),x);
1990 sin(2*nn*x)/(4*nn)+x/2;
1991
1992 integrate(cos(nn*x)*sin(mm*x),x);
1993 -cos(nn*x)^2/(2*nn);
1994
1995 integrate(sin(nn*x)*sin(mm*x),x,0,2*%pi);
1996 %pi;
1997
1998 integrate(cos(nn*x)*cos(mm*x),x, 0, 2*%pi);
1999 %pi;
2000
2001 (reset(), 0);
2002 0;