Три или четыре способа извлечь корень

Если коротко, то разница в реализации. И очевидно в скорости выполнения, входных и выходных параметрах. А вообще-то интересно заглянуть, что там во встроенных функциях.

pow(x, n)

 @overload
    def pow(base: int, exp: int, mod: None = ...) -> Any: ...  # returns int or float depending on whether exp is non-negative
    @overload
    def pow(base: int, exp: int, mod: int) -> int: ...
    @overload
    def pow(base: float, exp: float, mod: None = ...) -> float: ...
    @overload
    def pow(base: _SupportsPow2[_E, _T_co], exp: _E) -> _T_co: ...
    @overload
    def pow(base: _SupportsPow3[_E, _M, _T_co], exp: _E, mod: _M) -> _T_co: ...

math.sqrt()

def sqrt(__x: SupportsFloat) -> float: ...

math.pow()

def pow(__x: SupportsFloat, __y: SupportsFloat) -> float: ...

И немного документации math.pow(). Здесь хотелось бы перевести одно заключение из документации:

В отличие от встроенного оператора **, math.pow () преобразует оба своих аргумента в тип float. Используйте ** или встроенную функцию pow () для вычисления точных целочисленных степеней.

А так же можно пользоваться и cmath модулем, для комплексных чисел и не только. Он возращает комплексное число, но при получении квадратного корня из 4, например, мнимая часть будет равняться 0j.

def sqrt(__z: _C) -> complex: ...

Можно сделать вывод, что под каждую конкретную задачу может быть подобран оптимальная функция для вычислений. Думаю, что тема достаточно обширна и можно исследовать, например, скорости для вычисления в том или ином случае.

Все вышеописанное для версии python 3.х.x.

upd

И все таки я решил посмотреть глубже. Не уверен в верности моей попытки, но все же (скажите, где я ошибся). Я решил исследовать две built-in функции: pow() и math.pow(). Найти их исходный код можно pow и math.pow.

Для pow() были найдены данные строки

/*[clinic input]
pow as builtin_pow
    base: object
    exp: object
    mod: object = None
Equivalent to base**exp with 2 arguments or base**exp % mod with 3 arguments
Some types, such as ints, are able to use a more efficient algorithm when
invoked using the three argument form.
[clinic start generated code]*/

static PyObject *
builtin_pow_impl(PyObject *module, PyObject *base, PyObject *exp,
                 PyObject *mod)
/*[clinic end generated code: output=3ca1538221bbf15f input=435dbd48a12efb23]*/
{
    return PyNumber_Power(base, exp, mod);
}

Что это может означать? В данной документации можно найти, как это работает. А конкретно в документации написано, что используется препроцессор для CPython и в данной конструкции "магическим образом преобразуются в переменные С с типом данных, которые нам нужны"(вольный перевод.). Функция PyNumber_Power() импортируется из DLL библиотеки (могу ошибаться, т.к. нашел только вызов в python3dll.c-файле, строка 409, EXPORT_FUNC(PyNumber_Power)). Вот в принципе и все.

Для math.pow() реализация немного интересней:

/*[clinic input]
math.pow
    x: double
    y: double
    /
Return x**y (x to the power of y).
[clinic start generated code]*/

static PyObject *
math_pow_impl(PyObject *module, double x, double y)
/*[clinic end generated code: output=fff93e65abccd6b0 input=c26f1f6075088bfd]*/
{
    double r;
    int odd_y;

    /* deal directly with IEEE specials, to cope with problems on various
       platforms whose semantics don't exactly match C99 */
    r = 0.; /* silence compiler warning */
    if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
        errno = 0;
        if (Py_IS_NAN(x))
            r = y == 0. ? 1. : x; /* NaN**0 = 1 */
        else if (Py_IS_NAN(y))
            r = x == 1. ? 1. : y; /* 1**NaN = 1 */
        else if (Py_IS_INFINITY(x)) {
            odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
            if (y > 0.)
                r = odd_y ? x : fabs(x);
            else if (y == 0.)
                r = 1.;
            else /* y < 0. */
                r = odd_y ? copysign(0., x) : 0.;
        }
        else if (Py_IS_INFINITY(y)) {
            if (fabs(x) == 1.0)
                r = 1.;
            else if (y > 0. && fabs(x) > 1.0)
                r = y;
            else if (y < 0. && fabs(x) < 1.0) {
                r = -y; /* result is +inf */
            }
            else
                r = 0.;
        }
    }
    else {
        /* let libm handle finite**finite */
        errno = 0;
        r = pow(x, y);
        /* a NaN result should arise only from (-ve)**(finite
           non-integer); in this case we want to raise ValueError. */
        if (!Py_IS_FINITE(r)) {
            if (Py_IS_NAN(r)) {
                errno = EDOM;
            }
            /*
               an infinite result here arises either from:
               (A) (+/-0.)**negative (-> divide-by-zero)
               (B) overflow of x**y with x and y finite
            */
            else if (Py_IS_INFINITY(r)) {
                if (x == 0.)
                    errno = EDOM;
                else
                    errno = ERANGE;
            }
        }
    }

    if (errno && is_error(r))
        return NULL;
    else
        return PyFloat_FromDouble(r);
}

