10

Мне не раз приходилось встречать рекомендации по типу "не пользуйтесь функцией std::pow". В связи с чем многие дают эту рекомендацию?

void MyPow(int& a, int n)
{
    int c = a;
    for (size_t i = 0; i < n; i++)
        a *= c;
}

Данная функция будет быстрее, чем стандартный std::pow? Или данную рекомендацию дают из-за того, что std::pow реализован не бинарно, а линейно (вот этого я кстати не знаю бинарный он или нет)?

1
  • потому что в стандартной бибилиотеке С++ отсутствует функция pow для целых чисел? 20 июн 2020 в 6:12

4 ответа 4

29
+50

Вынужден извиниться заранее - я хочу поговорить о применимости pow вообще, а не только для целочисленных значений.

Даже если вы хотите возводить в степень умножением - то не надо делать это столь прямолинейно, есть метод быстрого возведения в степень n со скоростью O(lg n).

Далее, как и любой совет, это совет, а не догма. Этот совет совершенно справедлив, например (и я его постоянно даю в таких ситуациях), когда начинают вычислять что-то типа pow(-1,n) (сами догадаетесь, как это вычислить быстро и точно?) или pow(x,2) - потому что даже в том же VC++ pow с целочисленной степенью в <cmath> реализована как (выбросил лишнее для понимания)

double pow(double _Xx, int _Yx) noexcept
{
    if (_Yx == 2) return (_Xx * _Xx);
    return pow(_Xx, static_cast<double>(_Yx));
}

Так какой смысл в постоянной проверке равно ли два - двум? :) При малых значениях степени также может оказаться, что непосредственное вычисление быстрее, чем вызов функции.

Если вы намекаете на шаблонную реализацию типа pow<int,int> - то, выбрасывая несущественное для понимания, в VC++ она выглядит так:

template<class _Ty1, class _Ty2,
    class = enable_if_t<is_arithmetic_v<_Ty1> && is_arithmetic_v<_Ty2>>> 
    _Common_float_type_t<_Ty1, _Ty2> pow(const _Ty1 _Left, const _Ty2 _Right)
    {   // bring mixed types to a common type
    using _Common = _Common_float_type_t<_Ty1, _Ty2>;
    return (pow(static_cast<_Common>(_Left), static_cast<_Common>(_Right)));
    }

Т.е. все равно сводится к обычной pow с плавающей точкой. Которая начинает выполнять ряд телодвижений по проверке аргументов и т.п., так что простая замена на exp(y*log(x)) работает несколько быстрее (впрочем, эта разница существенно зависит от используемой модели с плавающей запятой - у VC++ 2017 от практически равных при /fp:fast до разницы в 1.8 раза при /fp:precise). Кстати, думаю (точнее - знаю :), см. P.P.S.), если применить даже ваш линейный способ вычисления - он будет опережать стандартный до достаточно больших значений степени.

Точность при возведении в степень целочисленного значения также страдает, но об этом уже писали выше.

Словом, всякий инструмент хорош, когда правильно применен.

Еще одно замечание в связи с последней фразой - меня также бесит, когда начинают использовать pow для вычисления какого-нибудь ряда типа

введите сюда описание изображения

когда каждый член вычисляется возведением в степень, а не умножением на x, или когда так же в лоб вычисляют полином, игнорируя схему Горнера. Здесь применение pow глупо не потому, что она плоха, а потому, что здесь вообще не требуется возведение в степень!

P.S. А вообще, в программировании, как и во многих других областях деятельности - в том же кино масса примеров - сначала нечто начинают бездумно применять везде просто потому, что научились использовать это нечто. Потом приходит отрезвление - явный ведь перебор, может, вообще нужно отказаться от такой возможности?.. И только потом приходит понимание, что все хорошо в меру и на своем месте :) Но это так, отвлеченные размышления, не относящиеся к конкретно вопросу..

P.P.S. Не выдержал - заинтересовало, а в самом деле, когда будет быстрее использовать pow, чем просто линейное умножение? Набросал небольшой код, разово просчитал (VC++ 2017), построил график...

введите сюда описание изображения

Получается, где-то до 30 степени лучше просто множить, чем считать экспоненту от логарифма, и где-то до 50 - если использовать pow. Если использовать быстрое возведение в степень - то эта кривая на графике просто не видна, так как ее значение на всем диапазоне колеблется около 0.25-0.3 мс..

