解決済みですが...
目的はbitmap形式向けDCT後の、BGR->YCbCr変換
ご質問の処理は
・得たい値は浮動小数点数
・限界まで高速にしたい
・処理系がIEEE754を用いていること前提
・ビットパターンをバリバリに意識する
・この丸め論理の定義域・値域に特定の前提を置ける(NaNであったり負の値ではないなど)
といった前提を置きポータビリティーはある程度犠牲でいいぐらいの勢いなのかなと想像しました。
浮動小数点数演算がソフトウェアで実装されていた時代を経験している自分はつい「浮動小数点数演算は遅い」という意識がぬぐい切れず今だに整数演算でこちょこちょやった方が早いかもなんて考えが浮かびます。そして大抵の場合ハードウェア(プロセッサーに搭載されている浮動小数点演算ユニット的なもの)に素直に任せた方がずっと早いことを再認識するハメになります・・・orz
でもおもしろそうだったので実際に比較してみました。
C++
1 # include <cassert>
2 # include <cstdint>
3 # include <cinttypes>
4 # include <cmath>
5 # include <chrono>
6 # include <cstdio>
7
8 double f1 ( double y ) {
9 auto iy = static_cast < std :: int64_t > ( y ) ;
10 return ( iy + 128 ) >> 8 ;
11 }
12
13 double f2 ( double y ) {
14 int exp ;
15 double normalized = std :: frexp ( y + 128 , & exp ) ;
16 double y2 = std :: ldexp ( normalized , exp - 8 ) ;
17 return floor ( y2 ) ;
18 }
19
20 double f3 ( double y ) {
21 int exp ;
22 double normalized = std :: frexp ( y , & exp ) ;
23 double y2 = std :: ldexp ( normalized , exp - 8 ) ;
24 return y2 ;
25 }
26
27 double f4 ( double y ) {
28 assert ( ! std :: isnan ( y ) && ! std :: isinf ( y ) && y >= 0 ) ;
29
30 auto iy = * reinterpret_cast < std :: int64_t * > ( & y ) ;
31 // 8bit右シフト(256で除算)した結果の指数部を求める
32 // 非負を仮定しているので符号ビットのマスクは省略
33 auto mantissa = iy & INT64_C ( 0x000F 'FFFF' FFFF'FFFF ) ;
34 auto exp = ( int ) ( iy >> 52 ) - 8 ;
35 if ( exp <= 0x3ff - 1 ) {
36 // 値が0.5以下なら結果は0か1
37 return exp == 0x3ff - 1 ? 1.0 : 0.0 ;
38 } else if ( exp < 0x3ff + 52 ) {
39 // 仮数分に小数以下の値を含む
40 // exp 1.0の桁
41 // -------- ---------------
42 // 0x3ff 0x0010_..._0000
43 // 0x3ff+1 0x0008_..._0000
44 // 0x3ff+52 0x0000_..._0001
45 auto pos = 0x3ff + 51 - exp ; // 0.5のビット位置
46 auto one_half = INT64_C ( 1 ) << pos ;
47 mantissa += one_half ;
48 mantissa &= ~ ( ( one_half << 1 ) - 1 ) ; // 小数部のビットを落とす
49 if ( mantissa & INT64_C ( 0x0010 '0000' 0000 ' 0000 ) ) {
50 // 桁上がりしたら仮数部と指数部を補正
51 mantissa = ( mantissa & INT64_C ( 0x000F 'FFFF' FFFF'FFFF ) ) >> 1 ;
52 exp ++ ;
53 }
54 }
55 double result ;
56 * reinterpret_cast < int64_t * > ( & result ) = static_cast < uint64_t > ( exp ) << 52
57 | mantissa ;
58 return result ;
59 }
60
61 # define LOOP_COUNT 10000000
62 volatile double z ;
63
64 double timeit ( double ( * fp ) ( double ) , double arg ) {
65 # if ! defined ( NDEBUG )
66 assert ( false ) ; // デバッグモードで性能測定しないように...
67 # endif
68 auto t0 = std :: chrono :: system_clock :: now ( ) ;
69 for ( int i = 0 ; i < LOOP_COUNT ; i ++ ) {
70 z = fp ( arg ) ;
71 }
72 auto t1 = std :: chrono :: system_clock :: now ( ) ;
73 auto elapsed = std :: chrono :: duration_cast < std :: chrono :: nanoseconds > ( t1 - t0 ) . count ( ) ;
74 return elapsed / LOOP_COUNT ;
75 }
76
77 using namespace std ;
78
79 int main ( ) {
80 double ys [ ] = {
81 256 * 0.499 , 256 * 0.5 ,
82 ( INT64_C ( 0x000F 'FFFF' FFFF'FFFF ) ) * 256.0 ,
83 ( INT64_C ( 0x000F 'FFFF' FFFF'FFFF ) + 0.5 ) * 256.0 ,
84 } ;
85 for ( int i = 0 ; i < sizeof ( ys ) / sizeof ( ys [ 0 ] ) ; i ++ ) {
86 double y = ys [ i ] ;
87 printf ( "------------------------------------\n" ) ;
88 printf ( "y = %20lf\n" , y ) ;
89 printf ( " f1 -> %20lf %3.0lf ns\n" , f1 ( y ) , timeit ( f1 , y ) ) ;
90 printf ( " f2 -> %20lf %3.0lf ns\n" , f2 ( y ) , timeit ( f2 , y ) ) ;
91 printf ( " f3 -> %20lf %3.0lf ns\n" , f3 ( y ) , timeit ( f3 , y ) ) ;
92 printf ( " f4 -> %20lf %3.0lf ns\n" , f4 ( y ) , timeit ( f4 , y ) ) ;
93 }
94 return 0 ;
95 }
本件で注意が必要なのは「四捨五入の論理」だと思います。また浮動小数点数の演算を避けるためにfrexp, ldexpのようなものを用いたとき「どのくらいオーバーヘッドがあるか」を知っておくことも無駄ではないでしょう。
上記では4つぐらい方法を試してみてます。
f1は浮動小数点数を整数化して整数演算として四捨五入&256での除算を行う方法。
f2はfrexp/ldexpと浮動小数点数の加算、floorによる切り捨てを行ったもの。
f3は四捨五入の論理を放棄して単にfrexp/ldexpにより256で割るだけのもの。
f4はポータビリティー無視でIEEE754のフォーマットを直接意識したもの。
自分のPC(Intel core i5-8400 2.8GHz)を用いて実験してみますと
bash
1 $ g++ -DNDEBUG t.cpp
2 $ ./a.out
3 ------------------------------------
4 y = 127.744000
5 f1 - > 0.000000 2 ns
6 f2 - > 0.000000 17 ns
7 f3 - > 0.499000 13 ns
8 f4 - > 0.000000 3 ns
9 ------------------------------------
10 y = 128.000000
11 f1 - > 1.000000 2 ns
12 f2 - > 1.000000 18 ns
13 f3 - > 0.500000 13 ns
14 f4 - > 1.000000 3 ns
15 ------------------------------------
16 y = 1152921504606846720.000000
17 f1 - > 4503599627370495.000000 2 ns
18 f2 - > 4503599627370495.000000 20 ns
19 f3 - > 4503599627370495.000000 13 ns
20 f4 - > 4503599627370495.000000 7 ns
21 ------------------------------------
22 y = 1152921504606846848.000000
23 f1 - > 4503599627370496.000000 2 ns
24 f2 - > 4503599627370496.000000 17 ns
25 f3 - > 4503599627370495.500000 13 ns
26 f4 - > 4503599627370496.000000 8 ns
結論を言えば平易な論理(f1)が一番早かったです。
frexp/ldexpのアセンブリコードを見ますとfrexp/ldexpともに関数呼び出しになっていました。f1とf3を比較しますと、浮動小数点数<->整数変換より関数呼び出し(分岐の)オーバーヘッドの方が大きいということなのかも知れません。またfrexp/ldexpの関数呼び出しオーバーヘッドを避ける意図で書いた論理(f4)は、整数演算しかしていないとはいえ、そこそこの命令数ですし分岐も含まれているせいか平易な論理(f1)にはかないませんでした。
frexp/ldexpを用い頑張って論理を書こうとするより、やはりmaisumakunさん回答のアドバイス
単に「整数にキャストしてからシフト」だけでいい
が得策のようです。f1では整数型として得た結果を再度浮動小数点数にキャストしてますが、それでもハードウェアにまかせた方が早かったので。
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
2019/06/23 14:07
2019/06/23 15:18