Можно увидеть не только использования Argument Clinic-препроцессора, но и обработку некоторых крайних значений (бесконечность ** бесконечность, например) и обработку ошибок.

Можно копать дальше, не уверен есть ли необходимость. Но инетересно же) можно так же увидеть еще одну интересную реализацию для long_pow() здесь

К сожалению, в Google Colab нет функции math.isqrt, потому что там Python 3.7.11. А насчёт остальных функций можно проверить скорость так (небольшой цикл добавлен, чтобы не было кэширования, а то без цикла %timeit ругается, что возможно результаты были закэшированы, хотя они примерно так же соотносятся, я проверял, но лучше принять меры предосторожности):

import math
n = 100
%timeit [x ** 0.5 for x in range(n)]
%timeit [pow(x, 0.5) for x in range(n)]
%timeit [math.sqrt(x) for x in range(n)]

Вывод:

100000 loops, best of 5: 17.7 µs per loop
10000 loops, best of 5: 23.9 µs per loop
100000 loops, best of 5: 13.2 µs per loop

Самая быстрая - специализированная функция math.sqrt, что логично. Оператор ** чуть более медленный, а функция pow самая медленная, что, в общем, давно известно. Функцию pow имеет смысл применять только с третьим аргументом (остаток по модулю), тогда её использование может быть оправдано.
Также проверим, а выдают ли эти варианты одинаковые результаты. Хотя бы на каком-то списке натуральных чисел:

import math
n = 10000
y1 = [x ** 0.5 for x in range(n)]
y2 = [pow(x, 0.5) for x in range(n)]
y3 = [math.sqrt(x) for x in range(n)]
print(y1 == y2 == y3, sum(y1) == sum(y2) == sum(y3))

Вывод:

True True

Ещё можно посмотреть в какой код переводятся эти варианты и убедиться, что код разный:

dis.dis('x ** 0.5')

  1           0 LOAD_NAME                0 (x)
              2 LOAD_CONST               0 (0.5)
              4 BINARY_POWER
              6 RETURN_VALUE

dis.dis('pow(x, 0.5)')

  1           0 LOAD_NAME                0 (pow)
              2 LOAD_NAME                1 (x)
              4 LOAD_CONST               0 (0.5)
              6 CALL_FUNCTION            2
              8 RETURN_VALUE

dis.dis('math.sqrt(x)')

  1           0 LOAD_NAME                0 (math)
              2 LOAD_METHOD              1 (sqrt)
              4 LOAD_NAME                2 (x)
              6 CALL_METHOD              1
              8 RETURN_VALUE

Если скорость важна, то для извлечения корня видимо нужно использовать специально заточенную под это функцию math.sqrt. В общем же случае проще на все случаи жизни использовать оператор **. И в случае, когда нас интересует вариант с возведением в степень и дальнейшим взятием остатка по модулю - тут лучше всего функция pow.

P.S. Ну и, конечно, если нужно получить сразу много корней, то лучше использовать векторизованный Numpy для этого. Использование для задачи range/sqrt/sum векторизированных функций даёт 50-кратное ускорение:

import math
import numpy as np

n = 1000000
%timeit sum_py = sum(math.sqrt(x) for x in range(n))
%timeit sum_nu = np.sum(np.sqrt(np.arange(n)))

Вывод:

10 loops, best of 5: 157 ms per loop
100 loops, best of 5: 3.43 ms per loop

Проверка, что результат сошёлся (в данном случае - приблизительно):

import math
import numpy as np

n = 1000000
sum_py = sum(math.sqrt(x) for x in range(n))
sum_nu = np.sum(np.sqrt(np.arange(n)))
print(sum_py, abs(sum_py-sum_nu))

Вывод:

666666166.4588418 1.9550323486328125e-05

Сравнение немного с другой колокольни.

1. math.isqrt(x) не имеет отношения, ни к IEEE 754, ни к POSIX, ни к C/C++. Сделано специально для Python, вероятно - хорошо. Мало того, ни одна из альтернатив не поддерживает int без формального ограничения размера и не сможет что-то в духе: math.isqrt((2**129 + 1)**2) - math.isqrt((2**129)**2) == 1;
2. math.sqrt(x) по документации Python отсылка к типам double и IEEE 754. А реализовано средствами библиотеки C и/или POSIX. Поскольку все трое сходятся на требовании корректности округления для sqrt(), то применимо к целым аргументам до 2⁵³ (numpy.sqrt(numpy.float128(x)) до 2⁶³);
3. x**0.5, pow(x, 0.5) и math.pow(x, 0.5), аналогично п.2, но в стандарте C есть оговорка, что для pow() корректное округление желательно, но необязательно. В принципе, даже у весьма и весьма хреновой glibc, который напрочь игнорирует требования POSIX и IEEE 754, по документации, точность 1 ulp, т.е. если точный результат представим, то он и будет результатом. К примеру, для квадратов целых чисел до 2⁵³, но за полумотонность glibc умалчивает.

А самый быстрый, обычно п.2, тем более, что он часто имеет поддержку со стороны процессора. Но это, конечно, существенно заметно для аргументов массивов numpy.sqrt().