9
  • 8
    Браво! Побольше бы таких объяснений.
    – sneas
    20 июн 2020 в 6:33
  • 2
    Попробовал exp2(log2) - получилось еще быстрее, хотя у меня разница с pow не такая сильная: quick-bench.com/AM18GUV5_BUlA5rIPhB-M1VagvQ 21 июн 2020 в 19:09
  • @HolyBlackCat Я вот сейчас перезапустил бенчмарк по вашей ссылке, и получилось, что exp2(log2) даёт такое же время, как pow, а вот exp(log) при этом быстрее )) Фиг знает от чего это всё зависит, но результат явно неустойчивый )
    – CrazyElf
    21 окт 2020 в 7:04
  • pow и mul на графике сравниваются в районе 55. Самое маленькое целое число, которое имеет смысл возводить в степень - это двойка. А самый большой стандартный целочисленный тип - 64 бита. Итого pow обгонит mul только если вы считаете двойку в степенях от 55 до 63. Во всех остальных случаях или переполнение или mul быстрее. 4 дек 2022 в 21:23
  • @StanislavVolodarskiy Это если делать это целочисленно, но тогда уж 1<<n еще быстрее :). Но никто же не мешает возводить в степень число с плавающей точкой...
    – Harry
    9 дек 2022 в 5:01
2

В вашей реализации pow - ошибка: она дает неверный ответ для нулевого n (и для отрицательных). Не говоря уже о том, что она не поддерживает дробные n.

Стандартная pow - быстрее, она работает за константное время (независимое от n) и может быть реализована несколькими инструкциями процессора pow(x,y)==exp( log(x)*y ). Но! Она преобразует числа в double и обратно, что может привести к потере точности т.к. 64-bit int - содержит больше значащих цифр, чем мантисса в double.

6
  • точность много от чего зависит, регистры у современных x86 FPU 80-разрядные и точность сильно зависит от того, как долго оптимизатор держит значение в регистрах не выгружая в память. То есть с выключенной оптимизацией точность будет заведомо меньше 20 июн 2020 в 9:44
  • @PavelGridin, я думаю, что от того, выгружеы ли значения из регистров в память, может зависеть не точность, а скорость - ну, загрузит опять из памятив регситры. Потратится на несколько наносекунд больше. Но на точность это не должно влиять, иначе результат одного и того же вычисления был бы разным в разных случаях. Кстати, я думаю, сам совет "не пользуйтесь функцийей POW" тянется с тех далеких времен (борланд си и т.п.) когда сама функция pow уже была реализована через вычисления с экспонентой, а сопроцессоры были не на всех компах, и поэтому на типичном компьютере это тормозило.
    – S.H.
    20 июл 2020 в 0:54
  • @S.H., так и есть результат одного и того же вычисления в разных случаях разный 23 июл 2020 в 5:09
  • @PavelGridin, я читал несколько популярных статей (на тему сопроцессора, например habr.com/ru/post/503034). Но судя по этим стаьям там всё - детерминировано. Не могли бы Вы дать мне какую нибудь информацию о том, в каких случаях результат одного и того же вычисления в разных случаях разный? Спасибо.
    – S.H.
    23 июл 2020 в 15:18
  • надо взять что-нибудь интегральное, хороший пример БПФ написанное на C/C++, просто скомпилировать и сравнить результат для Release и Debug. Поэтому нельзя проверить правильность расчётного алгоритма простым бинарным сравнением, приходится сравнивать результаты статистическими методами 23 июл 2020 в 15:49
0

Мне не раз приходилось встречать рекомендации по типу "не пользуйтесь функцией Pow". В связи с чем многие дают эту рекомендацию?

Дурная рекомендация, дело нехитрое. Функцию pow() использовать можно, и даже обязательно нужно, как минимум, в юнит-тестах.

Я уж и не знаю, специально, с определённым умыслом, в вопросе дана именно такая, достаточно странная реализация функции MyPow(), или так просто случайно вышло. Но это хорошая иллюстрация почему pow() следует использовать. В самом деле, если взглянуть на результаты тестов:

#include <cassert>
#include <climits>
#include <cmath>
#include <iostream>
#include <type_traits>

using namespace std;

void MyPow(int& a, int n)
{
    int c = a;
    for (size_t i = 0; i < n; i++)
        a *= c;
}

template<typename I>
I ipow(I a, I b) {
    if(a == 1 || b == 0) {
        return 1;
    }
    
    if(is_signed<I>::value && b < 0) {
        return (a == 0 ? numeric_limits<I>::min() :
                a != -1 ? 0 : 
                b&1 ? -1 : 1);
    }
    
    I ret = 1;

    while(b) {
        if(b&1) {
            ret *= a;
        }
        a *= a;
        b >>= 1;
    }

    return ret;
}

#include <boost/core/demangle.hpp>
#include <typeinfo>

#ifdef COMPARE_POWL
    typedef long double dtype_t;
#else
    typedef double dtype_t;
#endif
static_assert(numeric_limits<dtype_t>::is_iec559, "Соответствие IEEE 754");

