Как реализована функция pow() в C++? Как ускорить её, не сильно теряя в точности?

Задался недавно данным вопросом из интереса и не нашёл нужной реализации. Настрочил кривой код на C++, стараясь не пользоваться встроенными функциями (из интереса).

Вот сам код:

double my_pow(double x, double n) // 'x' - основание, 'n' - показатель степени
{
    bool flip{ false }; // Не переворачиваем дробь при возвращении результата
    if (n < 0) // Если 'n' отрицательное число, то
    {
        n = -n; // 'n' становится положительным числом
        flip = true; // Переворачиваем дробь при возвращении результата
    }
    int fn = n; // целая часть показателя степени (будет использована в бинарном возведении в степень)
    double sn{ n - fn }; // дробная часть показателя степени (будет разложена в кучу квадратных корней)
    int c{}, s{}; // 'с' показывает сколько квадратных корней нужно взять, 's' - параметр для точности расчётов, если убрать, то возрастёт точность и снизится скорость

    // начало бинарного возведения в степень (обработка целой части показателя степени):
    double bp_x{ x }, res{ 1 }; 
    unsigned long long bp_n = fn;
    for (; bp_n != 0;) {
        if ((bp_n & 1) != 0) {
            res *= bp_x;
        }
        bp_x *= bp_x;
        bp_n >>= 1;
    }
    // конец бинарного возведения в степень

    // начало обработки дробной части степени
    for (; sn != 0 && s < 10;) // пока дробная часть не равна 0 или пока счётчик 's' меньше 10
    {
        sn *= 2; // умножаем дробную часть на 2
        c++; // увеличиваем счётчик квадратных корней на 1
        if (sn >= 1) // если после домножения дробная часть больше или равна единице, то
        {
            --sn; // вычитаем из дробной части 1
            ++s; // счётчик подобных вычитаний увеличиваем на 1
        
            // Тут начинается вычисление квадратных корней (обработка дробной части показателя степени)
            double x_sqrt{ x }; // копируется значение 'x'
            const double EPS = 1E-7; // константа, влияющая на точность вычисления квадратного корня (чем меньше, тем точнее)
            for (int i = 0; i < c; ++i) // метод выполняется 'c' количество раз
            {
                // Начало метода Ньютона для вычисления квадратного корня
                double xn = 1;
                for (;;) {
                    double nx = (xn + x_sqrt / xn) / 2;
                    if (((xn - nx < 0) ? -(xn - nx) : xn - nx) < EPS)  break;
                    xn = nx;
                }
                x_sqrt = xn;
                // Конец метода
            }
            // Заканчивается вычисление квадратных корней
            res *= x_sqrt; // домножаем 'res' на полученное число
        }
    }
    /*
    x ^ 0.5 - квадратный корень из 'x'

    Пример разложения дробной части в квадратные корни:
    0) 3 ^ 0.23
    умножим показатель степени на 2 и всё выражение возьмём в квадратный корень
    1) (3 ^ 0.46) ^ 0.5
    умножим показатель степени на 2 и всё выражение возьмём в квадратный корень
    3) ((3 ^ 0.92) ^ 0.5) ^ 0.5
    умножим показатель степени на 2 и всё выражение возьмём в квадратный корень
    4) (((3 ^ 1.84) ^ 0.5) ^ 0.5) ^ 0.5
    теперь показатель степени равен 1.84 >= 1, вынесем целую часть и получим
    5) (((3 ^ 0.5) ^ 0.5) ^ 0.5) * ((((3 ^ 0.84) ^ 0.5) ^ 0.5) ^ 0.5)
    умножим показатель степени на 2 и всё выражение возьмём в квадратный корень
    6) (((3 ^ 0.5) ^ 0.5) ^ 0.5)  *  (((((3 ^ 1.68) ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5)
    теперь показатель степени равен 1.68 >= 1, вынесем целую часть и получим
    7) (((3 ^ 0.5) ^ 0.5) ^ 0.5)  *  ((((3 ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5)  *  (((((3 ^ 0.68) ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5)
    умножим показатель степени на 2 и всё выражение возьмём в квадратный корень
    8) (((3 ^ 0.5) ^ 0.5) ^ 0.5)  *  ((((3 ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5)  *  ((((((3 ^ 1.36) ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5)
    теперь показатель степени равен 1.36 >= 1, вынесем целую часть и получим
    9) (((3 ^ 0.5) ^ 0.5) ^ 0.5)  *  ((((3 ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5)  *  (((((3 ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5)  *  ((((((3 ^ 0.36) ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5)
    Допустим, мы досчитали до нужной нам точности (лень дальше расписывать)
    Отбрасываем ((((((3 ^ 0.36) ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5) ^ 0.5) и вычисляем оставшееся выражение
    Ответ: 1.271657844...
    3 ^ 0.23 = 1.2874722841...
    Погрешность = ~1%
    Если увеличить количество итераций, то точность вычислений возрастёт.
    Расписал это, т.к. придумал такой способ, когда страдал фигнёй в калькуляторе, в интернете такого способа не нашёл.
    */
    return (flip) ? 1 / res : res; // Возвращаем результат
}

Данная функция в разы медленнее встроенной pow(), а ещё я не расписал тут некоторые случаи. Можно оптимизировать код, но мне кажется, что я попросту выбрал проигрышный алгоритм. Есть ли где-то +- нормальная реализация функции pow() или её исходный код (чтобы могла считать нецелые степени)?

UPD: Немного доработал код, теперь он не пересчитывает квадратные корни заново каждую итерацию, вместо этого он сохраняет последнее значение c в c_old (и потом начинает цикл с него) и копирует значение x в sqrt_x только один раз. Вот обновлённый код:

    double my_pow( double x, double n)
    {
        bool flip{ false };
        if (n < 0)
        {
            n = -n;
            flip = true;
        }
        int fn = n;
        double sn{ n - fn };
        int c{}, s{}, c_old{}; // добавил 'c_old', чтобы уменьшить количество вычислений
        double bp_x{ x }, res{ 1 }; 
        unsigned long long bp_n = fn;
        for (; bp_n != 0;)
        {
            if ((bp_n & 1) != 0)
            {
                res *= bp_x;
            }
            bp_x *= bp_x;
            bp_n >>= 1;
        }
        double x_sqrt{ x }; // определяем 'x_sqrt' и копируем в него значение 'x' только 1 раз
        for (; sn != 0 && s < 10;)
        {
            sn *= 2;
            c++;
            if (sn >= 1)
            {
                --sn;
                ++s;
                const double EPS = 1E-7;
                for (int i = c_old; i < c; ++i) // начинаем вычислять квадратные корни от последнего значения x_sqrt
                {
                    double xn = 1;
                    for (;;)
                    {
                        double nx = (xn + x_sqrt / xn) / 2;
                        if (((xn - nx < 0) ? -(xn - nx) : xn - nx) < EPS)  break;
                        xn = nx;
                    }
                    x_sqrt = xn;
                }
                c_old = c; // обновляем старое значение 'c'
                res *= x_sqrt;
            }
        }
        return (flip) ? 1 / res : res;
    }

На днях попробую решить данную задачу, используя ряды.