8
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

整数乗を自動判別するpowを作った

Last updated at Posted at 2018-09-03

#動機
数値計算でxの二乗などをするとき、よく#define SQUARE(x) ((x)*(x))などとし、std::powは(遅いので)使うなと言われる[要出典]
それもあるのだが、単に整数の整数乗(3^5など)を計算したいときも、いちいちdoubleにキャストしてpowして戻り値をintにキャストするのが馬鹿らしい。(c++に累乗演算子などなかった)

#コード(没)
pow(x,n)で、x,nがコンパイル時に分かっている場合にはconstexpr演算するのはもちろん、nが整数の場合には再帰的に掛け算してxの型のまま計算するコードを書いた。(なお、c++初心者でSFINAEとかc++11 constexprとか扱えないのでif constexprを使ったため、c++17以降でないと使えません)
hnn::pow(x,n)は、nの型が整数かによって自動的にふさわしいpowを呼んでくれます。

botu.cpp
#include <iostream>
#include <type_traits>

#include <cmath>

#if __has_include("sprout/math.hpp")
    #include "sprout/math.hpp"
    #define HAS_SPROUT_POW
#endif

namespace hnn{
template<class X,class N>
class pow{
    X val = 0;
public:
    constexpr pow(X x,N n){
        //Nが整数なら倍々イディオムでO(logN)で計算
        if constexpr (std::is_integral<N>::value){
            if (0 > n){
                val = 1 / hnn::pow(x, -n);
                return;
            }
	        val = 1;
	        while (0 < n) {
		      if ((n % 2) == 0) {
			     x *= x;
			     n /= 2;
		      }
		      else {
			     val *= x;
			     --n;
		      }
            }
        }
        else{
#ifdef HAS_SPROUT_POW
            val = sprout::pow(x,n);
#else
            val = std::pow(x,n);
#endif
        }
    }
    constexpr operator X() const{
        return val;
    }
};
}//namespace hnn{

int main(int argc, char* argv[]){
    constexpr auto a=hnn::pow(1.0000001,10000000);      //自然対数の近似値
    auto b=hnn::pow(2.0,3.1);
    constexpr auto c=hnn::pow(1.9,60);
    constexpr auto d=hnn::pow(1.9,60.0);
    using namespace std;
    cout<<a<<","<<b<<","<<c<<","<<d<<endl;
    if(argc>1){
        const int value = std::atoi(argv[1]);
        //int型のまま計算
        const int valueCubed=hnn::pow(value,3);
        cout<<"input cubed="<<valueCubed<<std::endl;
    }
    return 0;
}

#何が問題か
一見constexpr pow(X x,N n)のおかげでnがコンパイル時定数ならx*x*xのように展開してくれると思いましたが、そううまくは行きません。constexpr指定はx,n両方がコンパイル時定数の時にしか働かないので、xがコンパイル時に分からない場合はnが定数であってもアセンブラでのcall命令が残り、マクロの場合#define SQUARE(x) ((x)*(x))にパフォーマンスで負けてしまいます。nを再帰テンプレートで扱うときちんと展開されますが、今度はnがコンパイル時に分からない場合にコンパイルエラーになってしまいます・・・ぐぬぬ

#そのpowやばくないですか
xのn乗を計算するhnn::pow(x,n)で、

  1. nがコンパイル定数かつ整数ならx*x*x*x...に展開
  2. x,nが共にコンパイル定数ならsprout::powに投げる
  3. そうでないならstd::powに投げる
    このpow、機能多くないですか?そもそも、hnn::pow(x,n)の行を見た段階で、x,nの型を見に行かないと計算量やふるまいが分からないのは可読上良くありませんね。まるで邪悪なc形式キャストみたいだ

#顧客が本当に欲していたpow
元々の問題は#define SQUARE(x) ((x)*(x))みたいなマクロをインライン関数などに直したい!ってことでした。その際、
・xの型のままで計算する(std::powのように浮動小数へのキャストをしたくない)
・インライン関数でsquare(double x),cubic(double x)...などと並べるのはいかにも格好悪い(任意の整数乗に対応したい)
・パフォーマンスはマクロの場合を下回ってはならない
ということは、pow<N>(x)のNに整数をテンプレートで渡して、xはクラステンプレートと引数の型推論を使って任意の型に対応させればよいですね。

#コード(真)

test.cpp
#include <iostream>

namespace hnn{
template <class T,int N>
struct Impl{
    static constexpr T calc(T x){
        if constexpr(N < 0){
            return 1/x*Impl<T,-N>::calc(x);
        }
        else if constexpr(N&1){
            return x*Impl<T,N-1>::calc(x);
        }
        else if constexpr(N==0){
            (void)x;
            return 1.0;
        }
        else{
            return Impl<T,N/2>::calc(x*x);
        }
    }
};


template <int N,class T>
constexpr T pow(T x){
    return Impl<T,N>::calc(x);
}
}//namespace hnn{

int main(int argc, char* argv[]){
    //コンパイル時自然対数計算
    constexpr auto napier=hnn::pow<1234567>(1.0+1.0/1234567);
    //コンパイル時に分からない数の3乗
    const auto input = argc>1? std::atoi(argv[1]) :0;
    const auto cubic=hnn::pow<3>(input);
    //確認
    std::cout<<napier<<","<<cubic<<std::endl;
    return 0;
}

#結論
これで少なくてもgcc -O2ならちゃんとインライン展開してくれました。あとは既存コードでSQUARE(x)hnn::pow<2>(x)に変えるだけですね。c++に累乗演算子と小数剰余演算子はあっていいと思う

8
3
11

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
8
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?