Как извлечь кубический корень как можно точнее?

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

Пример

Код

import numpy as np 

cubes = [125, -64, 27, -8, 1] 
cube_roots = np.cbrt(cubes) 
print(cube_roots)

Вывод

[ 5. -4.  3. -2.  1.]

Возможно вы хотели просто решить проблему неточного извлечения корней из точных кубов, но получилось как получилось:

import struct
from fractions import Fraction


def difference(float_expected, float_root, degree):
    '''
    Возвращает разницу между степенью корня и ожидаемым значением
    '''

    fraction_expected = Fraction.from_float(float_expected)
    fraction_root = Fraction.from_float(float_root)
    return abs(fraction_expected - fraction_root ** degree)


def ternary_search(lower_bound, upper_bound, function):
    '''
    Находит целое значение x из промежутка [lower_bound, upper_bound),
    для которого значение function(x) минимально
    '''

    while lower_bound + 3 < upper_bound:
        mid1 = (lower_bound * 2 + upper_bound) // 3
        mid2 = (lower_bound + upper_bound * 2) // 3
        res1 = function(mid1)
        res2 = function(mid2)
        if res1 < res2:
            upper_bound = mid2
        else:
            lower_bound = mid1

    return min(range(lower_bound, upper_bound), key=function)


def create_float(exponent, mantissa):
    '''
    Создает число binary64 IEEE-754 по экспоненте и мантиссе
    '''
    
    number = (exponent << 52) + mantissa
    return struct.unpack('d', number.to_bytes(8, 'little'))[0]


def precise_root(float_value, degree):
    '''
    Возвращает максимально точное значение корня среди всех возможных
    значений binary64 IEEE-754 (используя только стандартный формат
    1-11-52, исключая все субнормальные и специальные значения)

    value - положительное число с плавающей запятой
    degree - положительное целое число больше единицы
    '''
    
    nearest_values = [float_value ** (1 / degree)]
    min_difference = difference(float_value, nearest_values[0], degree)
    for exponent in range(1023, 2047):
        calc_difference = lambda x: difference(float_value,
                                               create_float(exponent, x),
                                               degree)
        mantissa = ternary_search(0, 2 ** 52, calc_difference)
        diff = difference(float_value, create_float(exponent, mantissa), degree)
        if diff < min_difference:
            nearest_values = [create_float(exponent, mantissa)]
        elif diff == min_difference:
            nearest_values.append(create_float(exponent, mantissa))

    return nearest_values[0]


if __name__ == "__main__":
    value = 1000000001.0
    degree = 9
    root = precise_root(value, degree)
    diff = difference(value, root, degree)
    print(format(root, '.70f'))
    print(format(value ** (1 / degree), '.70f'))
    print(format(float(diff), '.70f'))
    print(format(float(difference(value, value ** (1 / degree), degree)), '.70f'))
    # 10.0000000011111112030448566656559705734252929687500000000000000000000000
    # 10.0000000011111094266880172654055058956146240234375000000000000000000000
    # 0.0000000831848155171970377866333596662828941958878203877247869968414307
    # 0.0000015155363413641129334422444352448167137481505051255226135253906250

Ну вот такая есть библиотека, например:

import gmpy2

print(gmpy2.root(19683, 3))

Вывод:

27.0

Начиная с Python 3.11 всё есть из коробки, поэтому стало не обязательно следовать совету @DGDays : использовать numpy.

1. В рамках точности стандартного типа float это выполняет функция math.cbrt(x). Если базовая система в точности удовлетворяет стандартам IEC 60559/IEEE 754/POSIX, то cbrt(x) является реализации функции rootx(x, n) для случая n=3, и возвращает правильно округлённый результат точного значения кубического корня. Правда, стандарт языка C (ISO/IEC 9899) предоставляет послабления реализациям в деле вычисления спецфункций, поэтому к этому вопросу следует подходить осторожно, с предварительным тестированием;

2. Модуль стандартной библиотеки decimal содержит реализацию десятичных действительных чисел с произвольной точностью. К сожалению, реализация IEEE 854 в нём неполна, но значение степени 1/3, как и операцию ** возведение в степень, можно вычислить с произвольной точностью.

Сравнение этих двух вариантов:

import decimal
import math
import sys

dig = sys.float_info.dig
bigdig = 3*dig
tcbrt = math.cbrt
#tcbrt = glibc_math_cbrt
#tcbrt = glibc_np_cbrt
sabdfc, sumdfc, maxdfc, mindfc = 0, 0, -math.inf, math.inf
verbose = 1

