「素数」を編集中 (節単位)

===高速な方法===
<source lang="c">
/* suppose that a<c */
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+=a;
			if(result>=c)result-=c;
		}
		a<<=1;
		if(a>=c)a-=c;
		b>>=1;
	}
	return result;
}

/* suppose that n<(1ll<<63) */
int isPrime(unsigned long long n) {
	const int smallPrimes[30]={
		 2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
		31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
		73, 79, 83, 89, 97,101,103,107,109,113
	};
	int i,ok;
	unsigned long long d,a;
	unsigned long long pownow,powcurrent,powresult;
	int r,s;

	if(n<2)return 0;
	for(i=0;i<30;i++) {
		if(n==smallPrimes[i])return 1;
		if(n%smallPrimes[i]==0)return 0;
	}

	s=0;d=n-1;
	while((d&1)==0){s++;d>>=1;}
	for(i=0;i<12;i++) {
		powresult=1;
		pownow=smallPrimes[i];
		powcurrent=d;
		while(powcurrent>0) {
			if(powcurrent&1)powresult=repeat_add(powresult,pownow,n);
			pownow=repeat_add(pownow,pownow,n);
			powcurrent>>=1;
		}
		if(powresult==1)continue;
		ok=1;
		for(r=0;r<s;r++) {
			if(powresult==n-1){ok=0;break;}
			if(powresult==1)break;
			powresult=repeat_add(powresult,powresult,n);
		}
		if(ok)return 0;
	}
	return 1;
}
</source>
処理は複雑だが、ある程度大きな数でもかなり速く判定してくれる。