tests/rtestsum.mac

   1 kill (all);
   2 done;
   3
   4 (assume_pos_save: assume_pos, assume_pos: true);
   5 true;
   6
   7 /* Verify that we at least get errors when passing the wrong number
   8  * of arguments to sum and product.
   9  *
  10  * The ev tests are significant.  In the 2 and 3 arg cases, we used
  11  * to get bogus results.  In the case of more than 4 args, the extra
  12  * args used to be effectively ignored.  (In the 0 and 1 arg cases,
  13  * we got obscure error messages about assigning to false.  The tests
  14  * below don't do any checking about *why* we got an error.)
  15  */
  16
  17 errcatch (sum ());
  18 [];
  19
  20 errcatch (ev (sum ()));
  21 [];
  22
  23 errcatch (sum (f (x)));
  24 [];
  25
  26 errcatch (ev (sum (f (x))));
  27 [];
  28
  29 errcatch (sum (f (x), x));
  30 [];
  31
  32 errcatch (ev (sum (f (x), x)));
  33 [];
  34
  35 errcatch (sum (f (x), x, 0));
  36 [];
  37
  38 errcatch (ev (sum (f (x), x, 0)));
  39 [];
  40
  41 errcatch (sum (f (x), x, 0, 1, 2));
  42 [];
  43
  44 errcatch (ev (sum (f (x), x, 0, 1, 2)));
  45 [];
  46
  47 errcatch (product ());
  48 [];
  49
  50 errcatch (ev (product ()));
  51 [];
  52
  53 errcatch (product (f (x)));
  54 [];
  55
  56 errcatch (ev (product (f (x))));
  57 [];
  58
  59 errcatch (product (f (x), x));
  60 [];
  61
  62 errcatch (ev (product (f (x), x)));
  63 [];
  64
  65 errcatch (product (f (x), x, 1));
  66 [];
  67
  68 errcatch (ev (product (f (x), x, 1)));
  69 [];
  70
  71 errcatch (product (f (x), x, 1, 2, 3));
  72 [];
  73
  74 errcatch (ev (product (f (x), x, 1, 2, 3)));
  75 [];
  76
  77 /* Known failure: concat (q, k) => gensym
  78  */
  79 sum (concat (q, k), k, a, b);
  80 (b - a + 1) * qk;
  81
  82 /* Known failure: concat (q, k) => gensym
  83  */
  84 product (concat (q, k), k, a, b);
  85 qk ^ (b - a + 1);
  86
  87 /* Known failure: causes asksign to ask an unanswerable question about index variable;
  88  * this is a bug in integrate and/or asksign.
  89 (assume (b > 0, a > 0, b > a), sum (integrate (1/x^k, x, a, b), k, 1, n));
  90 integrate (1/x, x, a, b) + sum (integrate (1/x^k, x, a, b), k, 2, n);
  91  */
  92
  93 /* Known failure: causes asksign to ask an unanswerable question about index variable;
  94  * this is a bug in integrate and/or asksign.
  95 product (integrate (1/x^k, x, a, b), k, 1, n);
  96 integrate (1/x, x, a, b) * product (integrate (1/x^k, x, a, b), k, 2, n);
  97  */
  98
  99 block ([prederror: false], sum(if k < 1 then a else b,k,0,n));
 100 '( 'sum(if k < 1 then a else b,k,0,n));
 101
 102 /* see SF bug report # 625278 */
 103
 104 'sum(binomial(2,2-k)-binomial(2,1-k),k,1,2),simpsum;
 105 2;
 106
 107 'sum(binomial(2,2-k)-binomial(2,1-k),k,1,2),sum;
 108 2;
 109
 110 'sum(binomial(x,2-k)-binomial(x,1-k),k,1,2),simpsum;
 111 x;
 112
 113 'sum(binomial(x,2-k)-binomial(x,1-k),k,1,2),sum;
 114 x;
 115
 116 sum (f(k) + 1, k, 1, n), simpsum;
 117 n + 'sum (f(k), k, 1, n);
 118
 119 /* further examples in the same vein: part of summand depends on index, part doesn't. */
 120
 121 sum (sin(g(j)) - %pi^j + cos(h(k)) - %e^l, j, 1, m), simpsum;
 122 m*cos(h(k)) - m*%e^l + sum (sin(g(j)), j, 1, m) - ((%pi^(m + 1) - %pi)/(%pi - 1));
 123
 124 sum (sum (f(j) + g(k), j, 1, m), k, 1, n), simpsum;
 125 m * 'sum (g(k), k, 1, n) + n * 'sum (f(j), j, 1, m);
 126
 127 sum (sum (f(j) + g(k) - %pi^l + 2, j, 1, m), k, 1, n), simpsum;
 128 m * 'sum (g(k), k, 1, n) + n * 'sum (f(j), j, 1, m) - m*n*%pi^l + 2*m*n;
 129
 130 sum (sum (sum (%pi + %phi - 1, i, 1, n_i), j, 1, n_j), k, 1, n_k);
 131 n_i * n_j * n_k * (%pi + %phi - 1);
 132
 133 sum (sum (sum (f(%pi, i) + g(%phi, j) - h(k), i, 1, n_i), j, 1, n_j), k, 1, n_k);
 134 'sum ('sum ('sum (f(%pi, i) + g(%phi, j) - h(k), i, 1, n_i), j, 1, n_j), k, 1, n_k);
 135
 136 /* see SF bug report # 649428 */
 137
 138 (factsum3 (mt, ej) := sum ((-1)^(k + 1)/(k*mt^k)*sum((1 - l)^k - l^k, l, 1, ej),k,1,inf), t% : taylor (factsum3 (mt, ej), [mt, 0, 3, asymp]));
 139 'sum(1-2*l,l,1,ej)/mt-('sum(1-2*l,l,1,ej)/(2*mt^2))+'sum(-2*l^3+3*l^2-3*l+1,l,1,ej)/(3*mt^3);
 140
 141 %, simpsum, ratsimp;
 142 -((6*ej^2*mt^2-3*ej^2*mt+ej^4+ej^2)/(6*mt^3));
 143
 144 /* Known failure: interaction of taylor and simpsum
 145  */
 146 taylor (factsum3 (mt, ej), [mt, 0, 3, asymp]), simpsum;
 147 -(6*ej^2*mt^2-3*ej^2*mt+ej^4+ej^2)/(6*mt^3);
 148
 149 /* see SF bug report # 740134 */
 150
 151 (foo: i^2, [sum ('foo, i, 0, 2), sum (foo, i, 0, 2), sum (foo, i, 0, n)]);
 152 [3*foo, 5, 'sum (i^2, i, 0, n)];
 153
 154 /* summand is a function which has a side effect --
 155  * outcome here differs from recommendation by Stavros in 740134
 156  */
 157 block ([L:[]], f(i) := (L: append (L, [i]), i), S: sum (f (i), i, 1, 3), [S, L]);
 158 [6, [1, 2, 3]];
 159
 160 sum (integrate (x^i ,x), i, 0, 2);
 161 x^3/3 + x^2/2 + x;
 162
 163 sum (integrate (1/(x^i + 1), x), i, 0, 1);
 164 log(x+1) + x/2;
 165
 166 (f[i](x) := x^i, g[i](x) := x^i, h[i](x) := x^i, 0);
 167 0;
 168
 169 (f[i], g[i](t), 0);
 170 0;
 171
 172 sum (f[i](x), i, 0, n);
 173 'sum (x^i, i, 0, n);
 174
 175 sum (g[i](x), i, 0, n);
 176 'sum (x^i, i, 0, n);
 177
 178 sum (h[i](x), i, 0, n);
 179 'sum (x^i, i, 0, n);
 180
 181 /* see SF bug report # 817521 */
 182
 183 (kill (foo), foo (n):= block([aaa,m], modedeclare([m,n], fixnum), m:n, [ aaa[1], aaa[n], aaa[m], aaa[1]+aaa[2], sum(aaa[i],i,1,m) ] ), foo (2));
 184 [aaa[1],aaa[2],aaa[2],aaa[2]+aaa[1], aaa[2]+aaa[1]];
 185
 186 (translate(foo), 0);
 187 0;
 188
 189 foo (2);
 190 [aaa[1],aaa[2],aaa[2],aaa[2]+aaa[1], aaa[2]+aaa[1]];
 191
 192 /* see SF bug report # 851765 */
 193
 194 (apply (forget, facts ()), assume (i >= a), sum (i, i, a, b), facts ());
 195 [i >= a];
 196
 197 is (i >= a);
 198 true;
 199
 200 (assume (i < 0), sum (abs (i), i, 1, a));
 201 'sum (i, i, 1, a);
 202
 203 (apply (forget, facts ()), 0);
 204 0;
 205
 206 /* see SF bug report # 1007094 */
 207
 208 (u: 'sum((-1/(8*i+6)-1/(8*i+5)-2/(8*i+4)+4/ (8*i+1))/16^i,i,1,a), 0);
 209 0;
 210
 211 u, a=2, simpsum;
 212 1618091/196035840;
 213
 214 %, numer;
 215 0.0082540570132481894$
 216
 217 u, a=2, simpsum, numer;
 218 0.0082540570132481911$
 219
 220 u, a=2, simpsum, bfloat;
 221 8.25405701324819B-3;
 222
 223 /* see SF bug report # 1192935 */
 224
 225 (f(x) := x^2, b: 1/10, Ak: b*f(k*b), sum (Ak, k, 0, 9));
 226 57/200;
 227
 228 /* When A(k) is defined as A(k) := b*f(k*b) instead of define(A(k), b*f(k*b)),
 229  * A(k) evaluates to k^2/n^3 even if n is bound to something;
 230  * this is a consequence of Maxima's evaluation policies:
 231  * rhs of := is not evaluated, and variables (b in this case)
 232  * are evaluated just once.
 233  */
 234 (b:1/n, define (A(k), b*f(k*b)), A(k));
 235 k^2/n^3;
 236
 237 sum (A(k), k, 0, n-1);
 238 ('sum (k^2, k, 0, n-1))/n^3;
 239
 240 sum (A(k), k, 0, n-1), simpsum;
 241 (n+2*(n-1)^3+3*(n-1)^2-1)/(6*n^3);
 242
 243 [''%, sum (A(k), k, 0, n-1)], n=10;
 244 [57/200, 57/200];
 245
 246 (B(k) := k*f(k*b), B(k));
 247 k^3/n^2;
 248
 249 sum (B(k), k, 0, n-1), simpsum;
 250 ((n-1)^4 + 2*(n-1)^3 + (n-1)^2) / (4*n^2);
 251
 252 /* see SF bug report # 1363411 */
 253
 254 (kill (f), 'sum(1+f(k),k,1,2)), simpsum;
 255 2 + f(1) + f(2);
 256
 257 is(equal(block([gcd:subres], ratsimp((-x^3+1)/(2*x^5+2) + (-sqrt(5)+5)/(20*x^2+(-10*sqrt(5)-10)*x+20))),
 258          block([gcd:spmod],  ratsimp((-x^3+1)/(2*x^5+2) + (-sqrt(5)+5)/(20*x^2+(-10*sqrt(5)-10)*x+20)))));
 259 true;
 260
 261 /* General tests from here on */
 262
 263 sum (x, x, p, p + 2);
 264 p + (p + 1) + (p + 2);
 265
 266 product (x, x, p, p + 2);
 267 p * (p + 1) * (p + 2);
 268
 269 sum (x, x, %i, %i + 2);
 270 %i + (%i + 1) + (%i + 2);
 271
 272 product (x, x, %i, %i + 2);
 273 %i * (%i + 1) * (%i + 2);
 274
 275 sum (imagpart (1 + %i*k), k, 1, n), simpsum;
 276 (n^2 + n)/2;
 277
 278 product (imagpart (1 + %i*k), k, 1, n), simpproduct;
 279 n!;
 280
 281 (g(x) := q(x^2), sum (g(k), k, 1, n));
 282 'sum (q (k^2), k, 1, n);
 283
 284 (g(x) := q(x^2), product (g(k), k, 1, n));
 285 'product (q (k^2), k, 1, n);
 286
 287 sum (integrate (x^k, x, 0, 1), k, 1, n);
 288 'sum (1/(k + 1), k, 1, n);
 289
 290 product (integrate (x^k, x, 0, 1), k, 1, n);
 291 'product (1/(k + 1), k, 1, n);
 292
 293 (kill(b), sum ('concat (q, k), k, a, b));
 294 'sum ('concat (q, k), k, a, b);
 295
 296 sum ('concat (q, k), k, a, b), nouns, a=1, b=5;
 297 q5 + q4 + q3 + q2 + q1;
 298
 299 product ('concat (q, k), k, a, b);
 300 'product ('concat (q, k), k, a, b);
 301
 302 product ('concat (q, k), k, a, b), nouns, a=1, b=5;
 303 q5 * q4 * q3 * q2 * q1;
 304
 305 (assume (n>= 1), sum (if k <= n then a else b, k, 1, n));
 306 a * n$
 307
 308 product (if k <= n then a else b, k, 1, n);
 309 a ^ n$
 310
 311 sum (sin (%pi * k), k, 0, nn);
 312 0$
 313
 314 product (cos (2 * %pi * k), k, 0, nn);
 315 1$
 316
 317 sum (lambda ([x], x^i), i, 1, 3);
 318 lambda ([x], x) + lambda ([x], x^2) + lambda ([x], x^3);
 319
 320 sum (k: k^2, k, 1, 3);
 321 14;
 322
 323 sum (a, k, inf, minf);
 324 0;
 325
 326 (assume (i < 0, j < 0, n > 0), sum(sum(abs(i) + abs(j),i,1,j),j,1,n));
 327 'sum ('sum (i + j, i, 1, j), j, 1, n);
 328
 329 sum (sum (sqrt (i + j), i, 1, j), j, 1, n), simpsum;
 330 'sum ('sum (sqrt (i + j), i, 1, j), j, 1, n);
 331
 332 sum (1/k, k, rat(1), rat(n));
 333 'sum (1/k, k, 1, n);
 334
 335 sum (taylor (x^i, x, 0, 5), i, 0, 5);
 336 ''(taylor (x^5 + x^4 + x^3 + x^2 + x + 1,x,0,5));
 337
 338 (i: 888, sum (lambda ([i], i^2), i, 1, 3));
 339 3 * lambda ([i], i^2);
 340
 341 sum (lambda ([i], i^2), i, 1, n);
 342 n * lambda ([i], i^2);
 343
 344 (f(x) := sum (x, i, 1, 3), f(-x));
 345 -3*x;
 346
 347 (ak : k^2,  g(a,n) := sum(a,k,1,n), g(ak, 5));
 348 55;
 349
 350 (translate (g), g(ak,5));
 351 55;
 352
 353 sum (integrate (x^i, x), i, 0, n);
 354 'sum (x^(i+1) / (i+1), i, 0, n);
 355
 356 (assume (k < 0), assume (0 < n, m > n), sum (sqrt (k^2), k, n, m));
 357 'sum (k, k, n, m);
 358
 359 (assume (i >= n_1, i <= n_2), sum (z(i), i, n_1, n_2));
 360 'sum (z(i), i, n_1, n_2);
 361
 362 block ([prederror: false], [is (i >= n_1), is (i <= n_2)]);
 363 [true, true];
 364
 365 /* some examples adapted from describe("sum") */
 366
 367 (old_facts: facts (), foo: 42, bar: 12 + exp(7), baz: 1/quux, 0);
 368 0;
 369
 370 sum (foo^2, foo, 1, 7);
 371 140;
 372
 373 sum (aa[bar], bar, 1, 7);
 374 aa[7] + aa[6] + aa[5] + aa[4] + aa[3] + aa[2] + aa[1];
 375
 376 sum (aa(baz), baz, 1, 7);
 377 aa(7) + aa(6) + aa(5) + aa(4) + aa(3) + aa(2) + aa(1);
 378
 379 sum (aa(foo), foo, 1, n);
 380 'sum (aa(foo), foo, 1, n);
 381
 382 sum (2^bar + bar^2, bar, 0, n), simpsum;
 383 2^(n + 1) + (2*n^3 + 3*n^2 + n)/6 - 1;
 384
 385 sum (1/3^baz, baz, 1, inf), simpsum;
 386 1/2;
 387
 388 sum (foo^2, foo, 1, 4) * sum (1/foo^2, foo, 1, inf), simpsum;
 389 5*%pi^2;
 390
 391 (gg(bar) := 1/bar^2, sum (gg(bar)^2, bar, 1, inf));
 392 sum (1/bar^4, bar, 1, inf);
 393
 394 sum (gg(baz)^2, baz, 1, inf), simpsum;
 395 %pi^4/90;
 396
 397 [is (foo = 42), is (bar = 12 + exp(7)), is (baz = 1/quux)];
 398 [true, true, true];
 399
 400 (new_facts: facts(), is (sort (old_facts) = sort (new_facts)));
 401 true;
 402
 403 /* redo all of the tests, this time with product */
 404
 405 product (lambda ([x], x^i), i, 1, 3);
 406 lambda ([x], x) * lambda ([x], x^2) * lambda ([x], x^3);
 407
 408 product (k: k^2, k, 1, 3);
 409 36;
 410
 411 product (a, k, inf, minf);
 412 1;
 413
 414 (assume (i < 0, j < 0, n > 0), product(product(abs(i) + abs(j),i,1,j),j,1,n));
 415 'product ('product (i + j, i, 1, j), j, 1, n);
 416
 417 product (1/k, k, rat(1), rat(n));
 418 'product (1/k, k, 1, n);
 419
 420 product (taylor (x^i, x, 0, 5), i, 0, 5);
 421 ''(taylor (x^15,x,0,15));
 422
 423 (i: 888, product (lambda ([i], i^2), i, 1, 3));
 424 (lambda ([i], i^2))^3;
 425
 426 product (lambda ([i], i^2), i, 1, n);
 427 lambda([i],i^2)^n;
 428
 429 (f(x) := product (x, i, 1, 3), f(-x));
 430 -x^3;
 431
 432 (f(x) := x^2, b: 1/10, Ak: b*f(k*b), product (Ak, k, 1, 9));
 433 (9!)^2/1000^9;
 434
 435 (A(k) := k*f(k*b), b:1/n, A(k));
 436 k^3/n^2;
 437
 438 product (A(k), k, 0, n-1);
 439 ('product (k^3/n^2, k, 0, n-1));
 440
 441 (assume (n > 0, m > n), ev (product (k, k, n, m), simpproduct));
 442 m!/(n - 1)!;
 443
 444 (ak : k^2,  g(a,n) := product(a,k,1,n), g(ak, 5));
 445 (5!)^2;
 446
 447 (translate (g), g(ak,5));
 448 (5!)^2;
 449
 450 /* multiplicand is a function which has a side effect */
 451
 452 block ([L:[]], f(i) := (L: append (L, [i]), i), S: product (f (i), i, 1, 3), [S, L]);
 453 [6, [1, 2, 3]];
 454
 455 product (integrate (x^i ,x), i, 0, 2);
 456 x^3/3 * x^2/2 * x;
 457
 458 product (integrate (1/(x^i + 1), x), i, 0, 1);
 459 log(x+1) * x/2;
 460
 461 (f[i](x) := x^i, g[i](x) := x^i, h[i](x) := x^i, 0);
 462 0;
 463
 464 /* reference f[i] and g[i] -- see 740134 for the effect this has on previous defn of sum */
 465
 466 (f[i], g[i](t), 0);
 467 0;
 468
 469 product (f[i](x), i, 0, n);
 470 'product (x^i, i, 0, n);
 471
 472 product (g[i](x), i, 0, n);
 473 'product (x^i, i, 0, n);
 474
 475 product (h[i](x), i, 0, n);
 476 'product (x^i, i, 0, n);
 477
 478 product (integrate (x^i, x), i, 0, n);
 479 'product (x^(i+1) / (i+1), i, 0, n);
 480
 481 (assume (k < 0), assume (0 < n, m > n), ev (product (sqrt (k^2), k, n, m), simpproduct));
 482 m!/(n - 1)!;
 483
 484 (assume (i >= n_1, i <= n_2), product (z(i), i, n_1, n_2));
 485 'product (z(i), i, n_1, n_2);
 486
 487 block ([prederror: false], [is (i >= n_1), is (i <= n_2)]);
 488 [true, true];
 489
 490 /* some examples adapted from describe("sum") */
 491
 492 (old_facts: facts (), foo: 42, bar: 12 + exp(7), baz: 1/quux, 0);
 493 0;
 494
 495 product (foo^2, foo, 1, 7);
 496 (7!)^2;
 497
 498 product (aa[bar], bar, 1, 7);
 499 aa[7] * aa[6] * aa[5] * aa[4] * aa[3] * aa[2] * aa[1];
 500
 501 product (aa(baz), baz, 1, 7);
 502 aa(7) * aa(6) * aa(5) * aa(4) * aa(3) * aa(2) * aa(1);
 503
 504 product (aa(foo), foo, 1, n);
 505 product (aa(foo), foo, 1, n);
 506
 507 product (2^bar + bar^2, bar, 0, n), simpsum;
 508 product (2^bar + bar^2, bar, 0, n);
 509
 510 product (1/3^baz, baz, 1, inf), simpsum;
 511 product (1/3^baz, baz, 1, inf);
 512
 513 product (foo^2, foo, 1, 4) * product (1/foo^2, foo, 1, inf), simpsum;
 514 (4!)^2 * 'product (1/'foo^2, 'foo, 1, inf);
 515
 516 (gg(bar) := 1/bar^2, product (gg(bar)^2, bar, 1, inf));
 517 'product (1/'bar^4, 'bar, 1, inf);
 518
 519 product (gg(baz)^2, baz, 1, inf), simpsum;
 520 'product (1/'baz^4, 'baz, 1, inf);
 521
 522 [is (foo = 42), is (bar = 12 + exp(7)), is (baz = 1/quux)];
 523 [true, true, true];
 524
 525 (new_facts: facts(), is (sort (old_facts) = sort (new_facts)));
 526 true;
 527
 528 /*---- new tests from bw ------*/
 529
 530 sum(0,i,a,b);
 531 0$
 532
 533 sum(0,i,minf,inf);
 534 0$
 535
 536 sum(0,i,42,inf);
 537 0$
 538
 539 sum(sqrt(3),i,0,inf);
 540 inf$
 541
 542 sum(sqrt(3),i,-1932,inf);
 543 inf$
 544
 545 sum(5,i,minf, inf);
 546 inf$
 547
 548 sum(-sqrt(3),i,0,inf);
 549 minf$
 550
 551 sum(sqrt(3),i,1,3);
 552 3 * sqrt(3)$
 553
 554 sum(sqrt(3),i,-1,1);
 555 3 * sqrt(3)$
 556
 557 sum(p(i + 1) - p(i),i,1,5);
 558 p(6) - p(1)$
 559
 560 product(1,i,a,b);
 561 1$
 562
 563 product(1,i,minf,inf);
 564 1$
 565
 566 product(1,i,-42,inf);
 567 1$
 568
 569 product(%pi,i,1,3);
 570 %pi^3$
 571
 572 product(%pi,i,-1,1);
 573 %pi^3$
 574
 575 product(%pi,i,0,inf);
 576 inf$
 577
 578 product(%pi - 3,i,0,inf);
 579 0$
 580
 581 product(%pi - 3,i,minf,inf);
 582 0$
 583
 584 product(-%pi,i,0,inf);
 585 infinity$
 586
 587 product(0,i,0,n);
 588 0$
 589
 590 product(-%pi,i,minf,inf);
 591 infinity$
 592
 593 product((1 + %i)/sqrt(2),k,1,inf);
 594 und$
 595
 596 product(-1,k,1,inf);
 597 und$
 598
 599 product(a,i,0,inf);
 600 'product(a,i,0,inf);
 601
 602 product(p(i + 1) / p(i),i,1,5);
 603 p(6) / p(1)$
 604
 605 sum(lambda([x],x^i),i,1,3);
 606 lambda([x],x) + lambda([x],x^2) + lambda([x],x^3);
 607
 608 product(lambda([x],x^i),i,1,3);
 609 lambda([x],x) * lambda([x],x^2) * lambda([x],x^3);
 610
 611 sum([1,k],k,1,5);
 612 [5,15];
 613
 614 sum([1,k],k,1,n),simpsum;
 615 [n, (n^2 + n)/2];
 616
 617 (old_facts : facts(),0);
 618 0$
 619
 620 errcatch(sum(1/k,k,-1,1));
 621 [];
 622
 623 (new_facts: facts(), is (sort (old_facts) = sort (new_facts)));
 624 true;
 625
 626 (old_facts : facts(),0);
 627 0$
 628
 629 errcatch(product(1/k,k,-1,1));
 630 [];
 631
 632 (new_facts: facts(), is (sort (old_facts) = sort (new_facts)));
 633 true;
 634
 635 sum(diff(log(x),x,k),k,1,2);
 636 1/x - 1/x^2$
 637
 638 product(diff(log(x),x,k),k,1,2);
 639 -1/x^3;
 640
 641 sum(x^k * at(diff(sin(x),x,k)/k!,x=0),k,0,5) - taylor(sin(x),x,0,5);
 642 ''(taylor(0,x,0,5));
 643
 644 product(%i-x,k,1,inf);
 645 'product(%i-x,k,1,inf);
 646
 647 sum(%i,k,0,-1);
 648 0$
 649
 650 sum(%i,k,inf,minf);
 651 0$
 652
 653 sum(%i,k, 0, minf);
 654 0$
 655
 656 product(%i,k,0,-1);
 657 1$
 658
 659 product(%i,k,inf,minf);
 660 1$
 661
 662 product(%i,k, 0, minf);
 663 1$
 664
 665 sum(p,p,p,p);
 666 p$
 667
 668 product(p,p,p,p);
 669 p$
 670
 671 sum(k,k,p,p+1);
 672 2*p+1$
 673
 674 sum(k,k,p,p-1);
 675 0$
 676
 677 sum(k,k,%i,%i+1);
 678 2*%i + 1$
 679
 680 product(k,k,%i,%i+1);
 681 %i * (%i + 1)$
 682
 683 product(k,k,p,p+1);
 684 p*(p+1)$
 685
 686 product(k,k,p,p-1);
 687 1$
 688
 689 sum(k,k,1,entier(4/3));
 690 1$
 691
 692 sum(k,k,entier(4/3),1);
 693 1$
 694
 695 sum(a,k,entier(4/3),entier(5/2));
 696 2*a$
 697
 698 product(k,k,1,entier(4/3));
 699 1$
 700
 701 product(k,k,entier(4/3),1);
 702 1$
 703
 704 product(a,k,entier(4/3),entier(5/2));
 705 a^2$
 706
 707 sum(q(k),k,n^2 + 1,-n^2 - 1); /* It's an empty sum */
 708 0$
 709
 710 product(q(k),k,n^2 + 1,-n^2 - 1);
 711 1$
 712
 713 sum(concat(q,k),k,1,3);
 714 q1 + q2 + q3$
 715
 716 (s : sum('concat(q,k),k,1,n),0);
 717 0$
 718
 719 (s : subst(n=3,s), ev(s,nouns));
 720 q1 + q2 + q3$
 721
 722 sum(product(concat(q,i),i,1,j),j,1,3);
 723 q1 + q1 * q2 + q1 * q2 * q3$
 724
 725 product(concat(q,k),k,1,3);
 726 q1 * q2 * q3$
 727
 728 (s : product('concat(q,k),k,1,n),0);
 729 0$
 730
 731 (s : subst(n=3,s), ev(s,nouns));
 732 q1 * q2 * q3$
 733
 734 (a : k^2,0);
 735 0$
 736
 737 sum(realpart(a + %i * k),k,1,3);
 738 14$
 739
 740 (assume(k < 0),0);
 741 0$
 742
 743 sum(sum(abs(k),k,1,n),n,1,m);
 744 sum(sum(k,k,1,n),n,1,m)$
 745
 746 sum(sum(k,k,1,i),i,1,n),simpsum,ratsimp;
 747 ''(ratsimp(n^3/6+n^2/2+n/3))$
 748
 749 sum(sum(abs(k),k,1,n) + sum(abs(k),k,1,n),n,1,m);
 750 2 * sum(sum(k,k,1,n),n,1,m)$
 751
 752 (forget(k < 0),0);
 753 0$
 754
 755 (remvalue(a,b),0);
 756 0$
 757
 758 sum(if k < 1 then a else b,k,0,2);
 759 a + 2*b$
 760
 761 block ([prederror : true], sum(if k < 1 then a else b,k,0,2));
 762 a + 2*b$
 763
 764 sum(1/3^k,k,1,inf),simpsum;
 765 1/2$
 766
 767 sum(3^-i,i,1,inf),simpsum;
 768 1/2$
 769
 770 (kill(x), assume(abs(x)-1 < 0));
 771 [abs(x) < 1];
 772
 773 sum(x^k,k,0,inf),simpsum;
 774 1/(1-x)$
 775
 776 (forget(abs(x) - 1 < 0),0);
 777 0$
 778
 779 /*--- A few challenging sums ----*/
 780
 781 (assume (n >= 1), sum (if k <= n then a else b, k, 1, n));
 782 a * n$
 783
 784 (declare(nn,integer), sum(sin(%pi * k),k,0,nn));
 785 0$
 786
 787 sum(sin(x),x,1/2,1/2 + nn);
 788 'sum(sin(x),x,1/2,1/2 + nn);
 789
 790 remove (nn, integer);
 791 done;
 792
 793 (assume_pos: assume_pos_save, apply (forget, facts ()), 0);
 794 0;
 795
 796 nusum ((n - 1)*(n - 1)!, n, 1, n);
 797 n! - 1;
 798
 799 nusum (n!, n, 1, n);
 800 'sum (n!, n, 1, n);
 801
 802 nusum ((2*m)!/(m!*(m + 1)!), m, 0, n);
 803 'sum ((2*m)!/(m!*(m + 1)!), m, 0, n);
 804
 805
 806 /*--- A test for sumcontract ---*/
 807
 808 sumcontract(sum(k,k,1,n) + sum(k^2, k, 1, n));
 809 sum(k^2+k,k,1,n);
 810
 811 /* Bug ID: 2770575 - rtestsum test 226
 812  * This test simplifies correctly after revision 1.12 of sumcon.lisp.
 813  */
 814 sumcontract(sum(k, k, 1, n) + sum(k, k, 1, n+1));
 815 n+2*sum(k,k,1,n)+1;
 816
 817 /* Tests for bug report [ 1497706 ] sum(1/k^2,k,2,inf), simpsum; */
 818
 819 /* sum (1/k^2, ...) not simplified unless specifically requested */
 820 (sum (1/k^2, k, 1, inf), [op (%%), args (%%)]);
 821 [''(nounify (sum)), [1/k^2, k, 1, inf]];
 822
 823 sum (1/k^2, k, 1, inf), simpsum;
 824 %pi^2/6;
 825
 826 /* this is the case cited in the bug report */
 827 sum (1/k^2, k, 2, inf), simpsum;
 828 %pi^2/6 - 1;
 829
 830 sum (1/k^2, k, 1, 9);
 831 9778141/6350400;
 832
 833 sum (1/k^2, k, 10, inf), simpsum;
 834 %pi^2/6 - 9778141/6350400;
 835
 836 sum (1/k^4, k, 1, 5);
 837 14001361/12960000;
 838
 839 sum (1/k^4, k, 6, inf), simpsum;
 840 %pi^4/90 - 14001361/12960000;
 841
 842 sum (1/k^7, k, 1, inf), simpsum;
 843 zeta (7);
 844
 845 sum (1/k^7, k, 1, 3);
 846 282251/279936;
 847
 848 sum (1/k^7, k, 4, inf), simpsum;
 849 zeta (7) - 282251/279936;
 850
 851 /* Tests for bug report [ 1550985 ] niceindices
 852  * Test the cases mentioned in the bug report & some others for good measure
 853  */
 854
 855 (reset (niceindicespref), 0);
 856 0;
 857
 858 niceindices (sum (1/kk, kk, 1, n) - sum (1/ii, ii, 1, n));
 859 0;
 860
 861 niceindices (sum (1/k, k, 1, n) - sum(1/i, i, 1, n));
 862 0;
 863
 864 niceindices (sum (F(foo), foo, 1, inf));
 865 'sum (F(i), i, 1, inf);
 866
 867 (niceindicespref : '[foo, bar, baz, quux], 0);
 868 0;
 869
 870 niceindices (sum (1/kk, kk, 1, n) - sum (1/ii, ii, 1, n));
 871 0;
 872
 873 niceindices (sum (1/k, k, 1, n) - sum(1/i, i, 1, n));
 874 0;
 875
 876 niceindices (sum (F(foo), foo, 1, inf));
 877 'sum (F(foo), foo, 1, inf);
 878
 879 niceindices (product (sum (product (sum (product (m*F(i) + cos(G(j))^k/sin(l), m, 1, n5), l, 1, n4), k, 1, n3), j, 1, n2), i, 1, n1));
 880 'product ('sum ('product ('sum ('product (foo*F(foo0) + cos(G(quux))^baz/sin(bar), foo, 1, n5), bar, 1, n4), baz, 1, n3), quux, 1, n2), foo0, 1, n1);
 881
 882 niceindices (sum (sum (foo/foo0, foo, 1, n), foo0, 1, m));
 883 ('sum (1/foo, foo, 1, m)) * 'sum (foo, foo, 1, n);
 884
 885 niceindices (sum (sum (1/(foo + foo0), foo, 1, n), foo0, 1, m));
 886 'sum ('sum (1/(foo + bar), foo, 1, n), bar, 1, m);
 887
 888 (reset (niceindicespref), 0);
 889 0;
 890
 891 /* Test assignments to gensumnum and genindex */
 892 errcatch (gensumnum:-5);
 893 [];
 894 errcatch(genindex:4);
 895 [];
 896 gensumnum:false;
 897 false;
 898 genindex:'zot;
 899 'zot;
 900
 901 (reset(gensumnum), gensumnum);
 902 0;
 903 (reset(genindex), genindex);
 904 'i;
 905
 906 /* Tests for bashindices, which used to not transform all indices in nested
 907  * sums and products (only the outermost)
 908  */
 909
 910 (reset (gensumnum), 0);
 911 0;
 912
 913 bashindices (sum (product (sum (F (l, m, n), l, 0, m), m, 0, n), n, 0, N));
 914 'sum ('product ('sum (F (j1, j2, j3), j1, 0, j2), j2, 0, j3), j3, 0, N);
 915
 916 (reset (gensumnum), 0);
 917 0;
 918
 919 /* Test for bug report [ 1552710 ] product(sum(f(i),i,1,inf),j,1,inf) => inf (wrong) */
 920
 921 /* Since sum(f(i), i, 1, inf) is free of j, another result
 922  * such as limit(sum(f(i), i, 1, inf)^n, n, inf) would be OK, too.
 923  * In any event, inf (observed in the bug report) is not OK.
 924  */
 925 (kill(f), block ([x : product (sum (f(i), i, 1, inf), j, 1, inf)], block ([simp : false], is (x = 'product ('sum (f(i), i, 1, inf), j, 1, inf)))));
 926
 927 true;
 928
 929 /* Test for bug reported by Amit Aronovitch to mailing list 2006/01/06
 930  */
 931
 932 (kill(f), f(x) := sum (x^n/n!, n, 0, inf), f(2*x));
 933 'sum (2^n*x^n/n!, n, 0, inf);
 934
 935
 936 /* Test for intosum */
 937
 938 intosum(2*sum(k,k,1,inf));
 939 sum(2*k,k,1,inf);
 940
 941 /* correct a bug in SUBST-IF-NOT-FREEOF (src/asum.lisp):
 942  * array flag in CAR was lost before, but it's OK now
 943  */
 944 (kill (a, i), x : a[i], sum (x, i, 1, 3));
 945 a[1] + a[2] + a[3];
 946
 947 /* Restore treatment of definite summations as it was before summation
 948  * revisions of December 2005: simplify definite summation to explicit
 949  * sum only if simpsum is true.
 950  */
 951
 952 (kill (f), sum (f(i), i, 1, 3));
 953 f(1) + f(2) + f(3);
 954
 955 'sum (f(i), i, 1, 3);
 956 'sum (f(i), i, 1, 3);
 957
 958 'sum (f(i), i, 1, 3), simpsum;
 959 f(1) + f(2) + f(3);
 960
 961 /* (analogous tests for product) */
 962
 963 product (f(i), i, 1, 3);
 964 f(1)*f(2)*f(3);
 965
 966 'product (f(i), i, 1, 3);
 967 'product (f(i), i, 1, 3);
 968
 969 'product (f(i), i, 1, 3), simpproduct;
 970 f(1)*f(2)*f(3);
 971
 972 /* Differentiate wrt a subscripted symbol in a summation.
 973  * Should work the same as for an ordinary symbol.
 974  */
 975 (kill (x, a, a1), 0);
 976 0;
 977
 978 diff ('sum (sin (a1 * x[i]), i, 1, n), a1);
 979 'sum (x[i] * cos (a1 * x[i]), i, 1, n);
 980
 981 diff ('sum (x[i]^a1, i, 1, n), a1);
 982 'sum (x[i]^a1 * log(x[i]), i, 1, n);
 983
 984 diff ('sum (sin (a[1] * x[i]), i, 1, n), a[1]);
 985 'sum (x[i] * cos (a[1] * x[i]), i, 1, n);
 986
 987 diff ('sum (x[i]^a[1], i, 1, n), a[1]);
 988 'sum (x[i]^a[1] * log(x[i]), i, 1, n);
 989
 990 /* Feature request 1848704: no way to convert log(product(...)) to sum(log(...)...)
 991  */
 992
 993 (kill (u, a, x, T), 0);
 994 0;
 995
 996 log (product (u[k], k, 1, n));
 997 log (product (u[k], k, 1, n));
 998
 999 log (product (u[k], k, 1, n)), logexpand=true;
1000 log (product (u[k], k, 1, n));
1001
1002 log (product (u[k], k, 1, n)), logexpand=all;sum (log (u[k]), k, 1, n);
1003
1004 log (product (u[k], k, 1, n)), logexpand=super;
1005 sum (log (u[k]), k, 1, n);
1006
1007 log (product (a*%e^x(i), i, 1, T));
1008 log (product (a*%e^x(i), i, 1, T));
1009
1010 log (product (a*%e^x(i), i, 1, T)), logexpand=all;
1011 /* Yields sum(const + ...) instead of const*T + sum(...)
1012  * since sum is not known to be linear by default.
1013  * Change this result if default ever changes.
1014  */
1015 sum (log(a) + x(i), i, 1, T);
1016
1017 /* bug reported to mailing list 2008-04-17:
1018  * simpsum causes two evaluations of summand, should be only one
1019  */
1020 (bar (x, y) := sum (F (i, x, y), i, 1, n),
1021  ev (bar (y, x), simpsum=false));
1022 'sum (F (i, y, x), i, 1, n);
1023
1024 ev (bar (y, x), simpsum=true);
1025 'sum (F (i, y, x), i, 1, n);
1026
1027 /* original problem */
1028
1029 (kill (all),
1030  geom(x,m,n):= (x^m-x^(n+1))/(1-x),
1031  B(s,sp):=-phi[1](s)*phi[1](s+sp)*phi[0](s+sp)/(1-phi[1](s+sp)*phi[0](s +sp))
1032  +(1-phi[1](s))*sum((phi[1](s)*phi[0](s))^n*(
1033      geom( (phi[0](s+sp)*phi[1](s+sp))/(phi[0](s)*phi[1](s)), 1, n-1 )
1034      -phi[1](s+sp)/phi[1](s)
1035      *geom((phi[0](s+sp)*phi[1](s+sp))/(phi[0](s)*phi[1](s)), 0, n-1))
1036    , n , 1, inf ),
1037  A(s):=phi[1](s)*(1+phi[0](s))/(1-phi[1](s)*phi[0](s)),
1038  C(s,sp):=A(s+sp)+B(s,sp)+B(sp,s),
1039  0);
1040 0;
1041
1042 (expr_1 : ev (C (s, sp), simpsum=false),
1043  expr_2 : ev (C (s, sp), simpsum=true),
1044  is (expr_1 = expr_2));
1045 true;
1046
1047 (expr_1 : ev (C (sp, s), simpsum=false),
1048  expr_2 : ev (C (sp, s), simpsum=true),
1049  is (expr_1 = expr_2));
1050 true;
1051
1052 is (sublis ([s=sp, sp=s], C (s, sp)) = C (s, sp)), simpsum;
1053 true;
1054
1055 is (sublis ([s=sp, sp=s], C (sp, s)) = C (sp, s)), simpsum;
1056 true;
1057
1058 is (C (sp, s) = C (s, sp)), simpsum;
1059 true;
1060
1061 /* Tests for fbino in combin.lisp */
1062 sum(binomial(2*n-k,k), k, 0, n), simpsum;
1063 fib(2*n+1);
1064
1065 sum(binomial(n-k,k), k, 1, n), simpsum;
1066 fib(n+1)-1;
1067
1068 (assume(n>1), 0);
1069 0;
1070
1071 sum(binomial(n-k,k+1), k, 1, n-1), simpsum;
1072 fib(n+2)-n-1;
1073
1074 sum(binomial(2*n,2*k), k, 0, n), simpsum;
1075 2^(2*n-1);
1076
1077 sum(binomial(n,2*k+1), k, 0, n), simpsum;
1078 2^(n-1);
1079
1080 sum(binomial(n,2*k), k, 0, n), simpsum;
1081 2^(n-1);
1082
1083 sum(binomial(n,2*k), k, 1, n), simpsum;
1084 2^(n-1)-1;
1085
1086 sum(binomial(n+1,k), k, 1, n), simpsum;
1087 2^(n+1)-2;
1088
1089 sum(binomial(n,k), k, -1, n+1), simpsum;
1090 2^n;
1091
1092 sum(binomial(t+i,t),i,0,k), simpsum;
1093 binomial(t+k+1,t+1);
1094
1095 /* Because of revision 1.43 of csimp2.lisp binomial(x,x) no longer simplifies
1096  * automatically to 1. The sign must be known not be negative.
1097  * Therefore, we add a fact about the sign of t.
1098  */
1099 (assume(t>0),done);
1100 done;
1101 sum(binomial(t+i,t),i,1,k), simpsum;
1102 binomial(t+k+1,t+1)-1;
1103 (forget(t>0),done);
1104 done;
1105
1106 sum(binomial(t+i,i),i,0,k), simpsum;
1107 binomial(t+k+1,t+1);
1108
1109 (forget(n>1), 0);
1110 0;
1111
1112 /* bug in assume database causes trouble for conditional in summation
1113  * see mailing list 2009-07-29 "Fourier Series using fourie.mac"
1114  */
1115
1116 /* If ever the assume stuff is revised, this next test can go away.
1117  * It's here to ensure the underlying problem is fixed.
1118  */
1119 ?dcomp (1, 1.0);
1120 zero;
1121
1122 (kill (all), assume (xx >= 1), 0);
1123 0;
1124
1125 is (xx > 1);
1126 unknown;
1127
1128 if xx > 1 then aa else bb;
1129 if xx > 1 then aa else bb;
1130
1131 sum (if xx > 1 then aa else bb, xx, 1, 3);
1132 2*aa + bb;
1133
1134 foo : 'sum (if xx > 1 then aa else bb, xx, 1, 3);
1135 'sum (if xx > 1 then aa else bb, xx, 1, 3);
1136
1137 ev (foo, nouns);
1138 2*aa + bb;
1139
1140 foo : sum (if equal (xx, 1) then aa else bb, xx, 1, nn);
1141 'sum (if equal (xx, 1) then aa else bb, xx, 1, nn);
1142
1143 ev (foo, nouns, nn=3);
1144 aa + 2*bb;
1145
1146 (reset(assume_pos),1);
1147 1;
1148
1149 /* SF bug #3605: "Variable confusion in function handling Taylor series" */
1150
1151 kill(x, y);
1152 done;
1153
1154 /* function definition copied verbatim from bug report */
1155
1156 (test(f) := block([i:taylorinfo(f)[1]], sum(coeff(f, i[1], x) * i[1]^x, x, 0, i[3])), 0);
1157 0;
1158
1159 ratdisrep (test(taylor(cos(y), y, 0, 4)));
1160 (y^4-12*y^2+24)/24;
1161
1162 ratdisrep (test(taylor(cos(x), x, 0, 4)));
1163 29/3;
1164
1165 /* declare summation index with different name in block to distinguish from global */
1166 (test1(f) := block([i:taylorinfo(f)[1], x%], sum(coeff(f, i[1], x%) * i[1]^x%, x%, 0, i[3])), 0);
1167 0;
1168
1169 ratdisrep (test1(taylor(cos(y), y, 0, 4)));
1170 (y^4-12*y^2+24)/24;
1171
1172 ratdisrep (test1(taylor(cos(x), x, 0, 4)));
1173 (x^4-12*x^2+24)/24;
1174
1175 /* Bug #2267: simpsum error */
1176
1177 sum((2^(k+i))/(2^(k+i)+2^(i)),i,0,k), simpsum;
1178 (2^k*(1 + k))/(1 + 2^k);
1179