template<typename I>
void test_ipow(int verbose = 1) {
    I start = 0;

    if(is_signed<I>::value) {
        start = -5;
    }

    I finish = 0;

    while(pow<dtype_t>(finish, finish) <= (dtype_t)numeric_limits<I>::max() &&
          pow<dtype_t>(finish, finish) 
                    <= pow<dtype_t>(2, 1 + numeric_limits<dtype_t>::digits)) {
        finish++;
    }
    --finish;

    if(verbose) {
        cout << "test_ipow<" 
            << boost::core::demangle(typeid(I).name()).c_str()
            << ">(): start " << start
            << " finish " << finish
            << '\n';
    }

    for(I a = start; a <= finish; a++) {
        for(I b = start; b <= finish; b++) {
            if(verbose >= 2) {
                cout << a << ", " << b
                    << " -> " << ipow(a, b) 
                    << ", " << pow<dtype_t>(a, b) 
                    << ", " << (I)pow<dtype_t>(a, b)
                    << '\n';
            }
            assert(ipow(a, b) == (I)pow<dtype_t>(a, b));
        }
    }

    const I check_start = 3;
    const I check_finish = (I)trunc(cbrt(numeric_limits<I>::max()));
    intmax_t total_power = 0;
    int max_diff = 0;
    int min_diff = 0;
    intmax_t sum_diff = 0;
    intmax_t cnt_diff = 0;

    for(I a = check_start; a <= check_finish; a++) {
        I fi = (I)trunc(log((dtype_t)numeric_limits<I>::max())/log(a));

        for(I b = check_start; b <= fi; b++) {
            dtype_t d = pow<dtype_t>(a, b);
            I i = ipow(a, b);
            I diff = i - (I)d;
            total_power++;
            if(diff) {
                int idiff = (int)diff;
                assert(abs(idiff) <= d*numeric_limits<decltype(d)>::epsilon());
                if(verbose) {
                    cnt_diff++;
                    sum_diff += idiff;
                    if(idiff < min_diff) {
                        min_diff = idiff;
                    }
                    if(idiff > max_diff) {
                        max_diff = idiff;
                    }
                }
            }
        }
    }
    if(verbose) {
        cout << "Total power: " << total_power
            << " Total diff: " << cnt_diff
            << " Mean diff: " << (dtype_t)sum_diff/cnt_diff 
            << " Min: " << min_diff
            << " Max: " << max_diff
            << '\n';
    }
}

int main(int ac, char *av[]) {
    for(int i = 0; i <= 3; i++) {
        int a = 2;
        MyPow(a, i);
        cout << "MyPow(2, " << i
            << ") -> " << a
            << " - " << pow(2, i)
            << (a == pow(2, i) ? ": норма\n" : ": ОШИБКА\n");
    }

    test_ipow<int>();
    test_ipow<unsigned>();
    test_ipow<int64_t>();
    test_ipow<uint64_t>();

    return 0;
}

Результат:

$ clang++-mp-15 -O3 -I/opt/local/include MyPow.cpp && ./a.out 
MyPow(2, 0) -> 2 - 1: ОШИБКА
MyPow(2, 1) -> 4 - 2: ОШИБКА
MyPow(2, 2) -> 8 - 4: ОШИБКА
MyPow(2, 3) -> 16 - 8: ОШИБКА
test_ipow<int>(): start -5 finish 9
Total power: 1669 Total diff: 0 Mean diff: nan Min: 0 Max: 0
test_ipow<unsigned int>(): start 0 finish 9
Total power: 2067 Total diff: 0 Mean diff: nan Min: 0 Max: 0
test_ipow<long long>(): start -5 finish 14
Total power: 2161089 Total diff: 1587677 Mean diff: 4.30641 Min: -515 Max: 520
test_ipow<unsigned long long>(): start 0 finish 14
Total power: 2717817 Total diff: 2108847 Mean diff: 3.56374 Min: -1032 Max: 1040

Можно увидеть, что аргументация против использования функции pow() так себе:

  1. Точность. Для типа int для которого предназначена MyPow() функция pow() абсолютно точна, а функция MyPow() всегда выдаёт странный результат. Функция powl() для Intel/AMD без SSE/AVX абсолютно точна до uint64_t;
  2. Понятность, доказательность, корректность интерфейса. Функция pow() - часть стандартов POSIX/C/C++, а MyPow(), прошу прощения, но это ж "колхоз на коленке". Любой, кто читает вызов MyPow(a, 0) или MyPow(a, 1) некоторое время будет считать, что это возведение в 0 и в 1 степень, но, через некоторое время будет обескуражен, что это совсем не так;
  3. Производительность. И в части производительности MyPow() не блещет, во многом, благодаря своему странному интерфейсу, т.е. если компилятор видит её определение и может обеспечить inline подстановку, то ещё более менее, для неслишком больших степеней даже сможет немного обогнать pow(), но при раздельной трансляции - беда;

