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 (6*%i*log(2)+6*%pi)/18;
 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 doens'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(abs(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 */
1356 /* that would require handling the discontinuity in
1357    gamma_incomplete(0, exp(%i*x))                               */
1358 integrate(exp(cos(x))*cos(sin(x)),x,0,2*%pi);
1359 'integrate(exp(cos(x))*cos(sin(x)),x,0,2*%pi);
1360
1361 /* Bug ID: 3039452 - integrate(sqrt(t^c)/(t*(b*t^c+a)),t) hangs
1362  */
1363 /* These assumptions are already present. Consider to delete it more early.
1364 assume(a>0,b>0);
1365 [a>0,b>0]; */
1366 integrate(sqrt(t^c)/(t*(b*t^c+a)),t);
1367 2*atan(sqrt(b)*%e^(c*log(t)/2)/sqrt(a))/(sqrt(a)*sqrt(b)*c);
1368 forget(a>0,b>0);
1369 [a>0,b>0];
1370
1371 /* Bug ID: 3045559 integrate(exp(-u^2), u, minf, x) => incorrect gamma_incomple
1372  */
1373 integrate(exp(-u^2),u,minf,x);
1374 sqrt(%pi)*erf(x)/2+sqrt(%pi)/2;
1375
1376 /* integration error - ID: 3085498 */
1377 integrate(sqrt(1-cos(t)),t,0,2*%pi);
1378 2^(5/2);
1379
1380 /* incorrect integration - ID: 3153533 */
1381 assume(x>2);
1382 [x>2];
1383 integrate(1/(x-z),z,1,x-1);
1384 log(x-1);
1385 forget(x>2);
1386 [x>2];
1387
1388 /* integrate(x^2*exp(x)/(1+exp(x))^2, x, minf, inf); - ID: 3158526 */
1389 integrate(x^2*exp(x)/(1+exp(x))^2, x, minf, inf);
1390 %pi^2/3;
1391
1392 /* integrate(abs(sin(x)),x,0,2*%pi) wrong result - ID: 3165872 */
1393 integrate(abs(sin(x)),x,0,2*%pi);
1394 'integrate(abs(sin(x)),x,0,2*%pi);
1395
1396 /* integrate(sin(2x)atan(sin(x)),x) - ID: 3199708 */
1397 integrate(sin(2*x)*atan(sin(x)),x,0,%pi/2);
1398 %pi/2 - 1;
1399
1400 /* wrong integration answer - ID: 2989983 */
1401 integrate(cos(w+T)/(1+1/2*cos(T))^2,T,0,2*%pi);
1402 -8*%pi*cos(w)/3^(3/2);
1403
1404 /* wrong sign for integral of e^(-1/x^2) - ID: 3211937 */
1405 integrate(%e^(-1/x^2),x,0,1);
1406 %e^-1*(sqrt(%pi)*%e*erf(1)+1)-sqrt(%pi);
1407
1408 /* definite integration introduces xy variable - ID: 3292707 */
1409 integrate((log((2-x)/2)+log(2))/(1-x), x, 0, 1);
1410 li[2](2)-(6*log(2)^2-6*log(-2)*log(2)+%pi^2)/6;
1411
1412 /* Incorrect integral of log(1+a/(x*t)^2) - ID: 3291160 */
1413 integrate(log(1+7/(x^2)),x,1,inf);
1414 (14*atan(sqrt(7))-sqrt(7)*log(8))/sqrt(7);
1415
1416 /* Wrong sign in exponential integral - ID: 3517785 */
1417 integrate(1/(%e^(2*t)*sqrt(1-1/%e^(2*t))), t, 0, inf);
1418 1;
1419
1420 /* integral from minf to inf of an absolute value combo fails - ID: 3313531 */
1421 integrate(abs(x - 1) + abs(x + 1) - 2*abs(x),x,minf,inf);
1422 2;
1423
1424 /* definite integration - ID: 3315476 */
1425 integrate(1/2^x, x, 0, 1);
1426 1/(2*log(2));
1427
1428 /* discontinuities in integral */
1429 /* error in integrating exp(-x)*sinh(sqrt(x)) - ID: 3292033 */
1430 integrate(exp(sqrt(x)-x), x, 0, inf);
1431 %e^(1/4)*gamma_incomplete(1,1/4)-%e^(1/4)*gamma_incomplete(1/2,1/4)/2
1432                                        +%e^(1/4)*(sqrt(%pi)+2)/2
1433                                        +%e^(1/4)*sqrt(%pi)/2-%e^(1/4);
1434 /* equivalent to:
1435    %e^(1/4)*gamma_incomplete(1,1/4)-%e^(1/4)*gamma_incomplete(1/2,1/4)/2
1436                                        +%e^(1/4)*sqrt(%pi);
1437 */
1438
1439 /* integrate(log(t)*log(t+1),t,0,1) gives Lisp error - ID: 3381012 */
1440 integrate(log(t)*log(t+1),t,0,1);
1441 -2*log(2)-%pi^2/12+2;
1442
1443 /* integrate(log(2*sin(t/2)),t,0,%pi) contains false - ID: 3381037 */
1444 integrate(log(2*sin(t/2)),t,0,%pi);
1445 2*((8*%i*li[2](%i)+8*%i*li[2](-%i)-4*%pi*log(2)+%i*%pi^2)/8
1446  +%pi*log(2)/2-%i*%pi^2/12);
1447 /* above result is really zero */
1448
1449 /* integrate error - ID: 3388801 */
1450 integrate( sqrt((x-474)^2/(107669-(x-474)^2)+1), x, 181, 474);
1451 sqrt(107669)*%i*log(2*sqrt(107669)*%i)
1452  - sqrt(107669)*%i*log(4*sqrt(5455)*%i-586);
1453
1454 /* integrate erf fails - ID: 3454370 */
1455 integrate( erf(x+a)-erf(x-a), x, minf, inf);
1456 4*a;
1457
1458 /* Wrong result for definite integral - ID: 3538167 */
1459 integrate(cos(3*x)/(5-4*cos(x)),x,0,2*%pi);
1460 %pi/12;
1461
1462 /* #2575 Integration error: integrate(sqrt(k-k*cos(2*x)), x) */
1463 integrate(sqrt(cos(x)+1), x, -%pi, %pi);
1464 2^(5/2);
1465
1466 /* #2776 Error when integrate sqrt */
1467 integrate ((1+exp(t))/sqrt(t+exp(t)), t, 0, 1);
1468 2*sqrt(%e+1)-2;
1469
1470 /* #2862 Incorrect result for integrate(u/(u+1)^2,u,0,inf); */
1471 errcatch(integrate(u/(u+1)^2,u,0,inf));
1472 [];
1473
1474 error;
1475 ["defint: integral is divergent."];
1476
1477
1478 /* #2903 Wrong integration over 1/(x+1) with bounds 0 to inf, since version 5.34.0 */
1479 errcatch(integrate(1/(x+1), x,0,inf));
1480 [];
1481
1482 error;
1483 ["defint: integral is divergent."];
1484
1485
1486 /* #2919 Definite integration broken: integrate(1/(x^2), x, -inf, inf) gives zero  */
1487 errcatch(integrate(1/(x^2), x, -inf, inf));
1488 [];
1489
1490 error;
1491 ["defint: integral is divergent."];
1492
1493 /* other tests mentioned in #2919 */
1494
1495 errcatch (integrate(1/x, x, 0, 1));
1496 [];
1497
1498 error;
1499 ["defint: integral is divergent."];
1500
1501 integrate(1/x, x, -1, 1);
1502 0;
1503 /* expected output: Principal Value -- dunno how to capture that */
1504
1505 errcatch (integrate(1/x, x, 0, inf));
1506 [];
1507
1508 error;
1509 ["defint: integral is divergent."];
1510
1511 errcatch (integrate(1/x, x, 1, inf));
1512 [];
1513
1514 error;
1515 ["defint: integral is divergent."];
1516
1517 integrate(1/x, x, -inf, inf);
1518 0;
1519 /* expected output: Principal Value -- dunno how to capture that */
1520
1521 errcatch (integrate(1/x^2, x, 0, 1));
1522 [];
1523
1524 error;
1525 ["defint: integral is divergent."];
1526
1527 errcatch (integrate(1/x^2, x, -1, 1));
1528 [];
1529
1530 error;
1531 ["defint: integral is divergent."];
1532
1533 errcatch (integrate(1/x^2, x, 0, inf));
1534 [];
1535
1536 error;
1537 ["defint: integral is divergent."];
1538
1539 /* #2987 Some divergent integrals give error, some don't */
1540 errcatch (integrate(exp(x)/x, x, 0, 1));
1541 [];
1542
1543 error;
1544 ["defint: integral is divergent."];
1545
1546 integrate(1/x^2, x, 1, inf);
1547 1;
1548
1549 /* SF # 2633: ev(integrate,numer) gives strange result */
1550 (foo : ev(integrate(sin(2*%pi*x),x,0,1),numer), 0);
1551 0;
1552
1553 freeof (false, foo);
1554 true;
1555
1556 float (foo);
1557 0.0;
1558
1559 /* Maxima used to call these integrals divergent even though it would yield 0
1560  * when both limits were inf.
1561  */
1562 integrate(exp(-x^2), x, minf, minf);
1563 0;
1564
1565 integrate(const, x, minf, minf);
1566 0;
1567
1568 /* Maxima used to give the bogus result "-false" for this integral even though
1569  * it would yield 0 when both limits were inf.
1570  */
1571 integrate(erf(x), x, minf, minf);
1572 0;
1573
1574 /* From the mailing list on 29 Aug 2015. Maxima was leaving integrate(0,x) in the result.
1575  */
1576 (reset(integration_constant_counter), integrate(x=0,x));
1577 x^2/2 = %c1;
1578
1579 /* SourceForge bug #3017 - "integrate_use_rootsof" leads to wrong result */
1580 %rnum:0;
1581 0;
1582 integrate(1/(x^(1/3)+x+1),x),integrate_use_rootsof:true;
1583 3*'lsum((%r1^2*log(x^(1/3)-%r1))/(3*%r1^2+1),%r1,
1584               rootsof(%r1^3+%r1+1,%r1));
1585
1586 /* Bug #2507: Strange non-evaluation of integral
1587  *
1588  * The problem was that the second integral was not evaluating
1589  * while the first one was.  The assume is also important because
1590  * the integral would evaluate if answered interactively.
1591  */
1592 (assume(n>0),0);
1593 0;
1594
1595 integrate(1/sqrt(1+x^2*n),x,1,2);
1596 asinh(2*sqrt(n))/sqrt(n)-asinh(sqrt(n))/sqrt(n);
1597
1598 integrate(1/sqrt(1+x^2/n),x,1,2);
1599 asinh(2/sqrt(n))*sqrt(n)-asinh(1/sqrt(n))*sqrt(n);
1600
1601 (forget(n>0),0);
1602 0;
1603
1604 /* mailing list 2015-11-23: "Change in behavior for taylor/integrate - a bug?" */
1605
1606 (kill(f, x, c, deltax), trunc(integrate(taylor(f(x), x, c, 6), x, c-deltax, c+deltax)));
1607 (deltax^7*('at('diff(f(x),x,6),x = c))+42*deltax^5*('at('diff(f(x),x,4),x = c))
1608                                       +840*deltax^3*('at('diff(f(x),x,2),x = c))+5040*f(c)*deltax)
1609                                        /2520$
1610
1611 integrate (taylor (f(x), x, 0, 1), x, a, b);
1612 ((b^2-a^2)*('at('diff(f(x),x,1),x = 0))+2*f(0)*b-2*f(0)*a)/2$
1613
1614 /* mailing list 2015-11-08: "un-nerving behavior of numer... yikes!" */
1615
1616 /* disable this test pending resolution of SF bug #3097 !!
1617 integrate(x*sin(%pi/2*x), x,0,1 ),    numer;
1618 0.4052847345693511;
1619  */
1620
1621 integrate(x*sin(1.5*x), x,0,1 ),    numer;
1622 0.3961729707122222;
1623
1624 integrate(x*cos(1.5*x), x,0,1 ),    numer;
1625 0.2519909695883491;
1626
1627 integrate(x*(cos(2.5*x) + 1.25), x, 0, 1), numer;
1628 0.5762058791540733;
1629
1630 integrate(x*(cos(2.5*x^2 - 1) + 1.25), x, 0, 1), numer;
1631 0.9927931942823904;
1632
1633 integrate (sin(4*asin(x)),x, 0.5, 0.75), numer;
1634 0.09667085773175951;
1635
1636 /* disable these next few tests, since these fail (with a floating point error)
1637  * in the absence of float-to-bigfloat promotion, which was recently reverted.
1638  *
1639 /* These next 4 examples yield big messy rational expressions.
1640  * It's maybe not so great to have to squash them (with bfloat and rectform)
1641  * but at least we get a result; these examples cause a floating-point error
1642  * in previous versions.
1643  */
1644 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)));
1645 (- 6.358161510436547b-2*%i) - 2.859595163448958b-2;
1646
1647 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)));
1648 6.096077978946073b-2*%i - 3.382504333513841b-2;
1649
1650 rectform (bfloat (ev (integrate (x^2*%e^(30*%i*x-200*%i)*cos(2*x),x, 1.25, 1.5), numer)));
1651 (- 6.35816151043802b-2*%i) - 2.859595163448362b-2;
1652
1653 rectform (bfloat (ev (integrate (x^2*%e^(30*%i*x-20*%i)*cos(2*x),x, 1.25, 1.5), numer)));
1654 6.096077978946119b-2*%i - 3.382504333513988b-2;
1655  *
1656  */
1657
1658
1659 /* bug #3114 */
1660 integrate(sec(x)/sqrt(tan(x)),x);
1661 'integrate(sec(x)/sqrt(tan(x)),x);
1662
1663 /* bug #3115 */
1664 integrate((x*sqrt(x^3-4*x))/(x^2-4),x);
1665 'integrate((x*sqrt(x^3-4*x))/(x^2-4),x);
1666
1667 /* SF bug #3144: "stackoverflow in integral" */
1668
1669 integrate(sin(k*x)/x*erf(x^2),x,0,inf);
1670 'integrate(sin(k*x)/x*erf(x^2),x,0,inf);
1671
1672 /* Borwein integrals; see: https://en.wikipedia.org/wiki/Borwein_integral */
1673
1674 integrate(sin(x)/x,x,0,inf)$
1675 %pi/2 $
1676
1677 integrate((sin(x)/x) * (sin(x/3)/(x/3)),x,0,inf)$
1678 %pi/2 $
1679
1680 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)),x,0,inf)$
1681 %pi/2 $
1682
1683 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)) * (sin(x/7)/(x/7)),x,0,inf)$
1684 %pi/2 $
1685
1686 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)$
1687 %pi/2 $
1688
1689 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)) * (sin(x/7)/(x/7)) * (sin(x/9)/(x/9))
1690                      * (sin(x/11)/(x/11)),x,0,inf)$
1691 %pi/2 $
1692
1693 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)) * (sin(x/7)/(x/7)) * (sin(x/9)/(x/9))
1694                      * (sin(x/11)/(x/11)) * (sin(x/13)/(x/13)),x,0,inf)$
1695 %pi/2 $
1696
1697 integrate((sin(x)/x) * (sin(x/3)/(x/3)) * (sin(x/5)/(x/5)) * (sin(x/7)/(x/7)) * (sin(x/9)/(x/9))
1698                      * (sin(x/11)/(x/11)) * (sin(x/13)/(x/13)) * (sin(x/15)/(x/15)),x,0,inf)$
1699 467807924713440738696537864469*%pi/935615849440640907310521750000 $
1700
1701 /* Bug #2314: Two different results for an integral */
1702
1703 (assume (0 < x, x < 1, 0 < y, y < 1), 0)$
1704 0$
1705
1706 integrate (integrate (sqrt (4*x^2 - y^2), x, y, 1), y, 0, 1)$
1707 (2*sqrt(3)*%pi+9)/(2*3^(5/2))$
1708
1709 integrate (integrate (sqrt (4*x^2 - y^2), y, 0, x), x, 0, 1)$
1710 (2*%pi+3^(3/2))/18$
1711
1712 (forget (0 < x, x < 1, 0 < y, y < 1), 0)$
1713 0$
1714
1715 /* Bug #2295: Failure to evaluate definite integral */
1716
1717 integrate (1 / (sqrt (x) * (1 + sqrt (x))^2), x, 1, 4)$
1718 1/3$
1719
1720 /* SF bug #3274: "quotient by zero" error in integrate */
1721
1722 (expr_tau: (%e^-(t/tau3)*(%e^(t/tau3+t/tau1)-%e^t*tau3))
1723       /((tau1*%e^(t/tau1)+%e^(t/tau1))
1724        *%e^(t/tau2)*tau3+%e^(t/tau2+t/tau1)),
1725  I_tau : integrate(expr_tau, t));
1726 -(%e^((-t/tau2)-t/tau1)*(tau1*tau2*tau3^2*%e^(t-t/tau3)
1727                         +%e^(t/tau1)
1728                          *(tau1*(tau2^2*(tau3-1)-tau2*tau3)
1729                           -tau2^2*tau3)))
1730  /(tau1^2*(tau2*(tau3^2-tau3)-tau3^2)
1731   +tau1*((-tau3^2)-tau3-tau2)+tau2*((-tau3^2)-tau3))$
1732
1733 is(equal(diff(I_tau, t), expr_tau));
1734 true;
1735
1736 (expr_c : subst([tau1=c1, tau2=c2, tau3=c3], expr_tau),
1737  I_c : integrate (expr_c, t));
1738 (-(c1*c2*c3^2*%e^((-t/c3)-t/c2-t/c1+t))
1739  /(c1^2*(c2*(c3^2-c3)-c3^2)+c1*((-c3^2)-c3-c2)+c2*((-c3^2)-c3)))
1740  -(c2*c3*%e^-(t/c2))/(c1*c3^2+c3^2+c3)$
1741
1742 is(equal(diff(I_c, t), expr_c));
1743 true;
1744
1745 constantp (ratsimp (subst([c1=tau1, c2=tau2, c3=tau3], I_c) - I_tau));
1746 true;
1747
1748 /* following was incorrect (returned 0) at some point in recent past
1749  * ensure that it stays correct
1750  * from sage-support 2017-02-22
1751  */
1752
1753 errcatch (integrate(x/(1+x^2),x,0,inf));
1754 [];
1755
1756 /* Bug #3368: integrate('limit(...),...) internal error */
1757
1758 /* Original test cases from bug report: */
1759
1760 integrate('limit(x,x,0),x);
1761 x*'limit(x,x,0);
1762
1763 integrate('limit(x,x,y),x);
1764 x*'limit(x,x,y);
1765
1766 integrate('limit(x,x,x),x);
1767 'integrate('limit(x,x,x),x);
1768
1769 /* Some more: */
1770
1771 integrate('at(x,x=1), x);
1772 x*'at(x,x=1);
1773
1774 integrate('sum(x,x,1,1), x);
1775 x*'sum(x,x,1,1);
1776
1777 integrate('product(x,x,1,1), x);
1778 x*'product(x,x,1,1);