# A simple sum causes division by zero exception

I'm very new to SageMath and did some small experiments yesterday. This is my Code:

x, n, r, i = var('x, n, r, i')
f(x, n, r) = x * sum((1 + r) ^ i, i, 1, n)
f(1, 30, 0.0)


It works well if r != 0, but will raise a exception of "division by zero" if equal.

My ipynb file: nbviewer.jupyter.org/gist/7sDream/5db3cfd153269fd1a1cacaf0b60f69bb

It seems some optimizations of sum did not consider the range of variable r. But there is also no change if I added assume(r >= 0)before define f.

So how can I disable this sum optimizations? Or what is wrong with the way I use it?

edit retag close merge delete

Curiously :

sage: def ef(x, n, r): return sum([(1+r)^u for u in (1..n)])
sage: ef(1, 30, 0)
30
sage: ef(1, 30, 0.0)
30.0000000000000


Interestingly, if we define fc = fast_callable( x * sum((1 + r) ^ i, i, 1, n), vars=[x,n,r] ), then fc(1, 30, 0.0) evaluates to NaN.

1

Note that the "optimization" takes place at the creation of f :

sage: var("x, n, r, i")
(x, n, r, i)
sage: f(x,n,r,i)=x*sum((1+r)^i,i,1,n)
sage: f
(x, n, r, i) |--> ((r + 1)^(n + 1) - r - 1)*x/r


FWIW, sympy has an interesting translation :

>>> x, n, r, i = symbols("x, n, r, i")
>>> f = Lambda((x, n, r), x*summation((1+r)**i, (i, 1, n)))
>>> f
Lambda((x, n, r), x*Piecewise((n, Eq(r + 1, 1)), (-(r - (r + 1)**(n + 1) + 1)/r, True)))


Sort by » oldest newest most voted

Sage's implementation of sum is incorrect in this case :

sage: sum((1+r)^i, i, 1, n)
((r + 1)^(n + 1) - r - 1)/r


which is incorrect if r==0.

Compare this with sympy's implementation of summation :

>>> x, n, r, i = symbols("x, n, r, i")
>>> x*summation((1+r)**i, (i, 1, n))
x*Piecewise((n, Eq(r + 1, 1)), (-(r - (r + 1)**(n + 1) + 1)/r, True))


Interestingly, Mathematica falls in the same trap as sage :

In:= InputForm[Sum[(1+r)^i, {i, 1, n}]]
Out//InputForm= ((1 + r)*(-1 + (1 + r)^n))/r


The workaround is to use the algorithm="sympy" optional argument of sum :

sage: sum((1+r)^i, i, 1, n, algorithm="sympy")
cases(((r + 1 == 1, n), (1, ((r + 1)^(n + 1) - r - 1)/r)))
sage: f(x,n,r,i)=x*sum((1+r)^i,i,1,n, algorithm="sympy")
sage: f
(x, n, r, i) |--> x*cases(((r + 1 == 1, n), (1, ((r + 1)^(n + 1) - r - 1)/r)))
sage: f(1, 30, 0)
30
sage: f(1, 30, 0.0)
30


HTH,

EDIT : This seems to be previously unreported, and is now Trac#31418.

more