ニュートン法を使った2進多桁の整数平方根

多桁の整数値(s)の平方根 √s があったとして
ニュートン法(2次漸化式)を使って(s)の平方根の整数部分を求める方法について考えます。
与える引数(s)に基数のべき乗を掛けることで任意の精度に拡張できます。 __前回同様__ 高次漸化式にも応用できますが、ここでは2次漸化式を使っています。

√s を展開すると除算が必要になるので, (1/√s)を求め、これに(s)を乗じることとします。
(1/√s) の2次漸化式は、

X(i+1) = X(i) * (3 - s * X(i)^2) / 2

と表せます。
これは、

Bigint iSqrt(Bigint s){
  int n = s.length(); // bit数

  Bigint x = 1;
  x <<= (n / 2);  // Newton法の初期値として, 1 << (n/2) を与える。
  Bigint m(0);
  Bigint t3 = 3;
  t3 <<= n;

  while(m != x){
    m = x;
    x *= t3 - ((s * x * x) >> n);
    x >>= (n + 1);
  }
  x *= s;         // x = s / √s = √s
  x >>= n;

  // 誤差修正
  m = x + 1;
  if (s >= m * m){
     x = m;
  }
  return x;
}


と書けます。
考え方としては、計算各項を (* 2n) することで整数を一時的に固定小数点化するものです。
演算量を減らす変形を重ねているのでわかりにくいかもしれませんが、非常に短いコードです。
初期値については、平方根は (n/2)桁程度になるので、その逆数もまた 2n倍することを考えれば、
(n/2)桁程度であることから決定していますが、初期値が甘いのでループ回数はやや多いです。
ニュートン法を使った除算については、__こちら__をご覧ください。

「Newton法初期値の改良バージョン」はこちら
ニュートン法を使った2進多桁の整数平方根(その2)

(注)ここでのBigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。


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

ニュートン法を使った2進多桁の整数除算

除算 x = p / q があったとして
ニュートン法(2次漸化式)を使って(q)の逆数を求め、それに(p)を乗算する方法は多桁演算の定石です。
この場合、多桁の浮動小数点型に変換して演算することが多いでしょう。
ここでは、ニュートン法を整数演算のみで実装する方法について考えます。
高次漸化式にも応用できますが、ここでは2次漸化式を使っています。

1/q の2次漸化式は、

X(i+1) = X(i) * (2 - q * X(i))

と表せます。
ここで被除数のBit数をpb, 除数のBit数をqbとして、n = pb + qb とすれば


Bigint iDiv(Bigint p, Bigint q){
  int pb = p.length(); // 被除数のbit数
  int qb = q.length(); // 除数のbit数
  int n = pb + qb;

  Bigint m(0);
  Bigint x = 1;
  x <<= pb;      // Newton法の初期値として, 1 << pb を与える。
  Bigint c2 = 2;
  c2 <<= n;

  while(m != x){
    m = x;
    x *= c2 - q * x;
    x >>= n;
  }
  x *= p;        // x = p / q
  x >>= n;

  // 誤差修正
  if (p >= (x+1)*q){
      x += 1;
  }
  return x;
}


と書けます。
考え方としては、計算各項を (* 2n) することで整数を一時的に固定小数点化するものです。
演算量を減らす変形を重ねているのでわかりにくいかもしれませんが、非常に短いコードです。
但し、初期値が甘いのでループ回数はやや多いです。
被除数と除数の桁差が少ない場合は、回復法を使った2進多桁の除算が有利です。


「Newton法初期値の改良バージョン」はこちら
ニュートン法を使った2進多桁の整数除算(その2)

(注)ここでのBigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。


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

モンゴメリ乗算 サンプルコード

ウィキぺディアでモンゴメリ乗算を検索すると、わかり易く詳しい説明が載っている。
残念ながらこの説明全体をプログラミングしたものが見つからなかったのでC++ソースを書いてみた。


