math.c
Raw
1#include <cutils/math.h>
2
3uintmax_t isqrt(uintmax_t n)
4{
5 uintmax_t start = 1, end = n/2, ans = 0;
6 if(n == 0 || n == 1)
7 return n;
8
9 while(start <= end)
10 {
11 uintmax_t mid = (start + end) / 2;
12
13 if(mid*mid == n)
14 return mid;
15
16 if(mid*mid < n)
17 {
18 start = mid + 1;
19 ans = mid;
20 } else {
21 end = mid - 1;
22 }
23 }
24 return ans;
25}
26
27intmax_t ipow(intmax_t base, uintmax_t exp)
28{
29 intmax_t result = 1;
30 while (exp)
31 {
32 if (exp & 1)
33 result *= base;
34 exp >>= 1;
35 base *= base;
36 }
37
38 return result;
39}
40
41bool is_prime(uintmax_t n)
42{
43 uintmax_t i;
44
45 if (n < 2) return false;
46 if (n % 2 == 0 && n != 2) return false;
47
48 for(i = 3; i < isqrt(n); i+=2)
49 {
50 if (n % i == 0)
51 return false;
52 }
53 return true;
54}
55
56static void primesieve_set(byte* numbers, size_t index)
57{
58 numbers[index/8] |= (byte)(1 << (index%8));
59}
60
61static bool primesieve_get(byte* numbers, size_t index)
62{
63 return getbit((numbers[index/8] >> index%8), 0);
64}
65
66Bytearray* primesieve(uintmax_t limit)
67{
68 size_t i, j, length = limit/8+1;
69 byte* numbers;
70 Bytearray* bt;
71
72 numbers = calloc(length, 1);
73 if(!numbers)
74 return NULL;
75
76 bt = new_bytearray(sizeof(uintmax_t));
77 if(!numbers)
78 {
79 free(numbers);
80 return NULL;
81 }
82
83
84 for(i = 2; i < isqrt(limit); i++)
85 {
86 if(primesieve_get(numbers, i) == 0)
87 {
88 for(j = i*i; j < limit; j+=i)
89 {
90 primesieve_set(numbers, j);
91 }
92 }
93 }
94
95 for(i = 0; i < limit; i++)
96 {
97 if(primesieve_get(numbers, i) == 0)
98 bytearray_push(bt, &i);
99 }
100
101 free(numbers);
102 return bt;
103}
104