C++ Mathematical Expression Toolkit (ExprTk) release
Loading...
Searching...
No Matches
exprtk_miller_rabin_primality_test.cpp
Go to the documentation of this file.
1/*
2 **************************************************************
3 * C++ Mathematical Expression Toolkit Library *
4 * *
5 * ExprTk Miller-Rabin Probabilistic Primality Test *
6 * Author: Arash Partow (1999-2024) *
7 * URL: https://www.partow.net/programming/exprtk/index.html *
8 * *
9 * Copyright notice: *
10 * Free use of the Mathematical Expression Toolkit Library is *
11 * permitted under the guidelines and in accordance with the *
12 * most current version of the MIT License. *
13 * https://www.opensource.org/licenses/MIT *
14 * SPDX-License-Identifier: MIT *
15 * *
16 **************************************************************
17*/
18
19
20#include <algorithm>
21#include <array>
22#include <cstdio>
23#include <random>
24#include <string>
25
26#include "exprtk.hpp"
27
28
29template <typename T>
30struct random_uint final : public exprtk::ifunction<T>
31{
32 using exprtk::ifunction<T>::operator();
33
34 random_uint()
35 : exprtk::ifunction<T>(2)
36 {
37 std::random_device device;
38 std::array<unsigned int,std::mt19937::state_size> seed;
39 std::generate_n(seed.data(), seed.size(), std::ref(device));
40 std::seed_seq seq(std::begin(seed), std::end(seed));
41 generator.seed(seq);
42 }
43
44 inline T operator() (const T& r0,const T& r1) override
45 {
46 std::uniform_int_distribution<std::uint64_t>
47 distribution { static_cast<std::uint64_t>(r0), static_cast<std::uint64_t>(r1) };
48
49 return static_cast<T>(distribution(generator));
50 }
51
52 std::mt19937 generator;
53};
54
55template <typename T>
56void miller_rabin_probabilistic_primality_test()
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}
285
286int main()
287{
288 miller_rabin_probabilistic_primality_test<double>();
289 return 0;
290}
exprtk::function_compositor
Definition exprtk.hpp:42434
exprtk::ifunction::operator()
virtual T operator()()
Definition exprtk.hpp:19558
exprtk::ifunction::ifunction
ifunction(const std::size_t &pc)
Definition exprtk.hpp:19545
exprtk::symbol_table
Definition exprtk.hpp:19745
miller_rabin_probabilistic_primality_test
void miller_rabin_probabilistic_primality_test()
Definition exprtk_miller_rabin_primality_test.cpp:56
main
int main()
Definition exprtk_miller_rabin_primality_test.cpp:286
random_uint
Definition exprtk_miller_rabin_primality_test.cpp:31
random_uint::generator
std::mt19937 generator
Definition exprtk_miller_rabin_primality_test.cpp:52
random_uint::random_uint
random_uint()
Definition exprtk_miller_rabin_primality_test.cpp:34
Generated by doxygen 1.9.8