symbolic summation interface
Burcin Erocal
I plan to start implementing a user interface to symbolic summation
soon and I want to get some opinions on how this interface should be.
The most natural construct for summation, either of a list or symbolic
summation is of course "sum." Initially, I was thinking that it's a big
sin to override the python "sum" function, since I also use it very
often and consider it (performance) critical. However, I was encouraged
to think again about this after Mike's comment here:
http://trac.sagemath.org/sage_trac/ticket/3587#comment:3
For quick reference, here's the help for Python's sum function:
sum(sequence, start=0) -> value
Returns the sum of a sequence of numbers (NOT strings) plus the
value of parameter 'start'. When the sequence is empty, returns start.
We could easily extend this to check if the second parameter is a tuple
which defines a range (e.g., (x, 1, n) where the upper and lower
bounds are inclusive). If it is, try to solve the sum symbolically,
otherwise call the python sum function. (Actually, I recall that there
were plans to overwrite this function anyway with one that does
balanced summation if the argument is a list.)
So I propose the following:
sage: var('i,n')
(i, n)
sage: sum(2^i, (i, 0, n))
2^(n+1) - 1
sage: sum(1/i, (i, 1, n))
harmonic_number(1, n)
(As far as I can see, Sage doesn't have a construct for harmonic
numbers yet.)
These still will work of course:
sage: sum(range(5))
10
sage: sum(i for i in range(5))
10
sage: sum(range(5), 2)
12
Similarly, I suggest we extend prod the same way.
Any comments or objections?
I should also add that I don't like to differentiate between upper and
lowercase commands, so I don't think using "Sum" is an option. If
overriding "sum" is not accepted, "SymbolicSum" might be an alternative.
Thanks.
Burcin
William Stein
+1
Note that Robert Bradshaw wrote the current prod (in Cython), and it is
balanced and very fast.
> I should also add that I don't like to differentiate between upper and
> lowercase commands, so I don't think using "Sum" is an option.
crap.
> If overriding "sum" is not accepted, "SymbolicSum" might be an alternative.
>
>
> Thanks.
>
> Burcin
>
> >
>
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
Robert Bradshaw
This certainly seems like the best notation. Sums with tuple
endpoints are sometimes useful, e.g.
sage: sum(zip(range(5), range(5)), ())
(0, 0, 1, 1, 2, 2, 3, 3, 4, 4)
However, this is a much more obscure use case than what you're
proposing. We could make our sum should take an optional keyword
parameter "endpoint" so that this behavior is still possible (even
though we're changing the default behavior).
>> I should also add that I don't like to differentiate between upper
>> and
>> lowercase commands, so I don't think using "Sum" is an option.
>
> I agree. I'm strongly opposed to that Maple-style case-has-meaning
> crap.
Me too.
- Robert
Jason Grout
>
> (Actually, I recall that there
> were plans to overwrite this function anyway with one that does
> balanced summation if the argument is a list.)
FYI, see trac #2737
Jason
Jason Grout
> otherwise call the python sum function. (Actually, I recall that there
> were plans to overwrite this function anyway with one that does
> balanced summation if the argument is a list.)
>
See trac 2737
Jason
Franco Saliola
What about list comprehensions? Something like the following.
sage: var('i,n')
(i, n)
sage: sum(2^i for i in range(n+1))
2^(n+1) - 1
I ask because this seems like the natural/first thing a user would try.
Franco
--
Jason Grout
>>
>> Any comments or objections?
>
> What about list comprehensions? Something like the following.
>
> sage: var('i,n')
> (i, n)
> sage: sum(2^i for i in range(n+1))
> 2^(n+1) - 1
>
> I ask because this seems like the natural/first thing a user would try.
It's consistent:
sage: sum(i for i in range(10))
45
However, range would have to also be overridden here, since range only
accepts integers.
Or we could override the ellipsis notation:
sum(i for i in (1..oo))
sum(i for i in (1..n+1))
hmmm, interesting ideas.
Jason
Mike Hansen
they create a "black box" generator which gets passed into sum.
--Mike
Harald Schilly
> > sage: sum(2^i for i in range(n+1))
than one variable and relationships between them (e.g. for all i,j in
N where i>j and i mod 7 == 0 and so on) Therefore, I think a new kind
of object to represent the set over which the sum is made must be
created. This set could be created by various constructors and returns
"information" about itself (not a black box).
A = sumset(vars=[i,j], start=[0,0], constraint={ i>j , mod(i,7)==0},
intersect=sumset_B, union=sumset_C)
and then things like: values = A.getGenerator() ; A.isInfinite(i) ->
true ; j_vals_sumset = A.getSet(i=fixed_value)
But that's just a small idea, i don't know how a good interface would
look like and how it is done anywhere else (I assume all other
implmenetations directly connect the input of a sum-function with the
code to do the summation and don't use any interfaces at all)
Also, the expression in the sum needs to be transformed into certain
standard forms if possible.
h
Jason Grout
I agree that they probably aren't feasible. I was throwing them out for
the purposes of brainstorming for interfaces.
Jason
Ronan Paixão
> On Nov 14, 9:26 am, Jason Grout <jason-s...@creativetrax.com> wrote:
> > > sage: sum(2^i for i in range(n+1))
>
> > sum(i for i in (1..oo))
>
> They look nice, but i doubt if they could work. e.g. think of more
> than one variable and relationships between them (e.g. for all i,j in
> N where i>j and i mod 7 == 0 and so on) Therefore, I think a new kind
> of object to represent the set over which the sum is made must be
> created. This set could be created by various constructors and returns
> "information" about itself (not a black box).
> A = sumset(vars=[i,j], start=[0,0], constraint={ i>j , mod(i,7)==0},
> intersect=sumset_B, union=sumset_C)
> and then things like: values = A.getGenerator() ; A.isInfinite(i) ->
> true ; j_vals_sumset = A.getSet(i=fixed_value)
No need to complicate with another object.
Something like that would be nicer (if possible):
sum(sum([f(j%7) for j in (0..oo) if j>k]) for k in (0..oo)]))
using j and k, since i is parsed as the complex constant.
Ronan Paixão
Robert Bradshaw
Given that they're by far the most natural interface, we shouldn't
rule them out for sure. Consider
sage: import dis
sage: g = (i for i in range(10) if i % 3r)
sage: dis.dis(foo(10).gi_frame.f_code)
2 0 SETUP_LOOP 29 (to 32)
3 LOAD_GLOBAL 0 (range)
6 LOAD_FAST 0 (n)
9 CALL_FUNCTION 1
12 GET_ITER
>> 13 FOR_ITER 15 (to 31)
16 STORE_FAST 1 (i)
3 19 LOAD_FAST 1 (i)
22 LOAD_FAST 1 (i)
25 BINARY_MULTIPLY
26 YIELD_VALUE
27 POP_TOP
28 JUMP_ABSOLUTE 13
>> 31 POP_BLOCK
>> 32 LOAD_CONST 0 (None)
35 RETURN_VALUE
If we evaluate the conditional (12-19) with the loop variable set to
a symbolic we get the condition, and we evaluate the expression
(22-25) at a symbolic loop variable to get the generic term. The
inner iterator is local, and we can recognize range, xrange, and
recursively attempt to descend for others. Of course the [a..b]
notation would have to be extended to support symbolic endpoints, and
one would have to be careful to not botch up arbitrary generators.
- Robert
Robert Bradshaw
> On Nov 14, 2008, at 7:05 AM, Jason Grout wrote:
>
Check out the (still not 100% robust) code below.
sage: get_sum_components(n for n in range(10))
(n, None, <listiterator object at 0xd1b0cd0>)
sage: get_sum_components(n^2 for n in [1..10] if n != 3)
(n^2, n != 3, <listiterator object at 0xd1b0ff0>)
sage: get_sum_components(1/n^2 for n in (1..))
(1/n^2, None, <generator object at 0xd1abe18>)
sage: get_sum_components(n^2 for n in [1..100] if n != 10)
(n^2, n != 10, <listiterator object at 0xd1afef0>)
---------------------------------------------
import dis
def chop_code(c):
i = 0
all = []
while i < len(c):
if ord(c[i]) < dis.HAVE_ARGUMENT:
all.append(c[i])
else:
all.append(c[i:i+3])
i += len(all[-1])
return all
def get_sum_components(g):
frame = g.gi_frame
code = frame.f_code
ops = chop_code(code.co_code)
if (dis.opname[ord(ops[0][0])] != 'SETUP_LOOP'
or dis.opname[ord(ops[1][0])] != 'LOAD_FAST'
or dis.opname[ord(ops[2][0])] != 'FOR_ITER'):
raise TypeError, "Not a valid generator."
index_name = code.co_varnames[ord(ops[3][1])]
range_name = code.co_varnames[ord(ops[1][1])]
range_value = frame.f_locals[range_name]
cond = []
expr = []
for op in ops[4:]:
if dis.opname[ord(op[0])] == 'YIELD_VALUE':
break
elif dis.opname[ord(op[0])] == 'JUMP_IF_FALSE' and not cond:
cond = expr
expr = []
else:
expr.append(op)
if cond:
# there's an extra pop
expr = expr[1:]
index_var = var(index_name)
expr += [chr(dis.opmap['RETURN_VALUE'])]
def func(arg): pass
expr_code = type(code)(2, code.co_nlocals, code.co_stacksize,
32832, ''.join(expr),
code.co_consts, code.co_names,
code.co_varnames, code.co_filename, 'sum_expr', code.co_firstlineno,
'', code.co_freevars, code.co_cellvars)
func.func_code = expr_code
expr = func(None, index_var)
if cond:
cond += [chr(dis.opmap['RETURN_VALUE'])]
cond_code = type(code)(2, code.co_nlocals,
code.co_stacksize, 32832, ''.join(cond),
code.co_consts, code.co_names,
code.co_varnames, code.co_filename, 'sum_cond', code.co_firstlineno,
'', code.co_freevars, code.co_cellvars)
func.func_code = cond_code
cond = func(None, index_var)
else:
cond = None
return expr, cond, range_value