(64ビット*64ビット)%64ビットをオーバーフローせずに計算

[みけCAT]

unsigned long long型の数値a,b,cについて、(a*b)%cをオーバーフローしないように計算します。
このコードは、unsigned int型が32ビット、unsigned long long型が64ビットの環境を想定して作られています。
C言語で普通に書く方法を思いつけなかったので、インラインアセンブリで書きました。

利用(商用を含む)・改造などは自由です。転載するときは、出典(ここ)を明記してください。
プログラムのソースコードの一部としての再配布は、出典の記載があると嬉しいですが、必須ではありません。
このコードを利用したことによるいかなる損害についても、作者は責任を負いません。

コード:

/**
 * (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<<=1; */
			"shll $1,%3\n\t"
			"rcll $1,%2\n\t"
			/* if(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)<<32)|r2;
}

このコードは、日記にも投稿しました。
http://dixq.net/forum/blog.php?u=536&b=4024