with decimal.localcontext(prec=1000):
    d1_3 = decimal.Decimal(1)/3
    for b in (list(range(0, 100, 10)) +
              [n**3 for n in range(300)] +
              list(range(10**dig, 10**(dig+1), 10**dig)) +
              list(range(10**bigdig, 10**(bigdig+1), 10**bigdig))):
        for i in [b, -b, b*(1+0.000005), -b*(1-0.000005)]:
            fc = tcbrt(i)
            dc = (decimal.Decimal(abs(i))**d1_3).copy_sign(decimal.Decimal(i))
            dfc = (dc - decimal.Decimal(fc))/decimal.Decimal(math.ulp(fc))
            sabdfc += math.fabs(dfc)
            sumdfc += dfc
            maxdfc = max(maxdfc, dfc)
            mindfc = min(mindfc, dfc)
            w = ""
            if math.fabs(dfc) >= 1:
                w = "\033[91m ERROR >= 1 ulp \033[0m" 
            elif isinstance(i, int) and math.fabs(dfc) > 0.5:
                w = "\033[93m warning isinstance(i, int) and > 0.5 ulp\033[0m"
            elif verbose >= 2 and math.fabs(dfc) > 0.5:
                w = " warning > 0.5 ulp"
            if verbose >= 3 or w:
                print(f"{dfc: 7.3f} ulp for {i}{w}")
print(f"{sabdfc: .5g}, {sumdfc: .5g}, {maxdfc: .5g}, {mindfc: .5g}")

Результат для Python 3.12 Intel macOS:

 -0.518 ulp for 6000000000000000 warning isinstance(i, int) and > 0.5 ulp
  0.518 ulp for -6000000000000000 warning isinstance(i, int) and > 0.5 ulp
 -0.595 ulp for 7000000000000000000000000000000000000000000000 warning isinstance(i, int) and > 0.5 ulp
  0.595 ulp for -7000000000000000000000000000000000000000000000 warning isinstance(i, int) and > 0.5 ulp
 187.23,  7.0225,  0.67293, -0.63271

P.S.

ВАЖНОЕ ПРЕДУПРЕЖДЕНИЕ: оговорка в п.1, увы, оказалась часто встречающейся. Как отметили в комментариях, в glibc, увы, имеется грубая ошибка реализации double cbrt(x), срабатывающая для math.cbrt(27**3) (для math.cbrt(4**3) или math.cbrt(25**3) она не срабатывает), которую наследуют версии Python под большинство Linux. Естественно Python под Windows (распространенные сборки: с официального сайта, GitHub actions и т.п.), как и под macOS, FreeBSD, ... не имеют таковой ошибки.

Если под Linux нужна точность float, а пересобирать сам Python или только NumPy компиляторами c приличными библиотеками Intel, NAG и др. не вариант, то можно использовать:

import math

def glibc_math_cbrt(A):
    # Обходим ошибку glibc
    xk = math.cbrt(A)
    # DEBUG: import sys
    # DEBUG: import numpy
    # DEBUG: xk *= 1. + 1e7*sys.float_info.epsilon*(numpy.random.rand() - 0.5)
    # Дополнительная итерация метода Ньютона примерно,
    # как минимум, удваивает число значимых цифр для 
    # определённого ненулевого и не бесконечного результата.
    return xk + ((A/xk**2 - xk)/3) if math.isfinite(xk) and xk else xk
    
[glibc_math_cbrt(A) for A in (-19683, 0, -0., 19683, 
                                 -math.inf, math.inf, math.nan)]

Обновил код и результат сравнения (для glibc_np_cbrt аналогично, но чуть сложнее ввиду того, что аргумент может быть массивом).

P.P.S.

Насчёт совета @DGDays, поскольку в glibc аналогичная ошибка в long double cbrtl(x) всё таки была исправлена, то numpy.cbrt(numpy.float128(19683)) должен выдавать пристойный результат 27.0 и под Linux тоже.

P.P.P.S.

Забавно, что на нормальных системах (с нормальными библиотеками, по POSIX.1 и др.):

• math.cbrt(27**3) выдаёт 27.0 (ошибка 0 ulp);
• (27**3)**(1/3) выдаёт тоже самое, что и (27**3)**0.333333333333333315 - 26.999999999999996 (ошибка возведения в степень -0.390 ulp при ошибке показателя 0.333 ulp).

А на glibc, увы:

• math.cbrt(27**3) выдаёт 26.999999999999996 (ошибка 1 ulp);
• (27**3)**(1/3) выдаёт тоже самое, что и (27**3)**0.333333333333333315 - 26.999999999999996 (ошибка возведения в степень -0.390 ulp при ошибке показателя 0.333 ulp).