SymPy

Depending on your goal, it might be a good idea to delay the calculation of square roots for as long as possible. SymPy might help.

SymPy is a Python library for symbolic mathematics.

import sympy
sympy.sqrt(2)
# => sqrt(2)

This doesn't seem very useful at first.

But sympy can give more information than floats or Decimals:

sympy.sqrt(8) / sympy.sqrt(27)
# => 2*sqrt(6)/9

Also, no precision is lost. (√2)² is still an integer:

s = sympy.sqrt(2)
s**2
# => 2
type(s**2)
#=> <class 'sympy.core.numbers.Integer'>

In comparison, floats and Decimals would return a number which is very close to 2 but not equal to 2:

(2**0.5)**2
# => 2.0000000000000004

from decimal import Decimal
(Decimal('2')**Decimal('0.5'))**Decimal('2')
# => Decimal('1.999999999999999999999999999')

Sympy also understands more complex examples like the Gaussian integral:

from sympy import Symbol, integrate, pi, sqrt, exp, oo
x = Symbol('x')
integrate(exp(-x**2), (x, -oo, oo))
# => sqrt(pi)
integrate(exp(-x**2), (x, -oo, oo)) == sqrt(pi)
# => True

Finally, if a decimal representation is desired, it's possible to ask for more digits than will ever be needed:

sympy.N(sympy.sqrt(2), 1_000_000)
# => 1.4142135623730950488016...........2044193016904841204

NumPy

>>> import numpy as np
>>> np.sqrt(25)
5.0
>>> np.sqrt([2, 3, 4])
array([1.41421356, 1.73205081, 2.        ])

docs

Negative

For negative reals, it'll return nan, so np.lib.scimath.sqrt() is available for that case.

>>> a = np.array([4, -1, np.inf])
>>> np.sqrt(a)
<stdin>:1: RuntimeWarning: invalid value encountered in sqrt
array([ 2., nan, inf])
>>> np.lib.scimath.sqrt(a)
array([ 2.+0.j,  0.+1.j, inf+0.j])

Another option, of course, is to convert to complex first:

>>> a = a.astype(complex)
>>> np.sqrt(a)
array([ 2.+0.j,  0.+1.j, inf+0.j])