// モンゴメリ乗算剰余演算 c = (a*b) mod nを求める
// gcd(n,r) == 1 となること。(nは奇数であれば良い)
//
Bigint Montgomery(const Bigint& a, const Bigint& b, const Bigint& n) const{
  // n を法とする
  // Bit長を調べる(nb)
  int nb = n.length();

  // r: nより大きな2のべき
  Bigint r(1);
  r <<= nb;

  // r2 = r^2 mod n
  Bigint r2(r*r);
  r2 %= n;

  // n*n2 ≡ -1 mod r
  Bigint n2 = 0; /* 求めるN' */
  Bigint t = 0;
  Bigint vi = 1;
  int ni = nb;
  while (ni--){ /* Rのトップビットを除いたビット数分繰り返す */
    if ((t & 1) == 0){
      /* ゼロになっているビットがあったら、N'のその部分を1にする(NはRと互いに素なので必ず奇数)*/
      t += n; /* 掛け算だが、二進数一桁の掛け算なので実質は足し算 */
      n2 += vi; /* N'のその部分を1にする */
    }
    t >>= 1; /* 必ず端数が出るが切り捨てる */
    vi <<= 1; /* Rは2の冪なので、絶対端数は出ない */
  }
  // ここまでで、今後計算に必要になる、r2,n,nb,n2が得られる。
  // つまりnを法とするモンゴメリクラスを作成するなら
  // 引数nをコンストラクタとするクラスを作成し、
  // r2,n,nb,n2をメンバ変数とする。


  // aのモンゴメリ表現をam, bのモンゴメリ表現をbmとする
  // t は作業領域
  t = a * r2;
  Bigint am = t * n2;
  am.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  am *= n;
  am += t;
  am >>= nb; // 1/Rの意味
  if (am >= n) am -= n;

  t = b * r2;
  Bigint bm = t * n2;
  bm.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  bm *= n;
  bm += t;
  bm >>= nb; // 1/Rの意味
  if (bm >= n) bm -= n;

  // cm: am * bm のモンゴメリリダクション
  t = (am * bm);
  Bigint cm = t * n2;
  cm.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  cm *= n;
  cm += t;
  cm >>= nb; // 1/Rの意味
  if (cm >= n) cm -= n;

  // cmのモンゴメリリダクション
  t = cm;
  Bigint c = t * n2;
  c.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  c *= n;
  c += t;
  c >>= nb; // 1/Rの意味
  if (c >= n) c -= n;

  return c;
}



ここでは、できるだけwikipediaの説明と合致するように書きましたが、
実利用の際には、r2 = r^2 mod n の計算が重いので、
(これが早ければ、そもそもモンゴメリ演算の必要がない)
nを法とするモンゴメリクラスを作成するのが良いです。
つまりコード中にも書いてあるように、引数nをコンストラクタとしたモンゴメリクラスを作成し、r2,n,nb,n2はメンバ変数にして、モンゴメリリダクション, 乗算剰余演算, 加減算剰余演算, 冪剰余演算等は、モンゴメリクラスのMethodにします。

※ここで書いた乗算は「乗算剰余演算」の項の上段にある方法です。

(注)ここでのBigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。


Tags: モンゴメリ冪剰余サンプルコード     プログラムメモ
author : HUNDREDSOFT | - | -

回復法を使った2進多桁の除算(その2)

前回同様、除算 a = x / y があったとして
式を変形して、

a = (x * (2^(2k)-1)/y + 2^k) >> 2k  (k:=被除数のbit数)

とすれば、逆数を求める部分((2^(2k)-1)/y)の回復法もまた特殊な被除数のため、簡単なものに変形できる。

Bigint iDiv(Bigint x, Bigint y){
  int xb = x.length(); // 被除数のbit数
  int yb = y.length(); // 除数のbit数
  int n = xb * 2 - 1; // ループ回数

  Bigint a(0), m(0);
  while (n--){
    m |= 1;
    a <<= 1;
    if (m >= y){
      m -= y;
      a |= 1;
    }
    m <<= 1;
  }
  a *= x;
  m = 1;
  m <<= xb;
  a += m;
  return a >> (xb*2-1);
}

この方法は、被除数bit数の倍のループとなり、これが最大の減算回数になります。
しかし誤差が出ないので乗算は1回だけで済みます。
前回の方法はbit差に依存したループ回数なので多桁にも応用できますが、この方法だと被除数が大きいとループ回数が多くなってしまいます。

桁差が小さい時に有利な方法

ニュートン法を使った2進多桁の整数除算

(注)ここでのBigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。


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

回復法を使った2進多桁の除算

除算 a = x / y があったとして
回復法を使うと被除数のbit数分の加算・減算を繰り返し行うことになる
式を変形して、

