Ask Your Question

substitute() raises exception where python does not

asked 2016-03-09 08:23:25 -0600

stan gravatar image

updated 2016-03-09 09:32:23 -0600

In my usual workflow, I derive symbolic equations and then substitute their arguments by values stored in dictionaries to obtain numerical solutions. Now I found out that the substitute() command raises division by zero exception where standard python does not. This is a bit of a problem, as these exceptions prevent my code from evaluating all data sets once such an exception occurs and I don't want to wrap each line of code in a try-except construct. Does anyone know how the exception could be avoided? Here is an example:

var('y a x')
eq_y = y == a/x
vdict = {}
vdict[a] = 1.
vdict[x] = 0.


However, if I substitute vdict into eq_y, I get an exception:


Traceback (click to the left of this block for traceback)
ValueError: power::eval(): division by zero
edit retag flag offensive close merge delete


You forgot to declare 'y' as a variable...

vdelecroix gravatar imagevdelecroix ( 2016-03-09 09:08:14 -0600 )edit

3 answers

Sort by ยป oldest newest most voted

answered 2016-03-09 16:44:51 -0600

nbruin gravatar image

You can use fast_callable to compile a symbolic expression into something that evaluates more in the way of a python expression:

sage: fc=fast_callable(eq_y.rhs(),vars=[x,a])
sage: fc(0.0,1.0)
edit flag offensive delete link more

answered 2016-03-09 09:23:53 -0600

kcrisman gravatar image

updated 2016-03-09 09:24:25 -0600

I see. You are using a symbolic equation in the latter case, which will definitely barf at dividing by zero, but the things you are putting in are "real literals".

sage: type(a)
<type 'sage.symbolic.expression.Expression'>
sage: type(vdict[a])
<type 'sage.rings.real_mpfr.RealLiteral'>

As such, they obey 'usual' division rules for floating-point things like

sage: vdict[x]/vdict[x]

So if you don't want mathematically correct behavior, I suggest not using a symbolic variable!

sage: eq_y = y == z
sage: def mysubs(a=.1,x=0.,eq=eq_y):
....:     t = a/x
....:     return eq.subs(z=t)
sage: mysubs()
y == +infinity
sage: mysubs(.0,.0)
y == NaN

Indeed, we get

sage: vdict = {}
sage: vdict['a'] = 0.
sage: vdict['x'] = 0.
sage: mysubs(**vdict)
y == NaN
edit flag offensive delete link more


Thanks for this, but this involves re-writing the original symbolic equation inside the function, so it defeats my purpose of using a symbolic equation in the first place. Anyway, the explanations given above will help me formulate the real problem, which is how to catch such exceptions in a code while continuing with the evaluation. I'll post it shortly.

stan gravatar imagestan ( 2016-03-09 09:51:33 -0600 )edit

answered 2016-03-09 09:17:27 -0600


The problem comes from the fact that Sage floating point number behave differently than other numbers with respect to division by zero. Namely

sage: 1 / 0
Traceback (most recent call last):
ZeroDivisionError: rational division by zero
sage: 1. / 0.

I opened a discussion. But I think that the second example should raise the same error as the first one. See this sage-devel thread.

edit flag offensive delete link more


Incorrect, at least not if we want to stay sort of close to IEEE specs for floating point stuff. See e.g.

kcrisman gravatar imagekcrisman ( 2016-03-09 09:25:19 -0600 )edit

Of course, points this out: sage: float(1.0)/float(0.0) so I don't know what we should do... I'd go with IEEE on this, since I assume the MPFR people know what they are doing - ?

kcrisman gravatar imagekcrisman ( 2016-03-09 09:27:10 -0600 )edit

It is also Python choice

$ python -c "1. / 0."
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ZeroDivisionError: float division by zero
vdelecroix gravatar imagevdelecroix ( 2016-03-09 09:27:40 -0600 )edit

OK, thanks, I didn't realise that this is an inconsistency in sage (not throwing an exception for 1./0.). I agree that it should throw an exception, as it is not clear whether it should be +Infinity or -Infinity, since 0 is neither positive nor negative.

stan gravatar imagestan ( 2016-03-09 09:44:13 -0600 )edit

@vdelecroix - your example is what I was trying to show, float(1.0)gives native Python floats as in your example. Not sure why Python doesn't conform to IEEE, despite reading several sources. I believe Numpy behaves "correctly" in this case, see

kcrisman gravatar imagekcrisman ( 2016-03-10 08:33:28 -0600 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2016-03-09 08:23:25 -0600

Seen: 177 times

Last updated: Mar 09 '16