Почему не стоит пользоваться функцией std::pow для целых чисел?

Question

Мне не раз приходилось встречать рекомендации по типу "не пользуйтесь функцией 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 реализован не бинарно, а линейно (вот этого я кстати не знаю бинарный он или нет)?

Answer 1

Вынужден извиниться заранее - я хочу поговорить о применимости 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 мс..

Answer 2

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

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

Answer 3

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

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

Answer 4

Мне не раз приходилось встречать рекомендации по типу "не пользуйтесь функцией 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&lt;I&gt;::value &amp;&amp; b &lt; 0) {
    return (a == 0 ? numeric_limits&lt;I&gt;::min() :
            a != -1 ? 0 : 
            b&amp;1 ? -1 : 1);
}

I ret = 1;

while(b) {
    if(b&amp;1) {
        ret *= a;
    }
    a *= a;
    b &gt;&gt;= 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&lt;I&gt;::value) {
    start = -5;
}

I finish = 0;

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

if(verbose) {
    cout &lt;&lt; &quot;test_ipow&lt;&quot; 
        &lt;&lt; boost::core::demangle(typeid(I).name()).c_str()
        &lt;&lt; &quot;&gt;(): start &quot; &lt;&lt; start
        &lt;&lt; &quot; finish &quot; &lt;&lt; finish
        &lt;&lt; '\n';
}

for(I a = start; a &lt;= finish; a++) {
    for(I b = start; b &lt;= finish; b++) {
        if(verbose &gt;= 2) {
            cout &lt;&lt; a &lt;&lt; &quot;, &quot; &lt;&lt; b
                &lt;&lt; &quot; -&gt; &quot; &lt;&lt; ipow(a, b) 
                &lt;&lt; &quot;, &quot; &lt;&lt; pow&lt;dtype_t&gt;(a, b) 
                &lt;&lt; &quot;, &quot; &lt;&lt; (I)pow&lt;dtype_t&gt;(a, b)
                &lt;&lt; '\n';
        }
        assert(ipow(a, b) == (I)pow&lt;dtype_t&gt;(a, b));
    }
}

const I check_start = 3;
const I check_finish = (I)trunc(cbrt(numeric_limits&lt;I&gt;::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 &lt;= check_finish; a++) {
    I fi = (I)trunc(log((dtype_t)numeric_limits&lt;I&gt;::max())/log(a));

    for(I b = check_start; b &lt;= fi; b++) {
        dtype_t d = pow&lt;dtype_t&gt;(a, b);
        I i = ipow(a, b);
        I diff = i - (I)d;
        total_power++;
        if(diff) {
            int idiff = (int)diff;
            assert(abs(idiff) &lt;= d*numeric_limits&lt;decltype(d)&gt;::epsilon());
            if(verbose) {
                cnt_diff++;
                sum_diff += idiff;
                if(idiff &lt; min_diff) {
                    min_diff = idiff;
                }
                if(idiff &gt; max_diff) {
                    max_diff = idiff;
                }
            }
        }
    }
}
if(verbose) {
    cout &lt;&lt; &quot;Total power: &quot; &lt;&lt; total_power
        &lt;&lt; &quot; Total diff: &quot; &lt;&lt; cnt_diff
        &lt;&lt; &quot; Mean diff: &quot; &lt;&lt; (dtype_t)sum_diff/cnt_diff 
        &lt;&lt; &quot; Min: &quot; &lt;&lt; min_diff
        &lt;&lt; &quot; Max: &quot; &lt;&lt; max_diff
        &lt;&lt; '\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&lt;int&gt;();
test_ipow&lt;unsigned&gt;();
test_ipow&lt;int64_t&gt;();
test_ipow&lt;uint64_t&gt;();

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");.