多倍長整数でSteinのGCDを使う

次のような、32bit系の多倍長整数クラスがあったとして、


class Bigint
{
public:
	Bigint();
	...
	... 様々な定義(四則演算、operator=,operator<<=等)
	...
private:
	int   m_sz;   // 配列数
	int*  m_num;  // 配列(リトルエンディアン)
};

16進で次のように表せる整数があれば
0x111100002222000033330000

このクラスでは次のように格納されるものとする。(注1)

(m_sz=3)
|- m_num[0]-|- m_num[1]-|- m_num[2]-|
| 3333 0000 | 2222 0000 | 1111 0000 |


これを用いて、Stein's algorithmによるGCD計算を考える。


SteinのGCDについてはWikiPediaが詳しい。
http://en.wikipedia.org/wiki/Binary_GCD_algorithm

一般に ユークリッドの互除法(拡張) の方がループ回数は少ないが、
剰余計算がネックになり、応答性が悪くなる場合が多い。

但し、実装する乗算・剰余算のアルゴリズムや数値の桁数によっては、
ループ回数の少ないユークリッドの互除法(拡張)の方が早い場合もある。

Steinのアルゴリズムで遅いと思われるのは、Shift演算の回数に尽きる。
サンプルコードをそのまま使うと、いくら負荷の少ないShift演算と言えども、
多桁演算で処理回数が多くなると大きな遅延を招きます。
WikiPediaのSteinのサンプルコードで問題となるのは、


