2. 符号付き同士の飽和加算
3. 符号無しと符号付き同士の飽和加算
あまりに久しぶりなので自分でもどんな話題をしていたのか忘れてしまってる...
最近とあるソースをみて思ったことをだらだら書いてみます(チェックをしっかりしてないので中身はいい加減です).
float func1(float a, float *buf, int n){
int i;
float x = 0;
for(i = 0; i < n; i++){
x += buf[i] / a;
}
return x;
}
コンパイラの最適化を期待しているのかも知れませんが、普通毎回割算してくれます. なんだかなあ, です. 当たり前ですが
float func1a(float a, float *buf, int n){
int i;
float x = 0;
for(i = 0; i < n; i++){
x += buf[i];
}
return x / a;
}
とかくべきですよね.
float func2(float a, float *buf, int n){
int i;
float x = 1;
for(i = 0; i < n; i++){
x *= pow(buf[i], a);
}
return x;
}
こんなのもありました. 思わず天を仰ぎたくなります. func1()と同様
float func2a(float a, float *buf, int n){
int i;
float x = 1;
for(i = 0; i < n; i++){
x *= buf[i];
}
return pow(x, a);
}
とすべきでしょう. ...ところがうまくいきませんでした.
調べてみるとbuf[]の値が1010程度の大きい数でxがoverflowしてしまっていたのです.
aは0.1程度の小さい数だったので個別にpowしてる時はoverflowしなかったのです.
(注)実際のコードではdouble xでもはみでてました.
Intel系CPUのfloatの最大値は約1038, doubleの最大値は10308です.
さてどうしましょうか.
float func2b(float a, float *buf, int n){
int i;
float x = 1, y;
int r, q;
q = n >> 3; /* 8 で割った商 */
r = n & 7; /* 8 で割った余り */
for(i = 0; i < q; i++){
int j;
y = 1;
for(j = 0; j < 8; j++) y *= buf[i];
x *= pow(y, a);
}
if(r){
y = 1;
for(i = 0; i < r; i++) y *= buf[i];
x *= pow(y, a);
}
return x;
}
こんな風にoverflowしない範囲に分割してやる方法があります. でも分割数が4にしないで大丈夫かどうか不安だし, SIMD化しにくいです.
No. 003の方法でxをdoubleにし, 分割を16にすればpowのコストは1/16になりますからそれでもよいかもしれませんが, ここでは別の方法を考えてみます. ちなみにPIIIでの演算コストはだいたい以下のようです.
| 命令 | クロック数 |
|---|---|
| 足し算/かけ算 | 2程度 |
| 割り算 | 37 |
| log | 82 |
| exp | 142 |
| pow | 366 |
8回に1度powを使うのは毎回割算をする程度のコストになります. かなり大きいですね. logのコストが割と小さいのでxy=exp(y * log(x))ですから, これを使ってみます.
float func2c(float a, float *buf, int n){
int i;
float x = 0;
for(i = 0; i < n; i++){
x += log(buf[i]);
}
return exp(a * x);
}
これでpowの時の約4～5倍の速度になります. 次はこのlogの計算をなくすことを考えてみます.
その前にfunc2aでoverflowしてしまうのなら最初からある程度の数を掛けて調節するのはどうでしょうか.
#define ADJ (1.0e-10)
float func2d(float a, float *buf, int n){
int i;
double x, y, z;
x = 1;
for(i = 0; i < n; i++){
x *= (buf[i] * ADJ);
}
#if 0
/* もしn, aが固定ならz(=pow(ADJ, -n*a)は先に計算しておける */
#else
y = pow(ADJ, -a);
z = 1;
while(n){
if(n & 1) z *= y;
y *= y;
n >>= 1;
}
#endif
return pow(x, a) * z;
}
buf[i]の値が常に1010程度と分かっているならこの方法が一番簡単でかつ速いと思います. ただ私が見たケースではbuf[i]の取り得る範囲が105〜1015と広く, またn=576などの大きめの数字だったので func2b()と組み合わせないといけませんでした. なにより決め打ち数値が必要と言うのがかっこ悪い(^^;です.
浮動小数の形式のところで触れましたが, 浮動小数は指数部と仮数部に分かれています.
これを使ってみます.
x=2 s * tの場合,
(s:指数部, t:仮数部)
log 2x = s +
log 2t です.
今求めたいものはΣ logx iですからこれを変形すると,
(底の2は省略)
Σ log x i
=Σ s i
+Σ log ti
=Σ s i
+log Π ti
よって指数部を足していき, 仮数部を掛けていけば良いのです. 仮数部は1以上2未満ですからoverflowはあり得ません. 実際にはこんなコードになります.
/* func2cのlogを分解 */
float func2ca(float a, float *buf, int n){
int i;
int sum_s;
float mul_t, x;
sum_s = 0;
mul_t = 1;
for(i = 0; i < n; i++){
int s, *p;
volatile float t = buf[i]; /* buf[i]を破壊するのでコピー */
p = (int *)&t;
s = (*p >> 23) - 127; /* 指数部 */
*p = *p & ((1<<23)-1) | (127<<23);/* 仮数部 */
sum_s += s;
mul_t *= t;
}
x = sum_s * log(2) + log(mul_t); /* 底の変換 */
return exp(a * x);
}
指数部で127を引いているところは最後にまとめてsum_sから127*nを引けば計算を減らせます. shift 1回, 論理演算2回, 足し算 1回, かけ算 1回と読み書きです. これだとSIMDによる並列化も難しくないからそこそこ高速に出来ると思います. もし途中に0が入った場合, sum_sが極端に小さくなる(0の時は指数部は-127)のでそれには注意しておいた方がいいです.
この指数部と仮数部に分ける方法は巾指数が決まっているpowの近似計算にも使えそうです. x=2s*tのとき pow(x,a) を求めるにはpow(2s,a)とpow(t,a)をそれぞれ求めて掛ければよいのです. sの方はせいぜい256通りのテーブルを作ればよいので簡単です. aの方は1以上2未満の浮動小数なのでテーブルと一次補間でそこそこ精度が出せます.
#define BITLEN 10
#define SIZE (1<<BITLEN)
const float a = 0.1;
static float idxTbl[256], powTbl[SIZE];
/* 第0次近似によるpow(x, a)の計算 */
float pow_fast(float x){
static int initFlag = 1;
float t;
int s, *p;
if(initFlag){
int i;
for(i = 0; i < 256; i++){
idxTbl[i] = pow(2, (i - 127) * a);
}
for(i = 0; i < SIZE; i++){
powTbl[i] = pow(1 + i / (double)SIZE, a);
}
initFlag = 0;
}
t = x;
p = (int *)&t;
s = (*p >> 23); /* 指数部 */
*p = *p & ((1<<23)-1) | (127<<23);/* 仮数部 */
return idxTbl[s] * powTbl[(int)((t - 1) * SIZE)];
}
上の最後の1行,
powTbl[(int)((t - 1) * SIZE)];
は
powTbl[(int)(t * SIZE) - SIZE];
とすればpowTbl - SIZEは前もって計算できるので引き算が1つ減ります(^^). ついでにfloat 型でSIZE (=1<<BITLEN) 倍というのは 指数の部分にBITLENを出せば実現出来るのでかけ算も減らせます. すなわち
*p = *p & ((1<<23)-1) | ((127+BITLEN)<<23);/* 仮数部 */
return idxTbl[s] * powTbl[(int)(t) - SIZE)];
と等価です. 実はもっと簡単になることを昨日気がつきました. その話題にいく前に pow_fast() の中で一番時間がかかっている部分はどこか分かりますか(初期化部分というのは無し)?
実はなんとintへのcastの部分です. VC++で約78clk, gccでも50clk以上かかります. No. 004を見ても論外な遅さです. 何故こんなことになってるのでしょうか. これはVC++のintへのcast部分のルーチンです.
__ftol:
push ebp
mov ebp,esp
add esp,0FFFFFFF4h
wait
fnstcw word ptr [ebp-2]
wait
mov ax,word ptr [ebp-2]
or ah,0Ch
mov word ptr [ebp-4],ax
fldcw word ptr [ebp-4]
fistp qword ptr [ebp-0Ch]
fldcw word ptr [ebp-2]
mov eax,dword ptr [ebp-0Ch]
mov edx,dword ptr [ebp-8]
leave
ret
waitやleave命令など386, 486時代の名残り(?)の命令が見られます. それから殆んど使われないfnstw, fldcw. これらはFPUの状態を変更するための命令です. fist命令の結果を得る時は切り上げ, 切り下げ, 四捨五入を指定できるのですが, 制御レジスタはユーザが自由にいじれるのでどんな状態になっていても正しく動作するように毎回設定し直してあるのです.
厳密にするため仕方無いとしても, その様な命令の処理時間が大きいのは想像に難くありません. しかも普通はデフォルトの四捨五入モードから変更することは無いので無駄な処理です. 浮動小数が0以上の場合は
_f05: 0.5
fld dword [x]
fsub dword [_f05] ;0.5を引いて
fistp dword [int_x] ;四捨五入=>切捨て
で十分であり遥かに高速です. データの量子化や画像処理におけるドット打ちのときに知らないうちに使ってそうで結構恐いです. #ぱっと見ててもなんで遅いのか分かりにくい.
計さんからご指摘がありましたが, 完全に私は勘違いしてました.m(__)m
compilerが毎回0方向へ丸めるモードに移行していたのは
四捨五入モードでの `.5'の扱いがいつも切り上げであるのではないためでした.
具体的には
| 0 ← 0.5 | 1.5 → 2 ← 2.5 | 3.5 → 4 ← 4.5 | 5.5 → 6 ← 6.5 | ... |
|---|
というようになってます. これは.5でいつも切り上げだと平均して大きくなる傾向が出てしまうためです. だから常に切り上げだと思って書いた No.009 のコードは間違ってます. ただdskさんが書かれたようにループ内で何度も変更するのは無駄なので どうしても必要な部分の前後で切り換えるようにすべきです.
コントロールレジスタ(cw)の10, 11bitが丸め制御フィールド(RC)です. 値の意味は次の通り
| RC | 意味 |
|---|---|
| 00 | 四捨五入(但し.5の時は偶数方向に丸める) |
| 01 | 負の無限大に切り捨て |
| 10 | 正の無限大に切り上げ |
| 11 | 0方向へ切り捨て |
;int loadCW(void);
proc loadCW
push eax
fstcw word [esp]
pop eax
ret
;void changeCW(int cw);
proc changeCW
mov eax,[esp+4]
push eax
fldcw word [esp]
pop eax
ret
こんな関数を用意して
int cw;
cw = loadCW();
changeCW(cw | (3<<10)); /* 0方向へ切り捨て */
/* intへのcastにfistを使った演算ループ */
changeCW(cw); /* 元のモードに戻す */
の様にすればよいことになります.
余談ですが, 私は小さい時バスや電車の子供運賃が
「大人の半額. 但し10円未満の端数は四捨五入」と書いてあるのを見るたび,
「いつも切り上げになるんだからなんかずるい」
と思っていたのですが,
こんな四捨五入モードだったら損した気分にならずにすんだかもしれません:)
No. 007では仮数部の入力値が1以上2未満の小数なので, テーブル引きをするためSIZE倍とcastをしていたのでした. しかしよく考えたら仮数部自体が既に23bitの固定小数の形式とみなせます. 従って (指数部)→floatレジスタへ読み込み→intへのcastというのは2度手間です. 結局次の様にすれば良かったのです.
/* pow_fastの改良(テーブル初期化部分は除く) */
float pow_fast2(float x){
int s, t;
s = t = *(int *)&x;
s >>= 23; /* 指数部 */
t = (t & ((1<<23)-1)) >> (23 - BITLEN); /* 仮数部 */
return idxTbl[s] * powTbl[t];
}
次は具体的にsqrtの近似計算ルーチンを作ってみます.
流石に補間なしでは心許ないので一次補間をやってみます.
/* 一次補間ありの平方根 */
#define BITLEN 8
#define SIZE (1<<BITLEN)
static float idxTbl[256];
static float sqrtTbl1[SIZE];
static float sqrtTbl2[SIZE];
init_sqrt_fast(void){
int i;
for(i = 0; i < 256; i++){
idxTbl[i] = (float)pow(2, (i - 127) * 0.5);
}
for(i = 0; i < SIZE; i++){
sqrtTbl1[i] = (float)sqrt(1 + i / (double)SIZE);
}
for(i = 0; i < SIZE-1; i++){
sqrtTbl2[i] = (sqrtTbl1[i+1] - sqrtTbl1[i]) / (float)SIZE;
}
sqrtTbl2[SIZE-1] = ((float)sqrt(2) - sqrtTbl1[SIZE-1]) / (float)SIZE;
}
static float sqrt_fast(float x){
#if 0 /* 補完無し */
int s, t;
s = t = *(int *)&x;
s >>= 23;
t &= ((1<<23)-1);
t >>= (23-BITLEN);
return idxTbl[s] * sqrtTbl1[t];
#else /* 一次補完 */
int s, t, u;
s = t = u = *(int *)&x;
s >>= 23;
t &= ((1<<23)-1);
u &= ((1<<(23-BITLEN))-1);
t >>= (23-BITLEN);
u >>= (23-BITLEN*2);
return idxTbl[s] * (sqrtTbl1[t] + sqrtTbl2[t] * u);
#endif
}
素朴に書いてみました. sqrtTbl1[i]はその点での平方根そのものの値, sqrtTbl2[i]は(i,sqrtTbl1[i])と(i+1,sqrtTbl1[i+1])を結んだ線分の傾きです. 仮数部の上位BITLEN bitが補完無しの場合の近似, 次のBITLEN bitが補正項(u)です. 実際にどの程度のパフォーマンスなのか調べてみました.
ついでに他の命令も併せて表記しました(実測値なのでメモリ状態やCコンパイラのループ展開などで数値は変動します).
| 命令 | クロック数 | ||
|---|---|---|---|
| - | PIIIE | Athlon | K6-III(参考値) |
| 足し算/かけ算 | 2〜4 | 3 | 7 |
| 割り算 | 32 | 18 | 44 |
| Cの(int)cast | 60 | 38 | 71 |
| fistp使用cast | 17 | 16 | 33 |
| 平方根 | 58 | 24 | 48 |
| 今回の近似平方根 | 24 | 18 | 55 |
| sin | 98 | 88 | 79 |
| log | 82 | 123 | 150 |
| exp | 142 | 62 | 137 |
| pow | 300 | 300 | 450 |
一番速くなって欲しいK6-IIIでsqrt_fastが遅くなってるやん(;;. これは多分TOWNSに刺さってるからで普通のAT互換機上だともっと良い結果がでるはず. sqrt_powも速くなってるだろう(と信じたい). ただ3D Now!やSSEを使うと平方根は数クロックなんでね. 本当に速度が必要ならそれらを使うべきでしょう.
一次補完するテーブルの係数をいじって少し精度をあげられないかテストしてみました. 区間[a, b]上でf(x)=√xとg(x)=px+qの差の巾乗が最小になるようにpとqを決めるのです. つまりD=∫ab (f-g)2dxと置いて∂pD=∂qD=0となる(p,q)を求めます. いわゆる最小二乗法というやつです. 知人がMapleを持っていたのが借りてちょこっと計算してみると結果は
p = (4/5)(a2.5 + 5ab1.5 - 5ba1.5 - b2.5) / (a-b)3
q = (4/15)(-a3.5 + 10a2b1.5 - 10a1.5b2 + ab2.5 + b3.5) / (a-b)3
一瞬で答えがでてしまった. 思わず欲しくなりました(^^;. 私も買おうかな. それはさておき No. 012のinit_sqrt_fast()を次の様に修正して実験してみました. ついでにsqrtTbl1とsqrtTbl2を分けるとテーブル引きのときにメモリアクセスが分断されるのでくっつけて配置してます.
static float sqrtTbl[SIZE*2]; /* Tbl1とTbl2をインタリーブ */
init_sqrt_fast(void){
int i;
for(i = 0; i < 256; i++){
idxTbl[i] = (float)pow(2, (i - 127) * 0.5);
}
for(i = 0; i < SIZE; i++){
double a, b, sa, sb, t;
a = 1 + i / (double)SIZE;
b = a + 1 / (double)SIZE;
sa = sqrt(a);
sb = sqrt(b);
t = (sa + sb) * (sa + sb);
sqrtTbl[i*2] = (float)((b*sb + 2 * b*sa + 8 * a*sb + 4 * a*sa) * 4 / (15 * t));
t *= (sa + sb);
sqrtTbl[i*2+1] = (float)((a + 3 * sa*sb + b) * 4 / (5 * t * SIZE * SIZE));
}
}
で, 実際にテーブルを換えた事で誤差がどのぐらい変わったかというと... 残念ながら5%程度しかよくなりませんでした(^^;
sqrt_fastでせっかく最小二乗法を使ったのに効果がなかったので別のことに応用してみます. かなり前ですが 3D Now!命令を使ったlogルーチンに対して, もっとよい方法があるとのコメントがありました( BBS No.174).
つまり|z|≦3-2√2のとき
log((1+z) / (1-z)) = 2 * (z + z3/3 + z5/5 + z7/7)
で近似計算をしていたのですが,
log((1+z) / (1-z)) = 2 * (az + bz3 + cz5)
但し
a = 1.0+1.42606e-06 b = 1/3-0.000262124 c = 1/5+0.0117245
を使った方がよいということでした. 実際に
double log_fast(double x, double a, double b, double c){
double z, zz;
z = (x - sqrt(2)) / (x + sqrt(2));
zz = z*z;
return 0.5 * log(2) + 2 *z * (a + zz * (b+c*zz));
}
void calc_diff(double a, double b, double c){
double diff, x, s, t;
const int N = 1000;
int i;
diff = 0;
for(i = 0; i < N; i++){
x = i / (double)(N+1) + 1;
s = log(x);
t = log_fast(x, a, b, c);
diff += fabs(s - t);
}
printf("coef a=1%+e, b=1/3%+e, c=1/5%+e\n", a-1, b-1.0/3, c-1.0/5);
printf("ave. diff=%e\n", diff / N);
}
という実験用コードを使い[1,2]の区間での誤差を調べると確かによいのです. ただその時はどうやって求めたのか考える気力がなかったのでそのままほっておきました. で, これに最小二乗法を使ってみようと思い立ったのです. 結果は以下の通り
| - | 2項展開(5次まで) | 2項展開(7次まで) | T氏のやりかた | 今回の最小二乗法 |
|---|---|---|---|---|
| 1次の係数(a) - 1 | 0 | 0 | 1.426060e-06 | 3.104605e-007 |
| 3次の係数(b) - 1/3 | 0 | 0 | -2.621240e-04 | -9.440747e-005 |
| 5次の係数(c) - 1/5 | 0 | 0 | 1.172450e-02 | 6.987293e-003 |
| 7次の係数 - 1/7 | - | 0 | - | - |
| 平均誤差 | 1.648350e-07 | 3.027844e-09 | 5.465091e-08 | 1.089739e-08 |
| 最大誤差 | 1.279805e-06 | 2.933744e-08 | 9.203439e-08 | 4.923651e-08 |
なかなか良い結果を出しています. 同じ5次までしか計算しないのに2項展開に比べて1桁良くなってます. 正直こんなに変わるとは思いませんでした.
ではどうやって(a, b, c)を求めたのか? 少々難しい(力のある高校生か大学1回生程度の数学)ですか解説してみます. やりたいことは区間[-S, S] (S = 3 - 2√2)において
f(z) = log((1+z) / (1-z)) と g(z) = az + bz3 + cz5
の誤差が最小になるような(a, b, c)を求めることです. 誤差は|f(z) - g(z)|です. それを[-S, S]でかき集めればよいのですが, 絶対値は解析的に扱いにくいので代わりに二乗するのです.
かき集めるのは積分です. つまり
D = ∫-SS (f(z) - g(z))2dz
を求めて, そのDが最小となる(a, b, c)を求めるのです.
前回は知人にMapleでやってもらったのですが, 今回は手計算で行ないました.
#1回生の頃に比べて思いっきり計算力落ちてました. いや〜疲れた(^^;
全部展開してやるとz f(z), z3 f(z), z5 f(z)の積分が出てきますが, これは高校生でも出来ます. 積分結果をfa(z), fb(z), fc(z)と名前をつけておきます.
Sを固定すれば D は a, b, cそれぞれについて2次式です. 最小値を求めたいのでa, b, cについて微分します.
∂a D = ∂b D = ∂c D = 0
∂は偏微分といいます. ∂aは b, cを固定して a で微分せよという意味です. もともと2次式なので微分すれば1次式です. したがって3元1次連立方程式を解けばOK.
(1/3 1/5 1/7)
A=(1/5 1/7 1/9)
(1/7 1/9 1/11)
という行列なのでその逆行列を求めるのが多少楽です. こんな関係式になりました.
2 ( a bS2 cS4 ) A = (fa/S3 fb/S5 fc/S7)
A-1を求めて最後に S = 3 - 2√2をぶちこめば出来上がります. プログラムは こんな感じ.
割算が高速なSIMD利用を前提とした実際のlogのCなコードはこんな感じです.
float logC( float x )
{
float a, z, zz;
unsigned int *p;
int b;
const float A = 1.0 + 3.104605e-7; /* A*2 = 2.0000006209 */
const float B = 1.0/3 - 9.440747e-5; /* B*2 = 0.6664778517 */
const float C = 1.0/5 + 6.987293e-3; /* C*2 = 0.4139745860 */
p = (unsigned int *)&x;
b = (*p >> 23) - 127; /* suppose x_org > 0 */
*p= ( *p & 0x007FFFFF ) | 0x3F800000;
a = *(float *)p;
z = (a - SQRT2) / (a + SQRT2);
zz = z * z;
z = (b+0.5) * LOG_E_2 + z * ((A*2) + zz * ((B*2) + zz * (C*2)));
return z;
}
ただ FPU <=> 整数 のやりとりが必要なためなかなか使いにくいです. レジスタを共有している3D Now!系では問題無いのですがそうでないSSEでは絶対メモリを介在する必要があるためロスが大きいのです. しかも4つ並列にしないとおいしくないし. SSE2ではXMMレジスタ<=>MMXレジスタ間の命令が追加されましたがこれはSSEの段階でつけるべきでした. まったくIntelってのは... ところでmp3エンコーダでは音響心理解析の部分でかなり細かい演算をえんえんとしています. そのなかにはlogやpowやsinなど超越関数もかなりあってそれらをいかに減らすかがポイントになります.
心理解析部分(){
foreach(チャンネル[右, 左, まんなか, 両端]){
ごちゃごちゃっと300～400行;
}
}
gogo2のときはあっさりあきらめて手をつけていなかったのですがpetitでは気合いをいれて高速化してます. 具体的には上の各チャンネルごとに処理していた部分をデータ自体をインタリーブしてベクトル的に処理することにしたのです.
このプログラムはどちらが速いでしょうか.
#define SIZE 128
static void test1(const int src[SIZE][SIZE], int dst[SIZE][SIZE])
{
int i, j;
for (i = 0; i < SIZE; i++) {
for (j = 0; j < SIZE; j++) {
dst[j][i] = abs(src[j][i]);
}
}
}
static void test2(const int src[SIZE][SIZE], int dst[SIZE][SIZE])
{
int i, j;
for (i = 0; i < SIZE; i++) {
for (j = 0; j < SIZE; j++) {
dst[i][j] = abs(src[i][j]);
}
}
}
int main(void)
{
int src[SIZE][SIZE];
int dst[SIZE][SIZE];
test1(src, dst);
test2(src, dst);
return 0;
}
無論, 「連続アクセスしている test2 が速い」が答えなのですが, P4 でテストしたところなんと test1 の方が速かったのです. あれれ? と追試したのですが何度しても同じ. ちなみに P3 では test1 が 速いというまともな答えでした.
不思議に思って SIZE の大きさを変えてテストしてみると
| SIZE | test1(Kclk) | test2(Kclk) |
|---|---|---|
| 63 | 54.3 | 20.4 |
| 64 | 61.8 | 31.1 |
| 65 | 54.9 | 21.6 |
| 113 | 174 | 63.2 |
| 127 | 226 | 79.2 |
| 128 | 252 | 333 |
| 129 | 237 | 80.7 |
と SIZE = 128 のときだけ異常に遅い. ここまで来てはっと気が付いて src だけ static 変数にしてみたら果たして 80Kclk 程度と納得いく数値になりました. 種明かしをすると SIZE = 128 のときは src の大きさは 128 * 128 * 4 = 65Kbyte になります.
P4のマニュアルを見直すと
実行中の参照 (load または store) と次の参照 (load または store) のアドレスの下位 16bit が一致する場合は処理速度が1/3になる.
とかいてありました. まさにこれにひっかかっていたのです.
他には 2K, 16K, 32K である条件が満たされると同様のことが発生するようです. alignment をそろえるために テーブルサイズを 2 のべきにすることが多いので 気をつけないといけないですね.
P3 や Athlon でもそんな制約はあったような気がしますがここまで露骨な現象にあったことが無かったので殆ど気にしてませんでした. やっぱりマニュアルはしっかり読まんといかんのか(^^;
typedef unsigned char uint8;
typedef signed char sint8;
void add_us_s_and_clip(uint8 *p_us, sint8 *p_s, int count)
{
int i;
for (i = 0; i < count; i++) {
int a, b;
a = p_us[i];
b = p_s[i];
a += b;
if (a < 0) {
a = 0;
} else if (a > 255) {
a = 255;
}
p_us[i] = (uint8)a;
}
}
というコードの MMX 化を考えてみました. paddusb は符号無し同士の飽和加算だし, paddsb は符号付き同士の飽和加算なので使えません. 最初は us_a, s_b をともに word に符号拡張して paddsw をと思いましたが, 符号付き符号拡張は面倒だし後で byte 変換するのももどかしい. 0 ≦ a ≦ 255, -128 ≦ b ≦ 127 なので
の4パターンに分けて mask1 = (2番目が真) ? 255 : 0, mask2 = (4番目が真) ? 0 : 255 となるマスクを生成して (a + b) | mask1 & mask2 というふうに考えてみました. ...が泥沼(^^;
発想を変えるために当初使えないと思った paddusb と paddsb のグラフを考えてみることにしました. こんな感じ
色がついているところがクリッピングされているところでそれ以外は (a+b)&0xFF です.
これを見て最初はクリッピングパターンが違うのでどうしようもないと思ったのですが,
よく見ると2. と 3. は縦に2つに切って交換してつなぎあわせれば同じ形になります.
さてさて SJIS判定 と同じではないか, ということで左右をひっくり返すために
^ 0x80 をすればよい. すなわち
add_us_s_and_clip(a, b) = paddsb(a^0x80, b) ^ 0x80
なることが分かりました(^^)/.
%rep 20
%if 0
pxor mm0, mm0
paddq mm0, [esi]
movq mm1, mm0
pand mm1, mm7
movq mm1, mm0 ; [H:L]
psrad mm1, 31 ; [sign H:*]
punpckhdq mm0, mm1 ; [sign H:H]
%else
movq mm0, [esi]
movq mm1, mm0
pand mm1, mm7
movq mm1, mm0 ; [H:L]
psrad mm1, 31 ; [sign H:*]
punpckhdq mm0, mm1 ; [sign H:H]
%endif
%endrep
このコードで %if 0 と %if 1 とでどちらが速いのか測定してみました. 話の発端はPentium4 での MMX. 私のところでは P4 Willamette (stepping 15/0/7) 1.5GHz, mobile P4 (stepping 15/2/4) 1.7GHz のどちらでも上が下よりも 40clock つまり1単位あたり 2clock 遅い, つまり妥当な結果でした. 単純に考えるならこの2つの P4 の世代のどこかでバグっていたって感じでしょうか. それ以前にそもそも P4 での movq や movaps のレイテンシが恐ろしく遅いことからしてよく分からんのですが. その後の追試によると勘違いだったらしい.
発端はサポート情報のp.154なんですが, こんなことを考えてみました.
「r1, r2をランダムな53ビットの非負な整数として double な変数 x, y に対し
x = r1/253, y = r2/253 としたとき,
x2 + y2 == 1 となる (x, y) はどれぐらいあるのか?」
#x, y の取り方にそもそも作為的な条件があるのはご勘弁.
0 < y < x で考えるとして x が1に近づくと√(1 - x2) は小さくなるので y の有効桁数が少なり等式が成り立つ数が増えます. Count.javaで数えてみました. グラフにするとこんな感じ(横軸は i のとき x = 1-e*1016-i をあらわす. ただし e = 2-53. 分かりにくくてすいません).
見事に直線に乗ってます. となるとこの数値(0.5352など)は何? というのが興味の対象になりますが時間が無いのでまたいずれ.