a = (x * ((2^k)/y) >> k  (k:=被除数のbit数)

とすれば、逆数を求める部分((2^k)/y)の回復法は特殊な被除数のため、簡単なものに変形できる。

Bigint iDiv(Bigint x, Bigint y){
  Bigint a(0), m(1);
  int xb = x.length(); // 被除数のbit数
  int yb = y.length(); // 除数のbit数
  int n = xb - yb + 1; // ループ回数(bit差+1)

  m <<= yb-1;
  while (n--){
    a <<= 1;
    if (m >= y){
      m -= y;
      a |= 1;
    }
    m <<= 1;
  }
  a *= x;
  a >>= (xb - 1);

  // 誤差修正(ここでのaは,正解より1or2小さい場合がある)
  m = y * a;
  while(x > m){
    m += y;
    if (x > m)
      a += 1;
  }
  return a;
}


と書ける。
計算量は最大でもbit差程度の加減算と2回の乗算なので
bit差が数100のオーダーならNewton法を使うより有利かも。

ニュートン法を使った2進多桁の整数除算

(注)ここでのBigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。


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

PerformanceCounterを使ったマイクロ秒時計

PerformanceCounterを使ったマイクロ秒時計

Linuxには、gettimeofdayがありますが、
Windowsには、マイクロ秒まで計測可能なAPIは用意されていません。

PerformanceCounterを使えばかなり正確な時間間隔を知ることはできますが、
時計として利用することは簡単ではありません。
ある程度観測をしていると、ansi標準の_ftimeとはズレていきます。
このため、_ftimeをベースとしてミリ秒以下の部分をPerformanceCounterから得ることを考えましたが、秒の変わり目で時計が逆に進む場合があります。
そこで、ftimeとPerformanceCounterから得られる時間差を調節し、
時計が逆に進むことがないようなマイクロ秒時計を作ってみました。

usec.h

// (c) Hundredsoft Corporation 2010 All rights reserved.

#ifndef USEC_INCLUDED
#define USEC_INCLUDED
class usec
{
public:
	usec();
	~usec(){}
	__int64 get();
private:
	__int64 microfunc(__int64 a, __int64 b);
	__int64 usec::getstm();
	LARGE_INTEGER m_ti;
	LARGE_INTEGER m_tfi;
	__int64 m_msi;
	__int64 m_lastus;
	bool m_isSuported;
};
#endif



usec.cpp

// (c) Hundredsoft Corporation 2010 All rights reserved.

#include <windows.h>
#include <time.h>
#include <sys/timeb.h>
#include "usec.h"

usec::usec()
{
	m_isSuported = false;
	if (QueryPerformanceFrequency(&m_tfi)){
		if (QueryPerformanceCounter(&m_ti)){
			m_isSuported = true;
		}
	}
	m_msi = getstm();
}

__int64 usec::microfunc(__int64 a, __int64 b)
{
	__int64 ah, al, a0, m = 0;
	int i;
	ah = (a >> 44) & 0xffff;
	al = a << 20;
	for (i=0; i<128; i++){
		m <<= 1;
		if (ah < 0) m |= 1;
		ah <<= 1;
		if (al < 0) ah |= 1;
		al <<= 1;
		if (m >= b){
			m -= b;
			al |= 1;
		}
	}
	a0 = al;
	ah = (a0 >> 44) & 0xffff;
	al = a0 << 20;
	m = 0;
	for (i=0; i<128; i++){
		m <<= 1;
		if (ah < 0) m |= 1;
		ah <<= 1;
		if (al < 0) ah |= 1;
		al <<= 1;
		if (m >= 22634874){
			m -= 22634874;
			al |= 1;
		}
	}
	return a0 - al;
}

__int64 usec::getstm()
{
	_timeb  tb;
	_ftime(&tb);
	struct tm* stm = localtime(&tb.time);
	__int64 stms = (stm->tm_hour * 3600000 + stm->tm_min * 60000 + stm->tm_sec * 1000 + tb.millitm);
	return stms * 1000;
}

__int64 usec::get()
{
	__int64 utm;
	if (m_isSuported){
		LARGE_INTEGER t, tf;
		QueryPerformanceFrequency(&tf);
		QueryPerformanceCounter(&t);
		__int64 stms = getstm();
		__int64 diff = t.QuadPart-m_ti.QuadPart;
		if (diff > 0){
			utm = m_msi + microfunc(diff, tf.QuadPart);
			diff = utm - stms;
		}else{
			utm = stms;
			diff = 999999999;
		}
		__int64 absdiff = (diff < 0) ? -diff : diff;
		if (m_tfi.QuadPart != tf.QuadPart || absdiff > 10000){
			m_tfi = tf;
			m_ti = t;
			m_msi = stms;
			if (diff > 0){
				if (utm - m_lastus > 10 )
					utm = (m_lastus + utm) / 2;
				else
					utm = m_lastus + 2;
			}else{
				if (absdiff > 3600000000)
					utm = m_msi;
				else
					utm = (utm + m_msi) / 2;
			}
		}
		m_lastus = utm;
	}else{
		utm = getstm();
	}
	return utm;
}





テスト用コード
main.cpp

#include <tchar.h>
#include <windows.h>
#include <stdio.h>
#include "usec.h"

#include <time.h>
#include <sys/timeb.h>

int main(int argc, char**argv)
{
	usec ut;
	while(1){
		char work[256];
		__int64 utm = ut.get();

		_timeb  tb;
		_ftime(&tb);
		struct tm* stm = localtime(&tb.time);

		sprintf(work, "%02d:%02d:%02d.%06d [%02d.%03d]",
				(int)((utm / 3600000000) % 24),
				(int)((utm / 60000000) % 60),
				(int)((utm / 1000000) % 60),
				(int)(utm % 1000000),
				stm->tm_sec,
				tb.millitm
				);
		printf("%s\n", work);
		Sleep(300);
	}
}
usec::microfuncは、PerformanceCounterから得られる値(LARGE_INTEGER)から
マイクロ秒を得る128bit除算を行っています。

ある時間差(t0〜t1)があったとして
t0,t1はQueryPerformanceCounterによって得られる値とします。
t = t1 - t0、この時間間隔は、周波数:f(QueryPerformanceFrequencyによって得られる)として
t / f で得られますが、
整数演算でマイクロ秒を得たいので
t * 1000000 / f
となります。
このままdoubleで演算しても良いのですが整数演算に拘って、
t * (1024*1024 - 48576) / f
= 1024*1024 * t / f * (1 - 1 / 21.58629776021)
= 1024*1024 * t / f * (1 - 1024*1024 / 22634874)
A = (t << 20) / f
と置けば
≒A - (A << 20) / 22634874 [μSec]となります。
これは、64bitを超える演算になりますのでforで回して計算しています。

テスト用コードでは、usecクラスから得られたマイクロ秒時間と共に、
_ftimeによる「秒、ミリ秒」も表示しています。

usec::getでは、ftimeとPerformanceCounterから得られる時間差を調節しているので、
それほど精度の高い時計ではありませんが、ミリ秒単位の時計だとループ間隔を短くSleep(10)とかにすると、同じ時間が続いてしまいます。
この時計の場合は、ループ間隔(Sleep)を短くしても、同じ時間が並ぶようなことはありません。

※高分解能パフォーマンスカウンタがないPCでは、ミリ秒単位の出力となります。



usec.zipソースはこちら



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

自己署名でuiAccess=trueにする

自己署名を使って、uiAccess=trueを試してみた。

※Windows7 x64で確認しましたが、Vistaでは確認しておりません。


(1) VS2008でC++/CLRの適当なプロジェクトを作り、ビルド&起動確認。

wizardのままのFormだけで構いません。


(2) プロジェクトプロパティ→リンカ→マニフェストファイル

「UACによるUI保護のバイパス」を「はい」にする。


(3) ビルド&起動確認

「サーバーから紹介が返ってきました」と表示される。


(4) できたexeを "C:\Program Files\test2" にコピーする。

ここでは、"clitest4.exe" としています。

(5) 管理者でコマンドプロンプトを開き、

"C:\Program Files (x86)\Microsoft Visual Studio 9.0\Common7\Tools\vsvars32.bat"

(パスを入れるのが面倒なので)


(6) フォルダ移動

cd "C:\Program Files\test2"


(7) ルート証明書を作る

makecert -n "CN=testroot,O=aaa,C=JP,E=foo@hoge.jp" -sv testroot.pvk -r testroot.cer

パスワードを聞いてきます。

clip_1.gif















最後にもう一度尋ねられるので、先ほどのパスワードを入れます。


(8) ルート証明書を使ってデジタル署名用の証明書

makecert -n "CN=test,O=aaa,C=JP,E=foo@hoge.jp" -sv test.pvk -ic testroot.cer -iv testroot.pvk test.cer

先ほどと同じようにパスワードを聞いてきます。rootとは別のパスワードが良いです。
フォルダはこんな感じ。

clip_18.gif


























(9) testroot.cerをダブルクリックしてルート証明書をインストール

clip_6.gif

































信頼されたルート証明機関にストアします。

clip_7.gif





























(10) 今度は、test.cerをダブルクリックして証明書をインストール

clip_8.gif































自動で個人に入ると思います。

clip_9.gif































(11) mmcでスナップインの追加→証明書を追加して

 信頼されたルート証明機関にtestroot
 個人にtestが追加されていることを確認します。

clip_20.gif































clip_21.gif




























(12) Signtool

先ほどのコマンドプロンプトで

Signtool signwizard

と入れると、デジタル署名ウィザードが始まります。

exeを選びます。

clip_12.gif





























カスタムを選びます。

clip_13.gif





























次の画面で「ファイルから選択」をクリック。

clip_14.gif




























X.509証明書の"test.cer"を選びます。

clip_15.gif


























次の画面で秘密キーファイルに、test.pvkを選びます。


パスワード画面後、「ルート証明書を含む証明のパスのすべての証明書」を選んで次へ。

clip_16.gif


























最後の完了でもう一度パスワードを入れます。


(13) 確認

C:\Program Files\test2\clitest4.exe

を動かすと、普通に起動します。

MT.exeで確認しても良いですし、バイナリを直接見ても、
uiAccess="true"が設定されていることがわかります。

clip_22.gif





























いきなりexeをバイナリで見る人もいないとは思いますが....

ルート証明書がないとやっぱり「サーバーから紹介が返ってきました」と表示される。

この方法だとBuildの度にサインしなければなりませんが、
2回目以降は、証明書はできているので、Signtool signwizard だけで良いです。
入力するものが決まれば、コマンドプロンプトでSigntoolを動かせるので、
ビルド後のイベントに追加して自動化できます。

念のため書くと、uiAccess="true"にしたコードは動作はしますが、
・証明書のインストールが必要。
・Windowsソフトウェアロゴを得るには、MSの特別許可が必要。

ということなので、汎用ソフトというより開発ツールの類でしょうか。




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

requestedExecutionLevelのuiAccess

前回の記事で書いたrequestedExecutionLevelの
2番目のエレメントuiAccessについて

uiAccess="true" にすると上位権限へのメッセージ送信が可能になります。

つまり標準権限で起動したアプリでも、管理者権限で実行中のアプリに対して
ドラッグ&ドロップが可能になります。いわゆるUIPI問題の回避策です。

例えばExplorerで uiAccess="true" に設定できれば、管理者権限で起動したアプリにドラッグ&ドロップできるようになるので、こういった設定があればよいのですが、セキュリティ効果を半減させてしまうのでバランスが難しいです。

但し、uiAccess="true" には2つの条件があります。

(1) 署名付きのプログラムである。
(2) %WinDir% または、%ProgramFiles%フォルダーに保存されている。


単に MT.exe で変更しただけでは、条件を満たしていないので、
「サーバーから紹介が返ってきました」と表示されます。
(この文言が適当なのでしょうかね。)

オープンソースや自作ソフト・自社製品以外で(1)を満たすことはまず無理なので、一般のアプリケーションやMSのプログラムを uiAccess="true" に書き換えることは無理でしょう。

バイナリを直接いじるという手もありますが、ダイジェスト値が変わってしまうので、やはり認証されないと思います。

別プロセスでエクスプローラを管理者権限で立ち上げて、
そこからドラッグ&ドロップするという手はありますが、

ドラドロするようなアプリは、メニューにFileOpenがあるのが普通なので
そこまでするなら、FileOpenを使った方が簡単な気がする。

ドラッグ&ドロップだけでも何とかならないですかね。

----------------
【追記】
自己署名でuiAccess=true を試される方は、こちら




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