Понятно, что как любая библиотечная функция pow() является определённым компромиссом и в каждом конкретном случае может быть, в идеале, заменена на более лучший вариант. Но это ж только в идеале, в частности, при условии понимания автором кода численных методов и соответствующем документировании кода.

В части использования pow() для плавающих и комплексных чисел, можно отметить:

  • Почти всегда, pow(x, 2) и pow(x, 3) можно заменить на x*x и x*x*x, хотя я и не уверен в пользе такой замены для pow(x, 2);
  • Также не стоит использовать pow(x, 0.5) и pow(x, 1./3.), вместо sqrt(x) и cbrk(x), хотя для pow(x, 0.5) это вопрос только производительности;
  • В соседних ответах сравнивают производительности pow(x, y) с exp(y*log(x)) или с exp2(y*log2(x)), но такое сравнение имеет весьма ограниченную ценность, поскольку даже выражение exp2(y*log2(x)) всегда вычисляется с гораздо худшей точностью. Если диапазон данных ограничен, то можно рассмотреть замену pow(x, y) на комбинацию expm1() и log1p(), но это ж для тех, кто понимает.

Уточнение по результату замечания @wololo, добавлен static_assert(numeric_limits<dtype_t>::is_iec559, "Соответствие IEEE 754");.

22
  • 2
    Для типа int ... функция pow() абсолютно точна. Функция powl() ... абсолютно точна до uint64_t Ни стандарт арифметики с плавающей точкой IEEE 754, ни стандарты языков C и C++ не требуют абсолютно точного результата для функции pow, даже если точный результат представим типами double, int или uint64_t. Это сильнейший аргумент против использования данной функции в целочисленных выражениях. И реализации с неаккуратной имплементацией данной функции существуют: 1, 2.
    – wololo
    5 дек 2022 в 21:39
  • @wololo, насчёт IEEE 754, это Вы совершенно зря, функция pow() входит в раздел "9.2 Recommended correctly rounded functions", т.е. если не происходит исключения inexact (9.1.1), результат должен быть точен и правильно округлён. Что до ошибок, то известны даже ЦП у которых операция умножения с ошибкой, как, впрочем, и в некоторых версиях SunOS библиотечная арифметическая операция была с ошибкой (лично находил, Sun hot-fix делал). Тут, уж все под Богом, но и тестировать программы надо, а не абы как.
    – Serge3leo
    5 дек 2022 в 22:18
  • 2
    1) 9. Recommended operations. Clause 5 completely specifies the operations required for all supported arithmetic formats. This clause specifies additional operations, recommended for all supported arithmetic formats. 2) A conforming function shall return results correctly rounded for the applicable rounding direction for all operands in its domain. The preferred quantum is language-defined. 3) 9.2 Recommended correctly rounded functions ... names of the operations in Table 9.1 do not necessarily correspond to the names that any particular programming language would use.
    – wololo
    5 дек 2022 в 22:44
  • 2
    1) POSIX/C/C++, по-моему, холодно ссылаются на IEEE 754. Стандарт языка C ссылается вполне конкретно в Annex F (normative) IEC 60559 floating-point arithmetic. В разделе F.3 Operators and functions перечислены операторы и функции, выполняемые в соответствии с IEEE 754, и функции pow там нет, как и подавляющего большинства (если не всех) функций из таблицы рекомендованных в IEEE 754. Раздел F.10 Mathematics <math.h> повторно поясняет, что функции, напрямую специфицируемые IEEE 754, перечислены в F.3: 3 The functions that IEC 60559 specifies directly are identified in F.3.
    – wololo
    6 дек 2022 в 0:49
  • 2
    В черновике стандарта C23 пошли ещё дальше и всю тригонометрию, логарифмы, экспоненты и конечно же pow свели в одну таблицу и подписали: The C functions in the following table correspond to mathematical operations recommended by IEC 60559. However, correct rounding, which IEC 60559 specifies for its operations, is not required for the C functions in the table. 2) numeric_limits<dtype_t>::is_iec559 — это не обязывает функцию pow считаться правильно.
    – wololo
    6 дек 2022 в 0:50
0

Все просто: функция std::pow принимает в качестве параметра число с плавающей запятой, двойное или длинное двойное число и возвращает то же самое в качестве возвращаемого типа. Это вообще не целочисленно-совместимая функция. Благодаря неявным приведениям (превращению одного типа в другой) это все равно может работать.

Я бы посоветовал вам открыть сайт https://godbolt.org/ и поэкспериментировать с ним, изучая дизассемблированный код.

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.