tests/wester_problems/test_sums.mac

   1 /* Original version of this file copyright 1999 by Michael Wester,
   2  * and retrieved from http://www.math.unm.edu/~wester/demos/Sums/problems.macsyma
   3  * circa 2006-10-23.
   4  *
   5  * Released under the terms of the GNU General Public License, version 2,
   6  * per message dated 2007-06-03 from Michael Wester to Robert Dodier
   7  * (contained in the file wester-gpl-permission-message.txt).
   8  *
   9  * See: "A Critique of the Mathematical Abilities of CA Systems"
  10  * by Michael Wester, pp 25--60 in
  11  * "Computer Algebra Systems: A Practical Guide", edited by Michael J. Wester
  12  * and published by John Wiley and Sons, Chichester, United Kingdom, 1999.
  13  */
  14 /* ----------[ M a c s y m a ]---------- */
  15 /* ---------- Initialization ---------- */
  16 showtime: all$
  17 prederror: false$
  18 /* ---------- Sums ---------- */
  19 load(nusum1)$
  20 /* Simplify the sum below to sum(x[i]^2, i = 1..n) - sum(x[i], i = 1..n)^2/n */
  21 xbar: sum(x[j], j, 1, n) / n;
  22 sum((x[i] - xbar)^2, i, 1, n);
  23 niceindices(sumexpand(expand(%)));
  24 remvalue(xbar)$
  25 /* Derivation of the least squares fitting of data points (x[i], y[i]) to a
  26    line y = m x + b.  See G. Keady, ``Using Maple's linalg package with Zill
  27    and Cullen _Advanced Engineering Mathematics_, Part II: Vectors, Matrices
  28    and Vector Calculus'', University of Western Australia,
  29    ftp://maths.uwa.edu.au/pub/keady/ */
  30 f: sum((y[i] - m*x[i] - b)^2, i, 1, n)$
  31 solve([diff(f, m) = 0, diff(f, b) = 0], [m, b]);
  32 remvalue(f)$
  33 load("nusum1")$
  34 /* Indefinite sum => (-1)^n binomial(2 n, n).  See Herbert S, Wilf,
  35    ``IDENTITIES and their computer proofs'', University of Pennsylvania. */
  36 closedform(indefsum((-1)^k * binomial(2*n, k)^2, k));
  37 /* Check whether the full Gosper algorithm is implemented
  38    => 1/2^(n + 1) binomial(n, k - 1) */
  39 closedform(indefsum(binomial(n, k)/2^n - binomial(n + 1, k)/2^(n + 1), k));
  40 factcomb(makefact(%));
  41 /* Dixon's identity (check whether Zeilberger's algorithm is implemented).
  42    Note that the indefinite sum is equivalent to the definite
  43    sum(..., k = -min(a, b, c)..min(a, b, c)) => (a + b + c)!/(a! b! c!)
  44    [Wilf] */
  45 closedform(indefsum((-1)^k * binomial(a+b, a+k) * binomial(b+c, b+k)
  46                            * binomial(c+a, c+k), k));
  47 /* Telescoping sum => g(n + 1) - g(0) */
  48 closedform(sum(g(k + 1) - g(k), k, 0, n));
  49 /* => n^2 (n + 1)^2 / 4 */
  50 closedform(sum(k^3, k, 1, n));
  51 /* See Daniel I. A. Cohen, _Basic Techniques of Combinatorial Theory_, John
  52    Wiley and Sons, 1978, p. 60.  The following two sums can be derived directly
  53    from the binomial theorem:
  54    sum(k^2 * binomial(n, k) * x^k, k = 1..n) = n x (1 + n x) (1 + x)^(n - 2)
  55    => n (n + 1) 2^(n - 2)   [Cohen, p. 60] */
  56 closedform(sum(k^2 * binomial(n, k), k, 1, n));
  57 /* => [2^(n + 1) - 1]/(n + 1)   [Cohen, p. 83] */
  58 closedform(sum(binomial(n, k)/(k + 1), k, 0, n));
  59 /* Vandermonde's identity => binomial(n + m, r)   [Cohen, p. 31] */
  60 closedform(sum(binomial(n, k) * binomial(m, r - k), k, 0, r));
  61 /* => Fibonacci[2 n]   [Cohen, p. 88] */
  62 closedform(sum(binomial(n, k) * fib(k), k, 0, n));
  63 /* => Fibonacci[n] Fibonacci[n + 1]   [Cohen, p. 65] */
  64 closedform(sum(fib(k)^2, k, 1, n));
  65 factor(ratsimp(%));
  66 ratsimp(% - fibtophi(fib(n) * fib(n + 1)));
  67 /* => 1/2 cot(x/2) - cos([2 n + 1] x/2)/[2 sin(x/2)]
  68    See Konrad Knopp, _Theory and Application of Infinite Series_, Dover
  69    Publications, Inc., 1990, p. 480. */
  70 closedform(sum(sin(k*x), k, 1, n));
  71 factor(trigreduce(rectform(%)));
  72 halfangles(% - 1/2*cot(x/2) - cos((2*n + 1)*x/2)/(2*sin(x/2)));
  73 /* => sin(n x)^2/sin x   [Gradshteyn and Ryzhik 1.342(3)] */
  74 closedform(sum(sin((2*k - 1)*x), k, 1, n));
  75 /* => Fibonacci[n + 1]   [Cohen, p. 87] */
  76 closedform(sum(binomial(n - k, k), k, 0, floor(n/2)));
  77 /* => pi^2 / 6 + zeta(3) =~ 2.84699 */
  78 closedform(sum(1/k^2 + 1/k^3, k, 1, inf));
  79 /* sfloat(%); */
  80 bfloat(%);
  81 /* => pi^2/12 - 1/2 (log 2)^2   [Gradshteyn and Ryzhik 0.241(2)] */
  82 closedform(sum(1/(2^k*k^2), k, 1, inf));
  83 /* => pi/12 sqrt(3) - 1/4 log 3   [Knopp, p. 268] */
  84 closedform(sum(1/((3*k + 1)*(3*k + 2)*(3*k + 3)), k, 0, inf));
  85 /* => 1/2 (2^(n - 1) + 2^(n/2) cos(n pi/4))   [Gradshteyn and Ryzhik 0.153(1)]
  86    */
  87 closedform(sum(binomial(n, 4*k), k, 0, inf));
  88 /* => 1   [Knopp, p. 233] */
  89 closedform(sum(1/(sqrt(k*(k + 1)) * (sqrt(k) + sqrt(k + 1))), k, 1, inf));
  90 /* => 1/sqrt([1 - x y]^2 - 4 x^2)   (| x y | < 1 and -1 <= x < 1).
  91       From Evangelos A. Coutsias, Michael J. Wester and Alan S. Perelson, ``A
  92       Nucleation Theory of Cell Surface Capping'', draft. */
  93 'sum('sum(binomial(n, k)*binomial(n - k, n - 2*k)*x^n*y^(n - 2*k),
  94           k, 0, floor(n/2)),
  95      n, 0, inf);
  96 /* An equivalent summation to the above is: */
  97 'sum('sum(n!/(k!^2*(n - 2*k)!)*(x/y)^k*(x*y)^(n - k), n, 2*k, inf), k, 0, inf);
  98 closedform(%);
  99 /* => pi/2   [Knopp, p. 269] */
 100 'sum('product(k/(2*k - 1), k, 1, m), m, 2, inf);
 101 closedform(%);