みけCATのにっき（仮）
つれづれなるまゝに、日くらし、PCにむかひて、心に移りゆくよしなし事を、そこはかとなく書きつくれば、あやしうこそものぐるほしけれ。
（本当か！？）
出典

(64ビット*64ビット)%64ビットの計算

アバター
みけCAT
記事: 6734
登録日時: 15年前
住所: 千葉県
連絡を取る:
みけCAT に連絡する

(64ビット*64ビット)%64ビットの計算

投稿記事 by みけCAT » 12年前

C言語でunsigned long long型の変数a,b,cがあり、aとbの積をcで割った余りを求めたいとき、どうすればいいでしょうか？
普通に(a*b)%cと書くのでは、オーバーフローの危険があります。
また、もしa,b,cがunsigned int型だったらunsigned long long型にキャストして計算すればいいですが、
unsigned long long型の上位の整数型はなさそうです。
多倍長の実装もめんどくさいし、これだけのために多倍長ライブラリを使うのも嫌だ。
では、どうすればいいでしょうか？

まず、「繰り返し二乗法」を参考にした「繰り返し二倍法」を思いつきました。

CODE: 

unsigned long long repeat_add(
		unsigned long long a,unsigned long long b,unsigned long long c) {
	unsigned long long result=0;
	while(b>0) {
		if(b&1)result=(result+a)%c;
		a=(a>=1;
	}
	return result;
}
このコードでa,b,cがそれぞれLLONG_MAX以下(63ビット)なら計算できますが、64ビットだとオーバーフローしそうです。
それでも、素朴な方法の32ビットよりはかなり伸びました。

とりあえず、より高速に動作するようにしてみましょう。

CODE: 

unsigned long long repeat_add(
		unsigned long long a,unsigned long long b,unsigned long long c) {
	unsigned long long result=0;
	if(a>=c)a%=c;
	while(b>0) {
		if(b&1) {
			result+=a;
			if(result>=c)result-=c;
		}
		a+=a;
		if(a>=c)a-=c;
		b>>=1;
	}
	return result;
}
メインの処理では割り算が消え、より高速に動作するようになりました。
そして、割り算が消えたということは、インラインアセンブリで書けそうです。
実際に書いてみました。
unsigned long long型が64ビット、unsigned int型が32ビットであると仮定しています。

CODE: 

#include 

/**
 * (a*b)%cを計算する
 * @param a 掛けられる数
 * @param b 掛ける数
 * @param c あまりを取る時に割る数
 * @return (a*b)%cの計算結果
 */
unsigned long long repeat_add(
		unsigned long long a,unsigned long long b,unsigned long long c) {
	unsigned int a1,a2,b1,b2,c1,c2,r1,r2;
	if(a>=c)a%=c;
	a1=a>>32;a2=a;b1=b>>32;b2=b;c1=c>>32;c2=c;
	__asm__ volatile (
		"movl %6,%%ecx\n\t"
		"movl %7,%%edx\n\t"
		/* result=0; */
		"movl $0,%0\n\t"
		"movl $0,%1\n\t"
		/* while(b>0) { */
		"ra_while:\n\t"
		"movl %4,%%eax\n\t"
		"orl %5,%%eax\n\t"
		"jz ra_exit\n\t"
			/* if(b&1) { */
			"movl %5,%%eax\n\t"
			"andl $1,%%eax\n\t"
			"jz ra_noif1\n\t"
				/* result+=a; */
				"movl %2,%%eax\n\t"
				"movl %3,%%ebx\n\t"
				"addl %%ebx,%1\n\t"
				"adcl %%eax,%0\n\t"
				/* if(result>=c) */
				"jc ra_do_rmc\n\t"
				"cmpl %%ecx,%0\n\t"
				"ja ra_do_rmc\n\t"
				"jb ra_noif1\n\t"
				"cmpl %%edx,%1\n\t"
				"jae ra_do_rmc\n\t"
				"jmp ra_noif1\n\t"
				/* result-=c; */
				"ra_do_rmc:\n\t"
				"sub %%edx,%1\n\t"
				"sbb %%ecx,%0\n\t"
			/* } */
			"ra_noif1:\n\t"
			/* a=c) */
			"jc ra_do_amc\n\t"
			"cmpl %%ecx,%2\n\t"
			"ja ra_do_amc\n\t"
			"jb ra_dont_amc\n\t"
			"cmpl %%edx,%3\n\t"
			"jae ra_do_amc\n\t"
			"jmp ra_dont_amc\n\t"
			/* a-=c; */
			"ra_do_amc:\n\t"
			"subl %%edx,%3\n\t"
			"sbbl %%ecx,%2\n\t"
			/* b>>=1; */
			"ra_dont_amc:\n\t"
			"shrl $1,%4\n\t"
			"rcrl $1,%5\n\t"
		/* } */
		"jmp ra_while\n\t"
		"ra_exit:\n\t"
	: "=m"(r1),"=m"(r2)
	: "m"(a1),"m"(a2),"m"(b1),"m"(b2),"m"(c1),"m"(c2)
	: "%eax","%ebx","%ecx","%edx"
	);
	return (((long long)r1)) {
	my $now=$_;
	chomp($now);
	my ($as,$bs,$cs)=split(/ +/,$now);
	my $a=Math::BigInt->new($as);
	my $b=Math::BigInt->new($bs);
	my $c=Math::BigInt->new($cs);
	my $ans=($a*$b)%$c;
	print "$ans\n";
}
ランダムケースジェネレータでテストケースを生成し、C言語版とPerl版の出力を比較。

CODE: 

#include 
#include 
#include 

static unsigned int x=123456789;
static unsigned int y=362436069;
static unsigned int z=521288629;
static unsigned int w=88675123;

void seed(unsigned int a,unsigned int b,unsigned int c,unsigned int d) {
	if((a|b|c|d)==0)return;
	x=a;y=b;z=c;w=d;
}

unsigned int randint(void) {
	unsigned int t;
	t=(x^(x>19))^(t^(t>>8));
	return w;
}

unsigned long long randll(void) {
	return (((long long)randint())<<32)|randint();
}

int main(void) {
	int i;
	srand((unsigned int)time(NULL));
	seed(rand(),rand(),rand(),rand());
	for(i=0;i<1000;i++) {
		printf("%llu %llu %llu\n",randll(),randll(),randll());
	}
	for(i=0;i<1000;i++) {
		printf("%llu %llu %u\n",randll(),randll(),randint());
	}
	for(i=0;i<1000;i++) {
		printf("%u %u %llu\n",randint(),randint(),randll());
	}
	for(i=0;i<1000;i++) {
		printf("%u %u %u\n",randint(),randint(),randint());
	}
	return 0;
}
無事、C言語版の出力とPerl版の出力が一致したようです。

今回作成したコードは、サンプルを共有するコミュニティの方にも投稿しました。
http://dixq.net/forum/viewtopic.php?f=83&t=13816
ページトップ
コメントはまだありません。