# numpy division with RuntimeWarning: invalid value encountered in double_scalars

I wrote the following script:

``````import numpy

d = numpy.array([[1089, 1093]])
e = numpy.array([[1000, 4443]])
print res
``````

But I got this result and with the error occurred:

``````nan
C:UsersDesktoptest.py:16: RuntimeWarning: invalid value encountered in double_scalars
``````

It seems to be that the input element were too small that python turned them to be zeros, but indeed the division has its result.

How to solve this kind of problem?

You can’t solve it. Simply `answer1.sum()==0`, and you can’t perform a division by zero.

This happens because `answer1` is the exponential of 2 very large, negative numbers, so that the result is rounded to zero.

`nan` is returned in this case because of the division by zero.

Now to solve your problem you could:

• go for a library for high-precision mathematics, like mpmath. But that’s less fun.
• as an alternative to a bigger weapon, do some math manipulation, as detailed below.
• go for a tailored `scipy/numpy` function that does exactly what you want! Check out @Warren Weckesser answer.

Here I explain how to do some math manipulation that helps on this problem. We have that for the numerator:

``````exp(-x)+exp(-y) = exp(log(exp(-x)+exp(-y)))
= exp(log(exp(-x)*[1+exp(-y+x)]))
= exp(log(exp(-x) + log(1+exp(-y+x)))
= exp(-x + log(1+exp(-y+x)))
``````

where above `x=3* 1089` and `y=3* 1093`. Now, the argument of this exponential is

`-x + log(1+exp(-y+x)) = -x + 6.1441934777474324e-06`

For the denominator you could proceed similarly but obtain that `log(1+exp(-z+k))` is already rounded to `0`, so that the argument of the exponential function at the denominator is simply rounded to `-z=-3000`. You then have that your result is

``````exp(-x + log(1+exp(-y+x)))/exp(-z) = exp(-x+z+log(1+exp(-y+x))
= exp(-266.99999385580668)
``````

which is already extremely close to the result that you would get if you were to keep only the 2 leading terms (i.e. the first number `1089` in the numerator and the first number `1000` at the denominator):

``````exp(3*(1089-1000))=exp(-267)
``````

For the sake of it, let’s see how close we are from the solution of Wolfram alpha (link):

``````Log[(exp[-3*1089]+exp[-3*1093])/([exp[-3*1000]+exp[-3*4443])] -> -266.999993855806522267194565420933791813296828742310997510523
``````

The difference between this number and the exponent above is `+1.7053025658242404e-13`, so the approximation we made at the denominator was fine.

The final result is

``````'exp(-266.99999385580668) = 1.1050349147204485e-116
``````

``````1.105034914720621496.. × 10^-116 # Wolfram alpha.
``````

and again, it is safe to use numpy here too.

You can use `np.logaddexp` (which implements the idea in @gg349’s answer):

``````In : d = np.array([[1089, 1093]])

In : e = np.array([[1000, 4443]])

In : log_res
Out: -266.99999385580668

In : res = exp(log_res)

In : res
Out: 1.1050349147204485e-116
``````

Or you can use `scipy.special.logsumexp`:

``````In : from scipy.special import logsumexp

In : res = np.exp(logsumexp(-3*d) - logsumexp(-3*e))

In : res
Out: 1.1050349147204485e-116
``````