/* From here on, u is always odd. */
do {
  while ((v & 1) == 0)  /* Loop X */
     v >>= 1;


の部分。要は右端の0を除去しているだけである。
右端に並んでいる0の数を数える関数を用意して、Shift回数を減らことを目指す。

http://en.wikipedia.org/wiki/Count_trailing_zeros にもあるが、
intel系であれば、アセンブラ(bsf)が使える。


//
// 高速化のためゼロチェックを行っていない
//  ゼロだとアクセス違反の例外。
//
int Bigint::bsf() const{
	int* p = m_num;
	__asm{
		mov edx, [p]
		mov eax, edx
bsf_l2:
		bsf ecx, [eax]
		jnz bsf_l1
		add eax, 4
		jmp bsf_l2
bsf_l1:
		sub eax, edx
		shr eax, 2
		mov edx, BITPERBLOCK
		mul edx
		add eax, ecx
	}
}


BITPERBLOCKは、配列1つに含まれる多倍長整数のBit数である。
32bitであれば、32が使われる場合が多い。(注2)
(64bitの場合、bsfqを使う。)

(m_sz=2, BITPERBLOCK=32)の次のような多倍長整数があるとき

|- m_num[0]-|- m_num[1]-|
| 0000 0000 | ABCD 0000 |

bsf()は、32+16=48 を返す。


アセンブラを使わない次の方法もある。


int Bigint::bsf() const{
	static int table[32] = {0,1,10,2,11,14,22,3,30,12,15,17,19,23,26,4,31,9,13,21,29,16,18,25,8,20,28,24,7,27,6,5};

	int count = 0;
	for (int i=0; i<m_sz; i++){
		int x = m_num[i];
		if (x == 0){
			count += BITPERBLOCK;
		}else{
			count += table[((unsigned int)(x & -x) * 0x07C4ACDDUL) >> 27];
			break;
		}
	}
	return count;
}

// 64bit/Blockの場合
int Bigint::bsf64() const{
	static int table[64] = {
		0,1,59,2,60,40,54,3,61,32,49,41,55,19,35,4,62,52,30,33,50,12,14,42,56,16,27,20,36,23,44,5,
		63,58,39,53,31,48,18,34,51,29,11,13,15,26,22,43,57,38,47,17,28,10,25,21,37,46,9,24,45,8,7,6};

	int count = 0;
	for (int i=0; i<m_sz; i++){
		__int64 x = m_num[i];
		if (x == 0){
			count += BITPERBLOCK;
		}else{
			count += table[((unsigned __int64)(x & -x) * 0x03F566ED27179461ULL) >> 58];
			break;
		}
	}
	return count;
}

http://en.wikipedia.org/wiki/Count_trailing_zeros にもありますが、
M系列を使ったマジックです。
こちらの方が詳しいです。 http://d.hatena.ne.jp/siokoshou/20090704


すると、Steinのコードは次のように書ける。


//
// *this = GCD(*this, v)
//
void Bigint::gcd(const Bigint& v) {
	if (v == 0){
		return;
	}
	if (operator==(0)){
		operator= (v);
		return;
	}

	Bigint x[2] = {*this, v};

	if (x[0] < 0)	x[0].operator-();
	if (x[1] < 0)	x[1].operator-();

	int i0 = x[0].bsf();
	int i1 = x[1].bsf();
	int shift = (i0 < i1) ? i0 : i1;
	x[0] >>= i0;
	x[1] >>= shift;

	int id = 0;
	do {
		x[id^1] >>= x[id^1].bsf();

		if (x[id] > x[id^1]){
			id ^= 1;
		}
		x[id^1] -= x[id];
	} while (x[id^1] != 0);

	operator= (x[id] << shift);
}

これで30〜40%は、処理時間を短縮できました。

これこそがSteins;Gateの選択♪

このコードを実装したサンプルの16進・2進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。


(注1) 固定長ではないBigintクラスでは、伸長・縮退等のメモリの管理も必要になるので、これほど単純なものではありません。
(注2) 32bit(Block)全てを使うことが必ずしも有利ではありません。SSEを有効に使うなら桁あふれが発生しないよう 31bit 以下にします。


Tags: プログラムメモ
author : HUNDREDSOFT | - | -

カラツバ法で似非ゴースト暗算

ゴースト暗算より早い!...場合もある計算方法....(汗


ゴースト暗算というのを、少し前に知りました。
ゴースト暗算が筆算であるなら、乗算4回・加算3〜4回の演算になります。

例えば、48x53を考えてみると

【通常の筆算】
     4 8
   x 5 3
---------
     2 4
   1 2
   4 0
 2 0
---------
 2 5 4 4


【ゴースト暗算】
千の位と一の位を先に出して、残りを後から加算します。
通常の筆算の加算の順序を入れ替えていることがわかります。
(詳しくは書籍をご覧ください。魚の絵を使ったわかりやすいものです)


 (4x5)x10 + 8x5 = 240          ...千の位(2)
                   40x10 + 8x3
                  =  424         ...一の位(4)
                   42 + 4x3    ...正確には一の位と十の位で分解して加算します
                  =  54
-----------------------------
                    2544


【カラツバ】
それならばと、ゴーストっぽいカラツバ法を考えてみました。
とは言っても単純に 基数を10 にしただけですが。
同様に千の位と一の位を先に出して、残りを後から計算します。

 (4x5)x100 + (8x3) = 2024  上位桁同士の乗算と下位桁同士の乗算なので繰り上がりはない。並べるだけ。
 (4x5)     + (8x3) =  44   ,両絨棉分と下位部分の加算。20+24+2
------------------------------
                        46
 -(4-8)x(5-3)      =   8   カラツバ法のミソ。
------------------------------
                        54
------------------------------
                       2544

上の説明ではわかりにくいので、計算方法を簡単なパターンにします。

‐絨矛紊伐式矛綟瓜里魍櫃韻栃造戮泙后(→2024)
下3桁を足します(0+2+4→6)
上3桁を足します(2+0+2→4)
...ここで元の十の位と一の位の数「0,2」は忘れて、2464を記憶する。
こ櫃韻訖堯掛けられる数、それぞれ十の位と一の位の差を取って掛けます。(-4x2→-8)
ド箙翦薪(-8→+8)
と△任任た数に、イ凌瑤鯊します。(46+8→54)
併せて繋げると、2544の出来上がり。

いカラツバ法のミソで、たすき掛けの掛け算(4x3,8x5)は必要ありません。

いくつか計算サンプルを並べます。

これもカラツバの効果が良く出ています。

最も効果が出るのは、十と一の位が同じ数字の場合で、

簡単すぎて驚きます。

残念ながら全ての数字の並びで簡単ではありません。
次の例をご覧ください。

この例だと、桁上げと桁下げが発生するので、少し面倒です。

カラツバ法について詳しくは
http://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%A9%E3%83%84%E3%83%90%E6%B3%95
一般には、分割統治法として有名で、乗算の他にも除算や平方根の多倍長数計算で使われます。(注)
Karatsuba博士 もこんな使い方されるとは、思わなかったでしょう。


【まとめ】
幾つかの例では、カラツバ法の方が早いような気もしますが、
ゴースト暗算の良さは、単純化された小さな加算を組み合わせることで、小学校低学年でも計算できるところにあります。
パターンによって難易度があまり上下するのは、低学年では厳しいかもしれません。

カラツバ法の場合、乗算回数が減りますがキΔ減算になる場合もあるのでここが暗算上の鬼門になります。
例えば、乗数や被乗数に22や77など、同じ数字が並んでいる場合はきイいらなくなるので簡単ですが、最後の例などカラツバ暗算で解くのは減算による桁下げがあるので面倒です。
また、小さな桁数であれば人間にとっては、加算より乗算の方が早いので、
本来のカラツバ法の効果は、全くと言っていい程ありません。(注)

2桁2桁までしか通用しないとか「そろばん」の方が良いとか批判的な意見も散見しますが、
ゴースト暗算は、イメージ(魚パターン)で簡単にわかりやすく計算方法を解説したという点では、
今までなかった新しい試みであり、ちょっと真似してみたくなりました。


どうでも良いことですが、どこのページも「48 x 53」を参考にしているので、
48x53の答えだけは覚えちゃっている人がいるのではないだろうか。



Tags: プログラムメモ
続きを読む>>
author : HUNDREDSOFT | - | -

Windows7で簡単にマルチタッチを試す

Windows7で簡単にマルチタッチを試す方法として、 "Multi-Touch Vista" があります。
名前の通り、それほど新しいソフトではないのでご存じの方もいるでしょう。

2つのマウス入力を2つのタッチポイントに見立てて、
マルチタッチ環境を実現してくれるドライバとコントロールソフトです。

Virtual環境でこれを何とか動かせないか試したのですが、
VMware, VirtualPC, VirtualBox いずれも動かせませんでした。

ドライバの認識まではできて、コントロールパネルに
「ペンとタッチ」が追加されるところまでは行くのですが
Multitouch.Configuration.WPF.exeの"Configuer device"で
Block native windows...のチェックを入れると
マウスが一切効かなくなります。
やはり実環境が必要なのかもしれません。

実環境であれば、概ね問題なく使えるのですが、
スタートメニューやタスクマネジャーを開けるとちょっと困った問題が発生します。
.好拭璽肇瓮縫紂湿紊法▲織奪舛寮屬ぁのマークが表示されない。
▲織好マネージャー上で操作できない。

△砲弔い討蓮
Multitouch.Service.Console.exe,Multitouch.Driver.Console.exe
管理者として実行すれば可能ですが,量簑蠅あるので、
管理者権限で実行しているWindow上にタッチマークを表示することができずに、
管理者アプリのWindow背面にマークが表示されてしまいます。

"Multi-Touch Vista"は、コード公開されているので
multitouchvista-30935.zipをダウンロードして
Multitouch.InputProviders.slnをVS2010で開くと、
DebugCursor.csがマークを表示していることがわかります。

Z-Orderに関する問題ですが、常に最前面にユーザーWindowを表示する方法は、
WindowsOSのバージョンによっても変わります。

・Formの指定やCreateWindowでTopMostを指定する
・SetForegroundWindowsやAttachThreadInputの呼び出し
・BringWindowToTopの呼び出し

等がありますが、

Windows7でC#なら、Control.BringToFront()を使うのが簡単です。

DebugCursorを呼び出しているのは、DeviceStatusクラスですので
SyncUpdateLocation()に一行追加します。


void SyncUpdateLocation(object state)
{
   debugCursor.BringToFront(); //<-------------- この行を追加
   debugCursor.Location = new Point(Location.X - debugCursor.Width / 2, Location.Y - debugCursor.Height / 2);
}



これをBuildして、Outputフォルダに出力されたファイルを上書きすれば、
.好拭璽肇瓮縫紂湿紊法▲織奪舛寮屬ぁのマークが表示されない。
▲織好マネージャー上で操作できない。
の問題が解決します。




Tags: プログラムメモ VirtualPC関連 Windows7(x64)関連
author : HUNDREDSOFT | - | -