view example/Miller_Rabin/spe/Prime.cc @ 1244:cd50c48f45e7 draft real_matrix

fix
author Daichi Toma <amothic@gmail.com>
date Thu, 03 Nov 2011 20:40:17 +0900
parents 9d37fa6bc1da
children 163207b736c5
line wrap: on
line source

#include <stdio.h>
#include "SchedTask.h"
#include "Prime.h"
#include "Func.h"

typedef unsigned long long U64;

static U64 base[] = { 2, 3, 5, 7, 11 };
#define MRMAX 5

SchedDefineTask1(Prime, prime);

static U64
powMod(U64 base, U64 power, U64 mod)
{
	U64 result = 1;

	while (power > 0) {
		if ( power & 1 == 1 ) {
			result = (result * base) % mod;
		}
		base = ( base * base ) % mod;
		power >>= 1;
	}
	return result;
}

static bool
millarRabin(U64 a, U64 s, U64 d, U64 n) {
	U64 x = powMod(a, d, n);
	if ( x == 1) return true;
	for ( U64 r = 0; r < s; r++) {
		if ( x == n-1 ) return true;
		x = (x * x) % n;
	}
	return false;
}

static bool
isPrime(U64 n) {

	if ( n <= 1 ) return false;

	U64 d = n - 1;
	U64 s = 0;

	while( (d & 1) == 0 ) {
		d >>= 1;
		s++;
	}

	for ( int i = 0; i < MRMAX; i++ ) {
		if ( n == base[i] ) return true;
		if (!millarRabin(base[i], s, d, n)) return false;
	}
	return true;
}

	static int
prime(SchedTask *smanager, void *rbuf, void *wbuf)
{
	U64 start = (U64)smanager->get_param(0);	/* 素数判定の開始地点 */
	U64 end   = (U64)smanager->get_param(1);	/* 素数判定の終了地点 */
	U64 range = end - start;					/* 判定する範囲		  */
	U64 index_range = range >> 1;

	/* 判定結果を収める配列を受け取る */
	bool *output = (bool*)smanager->get_output(wbuf, 0);

	/* 初期化 */
	for (U64 i = 0; i < index_range; i++){
		output[i] = false;
	}

	for (U64 i = start + 1,index = 0; i < end ; i += 2, index++) {
		if (isPrime(i)) {
			output[index] = true;
		}
	}
	return 0;
}