C++ Mathematical Expression Toolkit (ExprTk) release
Loading...
Searching...
No Matches
Classes | Functions
exprtk_miller_rabin_primality_test.cpp File Reference
#include <algorithm>
#include <array>
#include <cstdio>
#include <random>
#include <string>
#include "exprtk.hpp"
Include dependency graph for exprtk_miller_rabin_primality_test.cpp:

Go to the source code of this file.

Classes

struct  random_uint< T >
 

Functions

template<typename T >
void miller_rabin_probabilistic_primality_test ()
 
int main ()
 

Function Documentation

◆ main()

int main ( )

Definition at line 286 of file exprtk_miller_rabin_primality_test.cpp.

287{
288 miller_rabin_probabilistic_primality_test<double>();
289 return 0;
290}

◆ miller_rabin_probabilistic_primality_test()

template<typename T >
void miller_rabin_probabilistic_primality_test ( )

Definition at line 56 of file exprtk_miller_rabin_primality_test.cpp.

57{
58 typedef exprtk::symbol_table<T> symbol_table_t;
59 typedef exprtk::expression<T> expression_t;
60 typedef exprtk::parser<T> parser_t;
61 typedef exprtk::function_compositor<T> compositor_t;
62 typedef typename compositor_t::function function_t;
63
64 T prime_numbers[] =
65 {
66 5083823, 5083879, 5083889, 5083909, 5083913, 5083927, 5083931, 5083957,
67 5083973, 5083993, 5084017, 5084033, 5084069, 5084099, 5084111, 5084113,
68 5084117, 5084129, 5084141, 5084147, 5084171, 5084179, 5084197, 5084213,
69 5084221, 5084243, 5084249, 5084263, 5084269, 5084333, 5084371, 5084383,
70 5084393, 5084399, 5084407, 5084423, 5084437, 5084447, 5084461, 5084473,
71 5084477, 5084489, 5084501, 5084509, 5084531, 5084537, 5084543, 5084549,
72 5084551, 5084557, 5084567, 5084593, 5084617, 5084627, 5084641, 5084671,
73 5084689, 5084693, 5084711, 5084731, 5084743, 5084749, 5084753, 5084803,
74 5084809, 5084813, 5084843, 5084851, 5084867, 5084869, 5084897, 5084903,
75 6782807, 6782821, 6782827, 6782833, 6782869, 6782903, 6782911, 6782921,
76 6782929, 6782933, 6782939, 6782959, 6782969, 6782977, 6783011, 6783031,
77 6783047, 6783067, 6783083, 6783089, 6783097, 6783107, 6783131, 6783143,
78 6783149, 6783151, 6783169, 6783193, 6783211, 6783223, 6783239, 6783241,
79 6783253, 6783263, 6783299, 6783317, 6783331, 6783353, 6783421, 6783433,
80 6783449, 6783457, 6783461, 6783467, 6783481, 6783521, 6783529, 6783533,
81 9723667, 9723691, 9723697, 9723709, 9723713, 9723757, 9723761, 9723767,
82 9723773, 9723797, 9723799, 9723803, 9723809, 9723821, 9723829, 9723841,
83 9723851, 9723853, 9723869, 9723887, 9723893, 9723907, 9723913, 9723937,
84 9723941, 9723943, 9723979, 9723997, 9724021, 9724061, 9724079, 9724087,
85 9724097, 9724103, 9724129, 9724151, 9724171, 9724201, 9724207, 9724213,
86 9724223, 9724241, 9724259, 9724277, 9724279, 9724283, 9724291, 9724327,
87 10619129, 10619131, 10619137, 10619159, 10619177, 10619179, 10619209, 10619243,
88 10619251, 10619263, 10619291, 10619303, 10619327, 10619341, 10619381, 10619393,
89 10619431, 10619489, 10619501, 10619503, 10619513, 10619563, 10619573, 10619591,
90 10619599, 10619617, 10619621, 10619641, 10619647, 10619669, 10619681, 10619711,
91 10619723, 10619729, 10619731, 10619737, 10619759, 10619773, 10619783, 10619801,
92 10619863, 10619893, 10619897, 10619899, 10619911, 10619923, 10619927, 10619957,
93 10619971, 10619989, 10620007, 10620023, 10620041, 10620067, 10620131, 10620133,
94 10620139, 10620143, 10620149, 10620151, 10620157, 10620163, 10620179, 10620199,
95 10748671, 10748687, 10748693, 10748701, 10748707, 10748713, 10748747, 10748767,
96 10748779, 10748783, 10748791, 10748807, 10748809, 10748827, 10748849, 10748863,
97 10748867, 10748869, 10748873, 10748893, 10748897, 10748917, 10748921, 10748923,
98 10748951, 10748953, 10748957, 10748963, 10748987, 10748999, 10749001, 10749023,
99 10749029, 10749059, 10749061, 10749071, 10749077, 10749107, 10749121, 10749143,
100 10749161, 10749163, 10749173, 10749191, 10749197, 10749199, 10749241, 10749251,
101 10749259, 10749283, 10749313, 10749317, 10749337, 10749367, 10749371, 10749373,
102 10749379, 10749391, 10749433, 10749437, 10749449, 10749461, 10749493, 10749581,
103 10749587, 10749601, 10749617, 10749631, 10749659, 10749677, 10749679, 10749703,
104 10749709, 10749719, 10749727, 10749743, 10749773, 10749779, 10749811, 10749821,
105 10749833, 10749853, 10749857, 10749863, 10749881, 10749883, 10749889, 10749917,
106 10749971, 10749983, 10749989, 10749997, 10750009, 10750021, 10750031, 10750037,
107 10750049, 10750057, 10750073, 10750099, 10750127, 10750141, 10750147, 10750189,
108 10750193, 10750249, 10750279, 10750303, 10750309, 10750319, 10750349, 10750373,
109 };
110
111 T composites[] =
112 {
113 199203677, 779234623, 843093203, 883543291, 1197162971,
114 1282615157, 1552390397, 1765737859, 1878769589, 1993904873,
115 2257133471, 2520523529, 2579094799, 2853450949, 2935025369,
116 3095780533, 3164132249, 3408963511, 4260042859, 4608613981,
117 4654875857, 5085931997, 7278175081, 7289187463, 9206112101
118 };
119
120 exprtk::rtl::io::println<T> println;
121 random_uint<T> random;
122
123 symbol_table_t symbol_table;
124
125 symbol_table.add_vector ("prime_numbers" , prime_numbers);
126 symbol_table.add_vector ("composite_numbers", composites );
127 symbol_table.add_function("println" , println );
128 symbol_table.add_function("random" , random );
129
130 compositor_t compositor(symbol_table);
131
132 compositor.add(function_t()
133 .name("modexp") // (base^exponent) mod m
134 .var("b") // base
135 .var("e") // exponent
136 .var("m") // modulo
137 .expression
138 (
139 " if (m == 1) "
140 " { "
141 " return [0]; "
142 " }; "
143 " "
144 " var result := 1; "
145 " b %= m; "
146 " "
147 " while (e > 0) "
148 " { "
149 " if ((e % 2) == 1) "
150 " { "
151 " result := (result * b) % m; "
152 " }; "
153 " "
154 " b := b^2 % m; "
155 " e := floor(e / 2); "
156 " }; "
157 " "
158 " result "
159 ));
160
161 compositor.add(function_t()
162 .name("is_probable_prime")
163 .var("n") // Number to test for primaility
164 .var("k") // Number of iterations to perform
165 .expression
166 (
167 " switch "
168 " { "
169 " case n <= 1 : return [false ]; "
170 " case frac(n) != 0 : return [false ]; "
171 " case n == 2 : return [true ]; "
172 " }; "
173 " "
174 " var primes[18] := "
175 " { "
176 " 2, 3, 5, 7, 11, 13, 17, 19, 23, "
177 " 29, 31, 37, 41, 43, 47, 53, 59, 61 "
178 " }; "
179 " "
180 " var upper_bound := min(n - 1, trunc(sqrt(n)) + 1); "
181 " "
182 " for (var i := 0; i < primes[]; i += 1) "
183 " { "
184 " if (primes[i] >= upper_bound) "
185 " return [true]; "
186 " else if ((n % primes[i]) == 0) "
187 " return [false]; "
188 " }; "
189 " "
190 " var s := n - 1; "
191 " var t := 0; "
192 " "
193 " while ((s % 2) == 0) "
194 " { "
195 " s := floor(s / 2); "
196 " t += 1; "
197 " }; "
198 " "
199 " for (var trials := 0; trials < k; trials += 1) "
200 " { "
201 " var a := random(2, n - 2); "
202 " var x := modexp(a, s, n); "
203 " "
204 " if ((x == 1) or (x == (n - 1))) "
205 " { "
206 " continue; "
207 " }; "
208 " "
209 " for (var r := 1; r < t; r += 1) "
210 " { "
211 " x := modexp(x, 2, n); "
212 " "
213 " if (x == 1) "
214 " return [false]; "
215 " else if (x == (n - 1)) "
216 " break; "
217 " }; "
218 " "
219 " if (x != (n - 1)) "
220 " { "
221 " return [false]; "
222 " }; "
223 " }; "
224 " "
225 " return [true]; "
226 ));
227
228 const std::string miller_rabin_primality_test_program =
229 " "
230 " println('Testing prime numbers...'); "
231 " for (var i := 0; i < prime_numbers[]; i += 1) "
232 " { "
233 " var n := prime_numbers[i]; "
234 " "
235 " if (is_probable_prime(n,1000)) "
236 " { "
237 " println(n, ' is correctly a prime number'); "
238 " } "
239 " else "
240 " { "
241 " println(n, ' is NOT a prime number (error)'); "
242 " } "
243 " }; "
244 " "
245 " println('Testing composite numbers...'); "
246 " for (var i := 0; i < composite_numbers[]; i += 1) "
247 " { "
248 " var n := composite_numbers[i]; "
249 " "
250 " if (not(is_probable_prime(n,1000))) "
251 " { "
252 " println(n, ' is correctly a composite'); "
253 " } "
254 " else "
255 " { "
256 " println(n, ' is a prime number (error)'); "
257 " } "
258 " }; "
259 " "
260 " for (var i := 0; i < prime_numbers[]; i += 1) "
261 " { "
262 " for (var j := i; j < prime_numbers[]; j += 1) "
263 " { "
264 " var n := prime_numbers[i] * prime_numbers[j]; "
265 " "
266 " if (not(is_probable_prime(n,1000))) "
267 " { "
268 " println(n, ' is correctly a composite'); "
269 " } "
270 " else "
271 " { "
272 " println(n, ' is a prime number (error)'); "
273 " } "
274 " } "
275 " }; ";
276
277 expression_t expression;
278 expression.register_symbol_table(symbol_table);
279
280 parser_t parser;
281 parser.compile(miller_rabin_primality_test_program,expression);
282
283 expression.value();
284}
exprtk::function_compositor
Definition exprtk.hpp:42434
exprtk::symbol_table
Definition exprtk.hpp:19745
random_uint
Definition exprtk_miller_rabin_primality_test.cpp:31
Generated by doxygen 1.9.8