Mercurial > hg > Game > Cerium

comparison example/Miller_Rabin/ppe/Prime.cc @ 1242:9d37fa6bc1da draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
add Miller_Rabin
author Daichi Toma <amothic@gmail.com>
date Tue, 01 Nov 2011 19:01:13 +0900
parents
children cd50c48f45e7
comparison
equal deleted inserted replaced
1241:b96b1604ab85 1242:9d37fa6bc1da
1 #include <stdio.h>
2 #include "SchedTask.h"
3 #include "Prime.h"
4 #include "Func.h"
5
6 typedef unsigned long long U64;
7
8 static U64 base[] = { 2, 3, 5, 7, 11 };
9 #define MRMAX 5
10
11 SchedDefineTask1(Prime, prime);
12
13 static U64
14 powMod(U64 base, U64 power, U64 mod)
15 {
16 U64 result = 1;
17
18 while (power > 0) {
19 if ( power & 1 == 1 ) {
20 result = (result * base) % mod;
21 }
22 base = ( base * base ) % mod;
23 power >>= 1;
24 }
25 return result;
26 }
27
28 static bool
29 millarRabin(U64 a, U64 s, U64 d, U64 n) {
30 U64 x = powMod(a, d, n);
31 if ( x == 1) return true;
32 for ( U64 r = 0; r < s; r++) {
33 if ( x == n-1 ) return true;
34 x = (x * x) % n;
35 }
36 return false;
37 }
38
39 static bool
40 isPrime(U64 n) {
41
42 if ( n <= 1 ) return false;
43
44 U64 d = n - 1;
45 U64 s = 0;
46
47 while( (d & 1) == 0 ) {
48 d >>= 1;
49 s++;
50 }
51
52 for ( int i = 0; i < MRMAX; i++ ) {
53 if ( n == base[i] ) return true;
54 if (!millarRabin(base[i], s, d, n)) return false;
55 }
56 return true;
57 }
58
59 static int
60 prime(SchedTask *smanager, void *rbuf, void *wbuf)
61 {
62 U64 start = (U64)smanager->get_param(0); /* 素数判定の開始地点 */
63 U64 end = (U64)smanager->get_param(1); /* 素数判定の終了地点 */
64 U64 range = end - start; /* 判定する範囲 */
65
66 /* 判定結果を収める配列を受け取る */
67 bool *output = (bool*)smanager->get_output(wbuf, 0);
68
69 /* 初期化 */
70 for (U64 i = 0; i < range; i++){
71 output[i] = true;
72 }
73
74 for (U64 i = start + 1,index = 0; index < range; i += 2, index++) {
75 if (!isPrime(i)) {
76 output[index] = false;
77 }
78 }
79 return 0;
80 }