1%% Copyright (C) 2014-2016, 2019 Colin B. Macdonald
2%%
3%% This file is part of OctSymPy.
4%%
5%% OctSymPy is free software; you can redistribute it and/or modify
6%% it under the terms of the GNU General Public License as published
7%% by the Free Software Foundation; either version 3 of the License,
8%% or (at your option) any later version.
9%%
10%% This software is distributed in the hope that it will be useful,
11%% but WITHOUT ANY WARRANTY; without even the implied warranty
12%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
13%% the GNU General Public License for more details.
14%%
15%% You should have received a copy of the GNU General Public
16%% License along with this software; see the file COPYING.
17%% If not, see <http://www.gnu.org/licenses/>.
18
19%% -*- texinfo -*-
20%% @documentencoding UTF-8
21%% @defmethod  @@sym symsum (@var{f}, @var{n}, @var{a}, @var{b})
22%% @defmethodx @@sym symsum (@var{f}, @var{n}, [@var{a} @var{b}])
23%% @defmethodx @@sym symsum (@var{f}, @var{a}, @var{b})
24%% @defmethodx @@sym symsum (@var{f}, [@var{a} @var{b}])
25%% @defmethodx @@sym symsum (@var{f}, @var{n})
26%% @defmethodx @@sym symsum (@var{f})
27%% Symbolic summation.
28%%
29%% The sum of the expression @var{f} for @var{n} from @var{a} to
30%% @var{b}.  When @var{n} is omitted it is determined using
31%% @code{symvar} and defaults to @code{x} if @var{f} is constant. The
32%% limits @var{a} and @var{b} default to 0 and @var{n} - 1
33%% respectively.
34%%
35%% @example
36%% @group
37%% syms n m
38%% symsum(1/n^2, n, 1, m)
39%%   @result{} (sym) harmonic(m, 2)
40%%
41%% symsum(exp(2*n)/sin(n), n, 2*m, 6*m)
42%%   @result{} (sym)
43%%           6⋅m
44%%           ____
45%%           ╲
46%%            ╲      2⋅n
47%%             ╲    ℯ
48%%             ╱   ──────
49%%            ╱    sin(n)
50%%           ╱
51%%           ‾‾‾‾
52%%         n = 2⋅m
53%% @end group
54%% @end example
55%%
56%% @seealso{@@sym/symprod, @@sym/sum}
57%% @end defmethod
58
59
60function S = symsum(f, n, a, b)
61
62  if (nargin > 4)
63    print_usage ();
64  end
65
66  idx1.type = '()';
67  idx1.subs = {1};
68  idx2.type = '()';
69  idx2.subs = {2};
70
71  if (nargin == 1)
72    n = symvar(f, 1);
73    if (isempty(n))
74      n = sym('x');
75    end
76    a = sym(0);
77    b = n - 1;
78  elseif (nargin == 2) && (length(n) == 2)
79    f = sym(f);
80    %a = n(1);  % issue #17
81    %b = n(2);
82    a = subsref(n, idx1);
83    b = subsref(n, idx2);
84    n = symvar(f, 1);
85    if (isempty(n))
86      n = sym('x');
87    end
88  elseif (nargin == 2)
89    f = sym(f);
90    n = sym(n);
91    a = sym(0);
92    b = n - 1;
93  elseif (nargin == 3) && (length(a) == 2)
94    f = sym(f);
95    n = sym(n);
96    %b = a(2);  % issue #17
97    %a = a(1);
98    b = subsref(a, idx2);
99    a = subsref(a, idx1);
100  elseif (nargin == 3)
101    f = sym(f);
102    b = a;
103    a = n;
104    n = symvar(f, 1);
105    if (isempty(n))
106      n = sym('x');
107    end
108  else
109    f = sym(f);
110    n = sym(n);
111    a = sym(a);
112    b = sym(b);
113  end
114
115  cmd = { '(f, n, a, b) = _ins'
116           'S = sp.summation(f, (n, a, b))'
117           'return S,' };
118
119  S = pycall_sympy__ (cmd, sym(f), sym(n), sym(a), sym(b));
120
121end
122
123
124%!error <Invalid> symsum (sym(1), 2, 3, 4, 5)
125
126%!test
127%! % finite sums
128%! syms n
129%! assert (isequal (symsum(n,n,1,10), 55))
130%! assert(isa(symsum(n,n,1,10), 'sym'))
131%! assert (isequal (symsum(n,n,sym(1),sym(10)), 55))
132%! assert (isequal (symsum(n,n,sym(1),sym(10)), 55))
133%! assert (isequal (symsum(1/n,n,1,10), sym(7381)/2520))
134
135%!test
136%! % negative limits
137%! syms n
138%! assert (isequal (symsum(n,n,-3,3), sym(0)))
139%! assert (isequal (symsum(n,n,-3,0), sym(-6)))
140%! assert (isequal (symsum(n,n,-3,-1), sym(-6)))
141
142%!test
143%! % one input
144%! syms n
145%! f = symsum (n);
146%! g = n^2/2 - n/2;
147%! assert (isequal (f, g))
148%! f = symsum (2*n);
149%! g = n^2 - n;
150%! assert (isequal (f, g))
151
152%!test
153%! % constant input
154%! f = symsum (sym(2));
155%! syms x
156%! g = 2*x;
157%! assert (isequal (f, g))
158
159%!test
160%! % two inputs
161%! syms n
162%! f = symsum (2*n, n);
163%! g = n^2 - n;
164%! assert (isequal (f, g))
165
166%!test
167%! % two inputs, second is range
168%! syms n
169%! f = symsum (n, [1 6]);
170%! g = 21;
171%! assert (isequal (f, g))
172%! f = symsum (n, [sym(1) 6]);
173%! g = 21;
174%! assert (isequal (f, g))
175%! f = symsum (2*n, [1 6]);
176%! g = 2*21;
177%! assert (isequal (f, g))
178
179%!test
180%! % three inputs, last is range
181%! syms n
182%! f = symsum (2*n, n, [1 4]);
183%! g = sym(20);
184%! assert (isequal (f, g))
185%! f = symsum (2*n, n, [sym(1) 4]);
186%! g = sym(20);
187%! assert (isequal (f, g))
188%! f = symsum (2, n, [sym(1) 4]);
189%! g = sym(8);
190%! assert (isequal (f, g))
191
192%!test
193%! % three inputs, no range
194%! syms n
195%! f = symsum (2*n, 1, 4);
196%! g = sym(20);
197%! assert (isequal (f, g))
198%! f = symsum (5, sym(1), 3);
199%! g = sym(15);
200%! assert (isequal (f, g))
201
202%!test
203%! % ok to use double's for arguments in infinite series
204%! syms n oo
205%! assert(isequal(symsum(1/n^2,n,1,oo), sym(pi)^2/6))
206%! assert(isequal(symsum(1/n^2,n,1,inf), sym(pi)^2/6))
207
208%!test
209%! % should be oo because 1 is real but seems to be
210%! % zoo/oo depending on sympy version
211%! syms n oo
212%! zoo = sym('zoo');
213%! assert (isequal (symsum(1/n,n,1,oo), oo) || ...
214%!         isequal (symsum(1/n,n,1,oo), zoo))
215