_primes_bitpacked.nx source
↩ module page · 86 lines · 2372 B
1// _primes_bitpacked.nx -- bit-packed Sieve of Eratosthenes.
2//
3// 1 bit per odd number (vs 1 byte in byte-sieve). Memory: SIEVE_SIZE/16
4// bytes vs SIEVE_SIZE/2 bytes. Cache-friendlier; better passes/sec.
5//
6// Bit i represents odd number (2*i + 3), as in byte sieve.
7// Byte b at bit-offset (i % 8): bit b of byte (i / 8).
8// Bit set = COMPOSITE.
9
10import "syscalls.nx"
11
12const SIEVE_SIZE: i64 = 1000000
13const N_PASSES: i64 = 100
14
15func bit_get(buf: *u8, i: i64) -> i64 {
16 let byte_idx: i64 = i >> 3
17 let bit_idx: i64 = i & 7
18 let b: i64 = buf[byte_idx] as i64
19 return (b >> bit_idx) & 1
20}
21
22func bit_set(buf: *u8, i: i64) -> i64 {
23 let byte_idx: i64 = i >> 3
24 let bit_idx: i64 = i & 7
25 let b: i64 = buf[byte_idx] as i64
26 let mask: i64 = 1 << bit_idx
27 buf[byte_idx] = (b | mask) as u8
28 return 0
29}
30
31func one_pass(buf: *u8, half: i64, byte_count: i64) -> i64 {
32 // Zero the buffer
33 var i: i64 = 0
34 while i < byte_count {
35 buf[i] = 0 as u8
36 i = i + 1
37 }
38 // Sieve
39 var factor: i64 = 3
40 while factor * factor <= SIEVE_SIZE {
41 let factor_idx: i64 = (factor - 3) / 2
42 var fi: i64 = factor_idx
43 var done: i64 = 0
44 while done == 0 {
45 if fi >= half { done = 1 }
46 if done == 0 {
47 if bit_get(buf, fi) == 0 { done = 1 }
48 if done == 0 { fi = fi + 1 }
49 }
50 }
51 if fi >= half { return -1 }
52 let cur_num: i64 = 2 * fi + 3
53 var mult: i64 = cur_num * cur_num
54 while mult <= SIEVE_SIZE {
55 let idx: i64 = (mult - 3) / 2
56 bit_set(buf, idx)
57 mult = mult + 2 * cur_num
58 }
59 factor = cur_num + 2
60 }
61 // Count: 2 is prime, plus all unset bits below half representing odd <= SIEVE_SIZE
62 var count: i64 = 1
63 var k: i64 = 0
64 while k < half {
65 if bit_get(buf, k) == 0 {
66 let n: i64 = 2 * k + 3
67 if n <= SIEVE_SIZE { count = count + 1 }
68 }
69 k = k + 1
70 }
71 return count
72}
73
74func main() -> i64 {
75 let half: i64 = SIEVE_SIZE / 2
76 let byte_count: i64 = (half + 7) / 8 // round up
77 let buf: *u8 = sys_mmap(byte_count)
78 var p: i64 = 0
79 var c: i64 = 0
80 while p < N_PASSES {
81 c = one_pass(buf, half, byte_count)
82 p = p + 1
83 }
84 if c != 78498 { return 2 }
85 return 0
86}