1 #include "cado.h" // IWYU pragma: keep
2 // IWYU pragma: no_include <ext/alloc_traits.h>
3 #include <cfloat>
4 #include <cstdlib>
5 #include <cstdio>
6 #include <cmath>
7 #include <climits> // for ULONG_MAX
8 #include <memory> // for allocator_traits<>::value_type
9 #include <vector> // for vector
10 #include <gmp.h> // for gmp_randstate_t, gmp_randclear
11 #include "cxx_mpz.hpp"
12 #include "decomp.h" // for decomp_t, decomp_create
13 #include "facul.hpp"
14 #include "facul_ecm.h"
15 #include "facul_fwd.hpp" // for facul_method_t
16 #include "fm.h" // for fm_t, fm_create, fm_free
17 #include "generate_factoring_method.hpp"
18 #include "generate_strategies.h"
19 #include "macros.h"
20 #include "modredc_15ul.h" // for MODREDC15UL_MAXBITS
21 #include "modredc_2ul2.h" // for MODREDC2UL2_MAXBITS
22 #include "modredc_ul.h" // for MODREDCUL_MAXBITS
23 #include "pm1.h"
24 #include "pp1.h"
25 #include "strategy.h" // for strategy_t, strategy_add_fm
26 #include "tab_decomp.h" // for tabular_decomp_t, tabular_d...
27 #include "tab_fm.h" // for tabular_fm_t, tabular_fm_ad...
28 #include "tab_strategy.h" // for tabular_strategy_t, tabular...
29 #include "timing.h" // microseconds
30
31 static double EPSILON_DBL = LDBL_EPSILON;
32
33 //convert type: from strategy_t to facul_strategy_t.
convert_strategy_to_facul_strategy(strategy_t * t)34 facul_strategy_t* convert_strategy_to_facul_strategy (strategy_t* t)
35 {
36 tabular_fm_t* tab_fm = strategy_get_tab_fm (t);
37 int nb_methods = tab_fm->index;
38
39 facul_strategy_t *strategy;
40 strategy = (facul_strategy_t*) malloc(sizeof(facul_strategy_t));
41 strategy->methods = (facul_method_t*) malloc((nb_methods+1) * sizeof(facul_method_t));
42 /*
43 without this second method, the function
44 facul_clear_strategy (strategy) failed!
45 */
46 strategy->methods[1].method = 0;
47 strategy->methods[1].plan = NULL;
48
49 strategy->lpb = ULONG_MAX;
50 strategy->assume_prime_thresh = 0;
51 strategy->BBB = 0;
52
53 for (int i = 0; i < nb_methods; i++)
54 {
55 fm_t* fm = tab_fm->tab[i];
56 int method = (int)fm->method[0];
57 ec_parameterization_t curve = (ec_parameterization_t) fm->method[1];
58 unsigned long B1 = fm->method[2];
59 unsigned long B2 = fm->method[3];
60 strategy->methods[i].method = method;
61
62 if (method == PM1_METHOD) {
63 strategy->methods[i].plan = malloc(sizeof(pm1_plan_t));
64 ASSERT(strategy->methods[i].plan != NULL);
65 pm1_make_plan((pm1_plan_t*) strategy->methods[i].plan, B1, B2, 0);
66 } else if (method == PP1_27_METHOD || method == PP1_65_METHOD) {
67 strategy->methods[i].plan = malloc(sizeof(pp1_plan_t));
68 ASSERT(strategy->methods[i].plan != NULL);
69 pp1_make_plan((pp1_plan_t*) strategy->methods[i].plan, B1, B2, 0);
70 } else if (method == EC_METHOD) {
71 unsigned long sigma;
72 if (curve == MONTY16) {
73 sigma = 1;
74 } else
75 sigma = 2 + rand();
76
77 strategy->methods[i].plan = malloc(sizeof(ecm_plan_t));
78 ASSERT((ecm_plan_t*) strategy->methods[i].plan != NULL);
79 ecm_make_plan((ecm_plan_t*) strategy->methods[i].plan, B1, B2, curve,
80 sigma, 0, 0);
81 } else {
82 exit(EXIT_FAILURE);
83 }
84 }
85 strategy->methods[nb_methods].method = 0;
86 strategy->methods[nb_methods].plan = NULL;
87
88 return strategy;
89 }
90
91 static int
select_random_index_dec(double sum_nb_elem,tabular_decomp_t * t,gmp_randstate_ptr rstate)92 select_random_index_dec(double sum_nb_elem, tabular_decomp_t* t, gmp_randstate_ptr rstate)
93 {
94 double alea = gmp_urandomm_ui(rstate, (unsigned long)sum_nb_elem);
95 int i = 0;
96 double bound = t->tab[0]->nb_elem;
97 int len = t->index;
98 while (i < (len-1) && (alea - bound) >= 1) {
99 i++;
100 bound += t->tab[i]->nb_elem;
101 }
102 return i;
103 }
104
105
106
107 /*
108 This function compute directly the probabity and the time to find a
109 factor in a good decomposition in a cofactor of length r.
110 */
111 weighted_success
bench_proba_time_st(gmp_randstate_t state,facul_strategy_t * strategy,tabular_decomp_t * init_tab,int r,int lpb)112 bench_proba_time_st(gmp_randstate_t state, facul_strategy_t* strategy,
113 tabular_decomp_t* init_tab, int r, int lpb)
114 {
115 int nb_succes_max = 10000, nb_succes = 0;
116 int nb_test = 0;
117 int nb_test_max = 10 * nb_succes_max; //20 for five percent
118 double time = 0;
119
120 double sum_dec = 0;
121 for (int i = 0; i < init_tab->index; i++)
122 sum_dec += init_tab->tab[i]->nb_elem;
123 //struct timeval st_test, end_test;
124
125 std::vector<cxx_mpz> f;
126
127 while (nb_succes < nb_succes_max && (nb_test<nb_test_max))
128 {
129 int index = select_random_index_dec(sum_dec, init_tab, state);
130 int len_p = init_tab->tab[index]->tab[0];
131
132 cxx_mpz N = generate_composite_integer(state, len_p, r);
133 f.clear();
134
135 double starttime = microseconds();
136 //gettimeofday (&st_test, NULL);
137 int nfound = facul(f, N, strategy);
138 //gettimeofday (&end_test, NULL);
139 double endtime = microseconds();
140
141 /* time += (end_test.tv_sec - st_test.tv_sec)*1000000L */
142 /* + end_test.tv_usec - st_test.tv_usec; */
143 time += endtime - starttime;
144 if (nfound != 0) {
145 int len_factor = mpz_sizeinbase(f[0], 2);
146 if (len_factor <= lpb && (r-len_factor) <= lpb)
147 nb_succes++;
148 }
149 nb_test++;
150 }
151
152 return weighted_success(nb_succes, time, nb_test);
153 }
154
155
156 int
get_nb_word(int r)157 get_nb_word (int r)
158 {
159 /*
160 We add 0.5 to the length of one word, because for our times the
161 lenght is inclusive. For example, if MODREDCUL_MAXBITS = 64
162 bits, a cofactor is in one word if is lenght is less OR equal to
163 64 bits. So, if you don't add 0.5 to MODREDCUL_MAXBITS, you
164 lost the equal and thus insert an error in your maths.
165 */
166 double half_word = (MODREDCUL_MAXBITS+0.5)/2.0;
167 int number_half_wd = floor(r /half_word);
168 int ind = (number_half_wd <2)? 0: number_half_wd - 1;
169 return ind;
170 }
171
172 /*
173 This function is like bench_time() but only compute the time for one
174 length. In fact, i remove the unnecessary computation to reduce the
175 cost of this binary.
176 */
bench_time_mini(gmp_randstate_t state,tabular_fm_t * fm,int r)177 void bench_time_mini(gmp_randstate_t state, tabular_fm_t * fm, int r)
178 {
179 unsigned long *param;
180 fm_t *elem;
181 int len = fm->index; //number of methods!
182 //precompute 4 arrays for our bench_time!
183 //{{
184 int len_n;
185 if (r <= MODREDCUL_MAXBITS)
186 len_n = MODREDCUL_MAXBITS;
187 else if (r <= MODREDC15UL_MAXBITS)
188 len_n = MODREDC15UL_MAXBITS;
189 else if (r <= MODREDC2UL2_MAXBITS)
190 len_n = MODREDC2UL2_MAXBITS;
191 else
192 len_n = MODREDC2UL2_MAXBITS +30;
193 int nb_test = 100000;
194 std::vector<cxx_mpz> N;
195 N.reserve(nb_test);
196 for (int i = 0; i < nb_test; i++) {
197 N.push_back(generate_prime_factor(state, len_n));
198 }
199 //}}
200 int ind = get_nb_word(r);
201 for (int i = 0; i < len; i++) {
202 elem = tabular_fm_get_fm(fm, i);
203 param = fm_get_method(elem);
204 unsigned long method = param[0];
205 ec_parameterization_t curve = (ec_parameterization_t) param[1];
206 unsigned long B1 = param[2];
207 unsigned long B2 = param[3];
208 if (B1 != 0 || B2 != 0) {
209 facul_strategy_t *st = generate_fm(method, curve, B1, B2);
210
211 double res[4];
212 ASSERT_ALWAYS(ind < 4);
213 res[ind] = bench_time_fm_onelength(st, N, nb_test);
214 fm_set_time(elem, res, 4);
215 facul_clear_strategy(st);
216 } else {
217 double time[4] = { 0, 0, 0, 0 };
218 fm_set_time(elem, time, 4);
219 }
220
221 }
222 }
223 /*
224 This function is like bench_proba() but only compute the necessary
225 probabilities. In fact, i remove the unnecessary computation to
226 reduce the cost of this binary.
227 */
228
bench_proba_mini(gmp_randstate_t state,tabular_fm_t * fm,int * val_p,int len_val_p,int len_p_min)229 void bench_proba_mini(gmp_randstate_t state, tabular_fm_t * fm, int* val_p,
230 int len_val_p, int len_p_min)
231 {
232 int len = fm->index; //number of methods!
233 int p_max = 100;
234 double *proba = (double*) calloc(p_max, sizeof(double));
235 ASSERT(proba != NULL);
236
237 unsigned long *param;
238 fm_t *elem;
239
240 //{{Will contain the precomputes of composite integer!
241 std::vector<std::vector<cxx_mpz>> N(len_val_p);
242 int nb_test_max = 10000;
243
244 //}}
245
246 for (int i = 0; i < len; i++) {
247 elem = tabular_fm_get_fm(fm, i);
248 param = fm_get_method(elem);
249 unsigned long method = param[0];
250 ec_parameterization_t curve = (ec_parameterization_t) param[1];
251 unsigned long B1 = param[2];
252 unsigned long B2 = param[3];
253
254 facul_strategy_t *st = generate_fm(method, curve, B1, B2);
255
256 int max_index = 0;
257 for (int j = 0; j < len_val_p; j++)
258 {
259 int ind_proba = val_p[j] - len_p_min;
260 if (B1 == 0 && B2 == 0)
261 proba[ind_proba] = 0;
262 else {
263 int len_p = len_p_min + ind_proba;
264 int len_n = 60 + len_p;
265 proba[ind_proba] = bench_proba_fm(st, state, len_p,
266 len_n, N[j],
267 nb_test_max);
268 }
269 if (ind_proba > max_index)
270 max_index = ind_proba;
271 }
272 fm_set_proba(elem, proba, max_index+1, len_p_min);
273 //free
274 facul_clear_strategy(st);
275 }
276 free(proba);
277 }
278
279
280 /************************************************************************/
281 /* MAIN */
282 /************************************************************************/
283
284
285 // coverity[root_function]
main()286 int main()
287 {
288 gmp_randstate_t state;
289 gmp_randinit_default(state);
290
291 int fbb = 15;//25;
292 int lpb = 19;//29;
293 int r = 35;//55;
294
295 //input decomp!
296 tabular_decomp_t* init_tab = tabular_decomp_create ();
297 int tmp1[2] = {17,18};
298 decomp_t* el1 = decomp_create (2, tmp1, 100);
299 tabular_decomp_add_decomp (init_tab, el1);
300 decomp_free (el1);
301 int tmp2[2] = {15,20};
302 decomp_t* el2 = decomp_create (2, tmp2, 200);
303 tabular_decomp_add_decomp (init_tab, el2);
304 decomp_free (el2);
305
306 int len_val_factor = 4;
307 int val_factor[4] = {15, 17, 18, 20};
308
309 //fm
310 fm_t* pm1 = fm_create ();
311 unsigned long elem1[4] = {PM1_METHOD, 0, 50, 500};
312 fm_set_method (pm1, elem1, 4);
313 fm_t* pp1 = fm_create ();
314 unsigned long elem2[4] = {PP1_65_METHOD, 0, 70, 700};
315 fm_set_method (pp1, elem2, 4);
316 fm_t* ecm = fm_create ();
317 unsigned long elem3[4] = {EC_METHOD, BRENT12, 80, 1000};
318 fm_set_method (ecm, elem3, 4);
319 tabular_fm_t* tab = tabular_fm_create();
320 tabular_fm_add_fm (tab, pm1);
321 tabular_fm_add_fm (tab, pp1);
322 tabular_fm_add_fm (tab, ecm);
323 //bench our method
324 bench_proba_mini(state, tab, val_factor, len_val_factor, fbb);
325 //bench_time_mini (state, tab, r);
326 /* printf ("mini\n"); */
327 /* tabular_fm_print (tab); */
328
329 /* bench_proba (state, tab, fbb); */
330 /* bench_time(state, tab); */
331 /* printf ("\n all \n"); */
332 /* tabular_fm_print (tab); */
333 //generate some examples of strategies!
334 strategy_t* strat1 =strategy_create ();
335 strategy_add_fm (strat1, tab->tab[0]);//pm1
336 strategy_add_fm (strat1, tab->tab[1]);//pp1
337 strategy_add_fm (strat1, tab->tab[2]);//ecm
338 double prob1 = compute_proba_strategy(init_tab, strat1, fbb, lpb);
339 //double time1 = compute_time_strategy(init_tab, strat1, r);
340
341 //Create our strategy to use facul().
342 facul_strategy_t* st = convert_strategy_to_facul_strategy (strat1);
343 //bench our strategies!
344 weighted_success res = bench_proba_time_st(state, st, init_tab, r, lpb);
345 double prob2 = res.prob;
346 //double time2 = res[1];
347
348 printf ("prob1 = %lf, prob2 = %lf\n", prob1, prob2);
349 // printf ("time1 = %lf, time2 = %lf\n", time1, time2);
350
351 double precision_p = 0.05;
352 if ((prob1 - prob2) > precision_p || (prob2 - prob1) > precision_p)
353 {
354 fprintf (stderr, "error with the test(1)\n");
355 return EXIT_FAILURE;
356 }
357
358 //{{test the function: bench_proba_time_pset()
359 int c = tab->tab[0]->method[3]/tab->tab[0]->method[2];
360 int param[6] = {(int) tab->tab[0]->method[2], (int) tab->tab[0]->method[2], 1, c,c,1};
361
362 tabular_fm_t* tmp =
363 bench_proba_time_pset (tab->tab[0]->method[0], (ec_parameterization_t) tab->tab[0]->method[1], state,
364 17, 20, 18*3, param);
365
366 /*check this probability: the probability to find a prime number
367 of length 17 or 18 bits must be more than this to find 18 bits and
368 less than this for 18 bits.*/
369 fm_t* pm1_bis = tmp->tab[1];
370 /* fm_print (tab->tab[0]); */
371 /* fm_print (tmp->tab[1]); */
372 if (pm1_bis->proba[0] < tab->tab[0]->proba[20-fbb] ||
373 pm1_bis->proba[0] > tab->tab[0]->proba[17-fbb])
374 {
375 fprintf (stderr, "error with the test(2)\n");
376 return EXIT_FAILURE;
377 }
378 tabular_fm_free (tmp);
379
380 //}}
381
382 //test the generation of our strategy for one pair of cofactor!!!
383 //{{
384 strategy_set_proba(strat1, 0.4);
385 strategy_set_time(strat1, 40);
386 tabular_strategy_t * strat_r0 = tabular_strategy_create ();
387 tabular_strategy_t * strat_r1 = tabular_strategy_create ();
388
389 //firstly, create the zero strategy!
390 strategy_t* zero_st = strategy_create();
391 unsigned long tab0[4] = {PM1_METHOD, 0, 0, 0};
392 fm_t* zero_fm = fm_create();
393 fm_set_method (zero_fm, tab0, 4);
394 strategy_add_fm (zero_st, zero_fm);
395 strategy_set_proba(zero_st, 0.0);
396
397 tabular_strategy_add_strategy (strat_r0, zero_st);
398 tabular_strategy_add_strategy (strat_r0, strat1);
399 tabular_strategy_add_strategy (strat_r1, zero_st);
400 tabular_strategy_t * res2 = generate_strategy_r0_r1(strat_r0, strat_r1);
401
402 //check proba + time!
403 if (!(res2->index == 1 && res2->tab[0]->proba < EPSILON_DBL
404 && res2->tab[0]->time < EPSILON_DBL))
405 {
406 fprintf (stderr, "error with the test(3)\n");
407 return EXIT_FAILURE;
408 }
409 strategy_set_proba(strat_r1->tab[0], 1.0);
410 tabular_strategy_t * res3 = generate_strategy_r0_r1(strat_r0, strat_r1);
411
412 //check proba + time!
413 if (!(res3->index == 2 && res3->tab[0]->proba < EPSILON_DBL
414 && res3->tab[0]->time < EPSILON_DBL &&
415 (res3->tab[1]->proba - strat1->proba) < EPSILON_DBL &&
416 (res3->tab[1]->time - strat1->time) < EPSILON_DBL))
417 {
418 fprintf (stderr, "error with the test(3)\n");
419 return EXIT_FAILURE;
420 }
421
422 //}}
423 //free
424 facul_clear_strategy(st);
425 strategy_free (strat1);
426 tabular_fm_free (tab);
427 fm_free (pm1);
428 fm_free (pp1);
429 fm_free (ecm);
430 gmp_randclear (state);
431 tabular_decomp_free (init_tab);
432
433 fm_free (zero_fm);
434 strategy_free (zero_st);
435 tabular_strategy_free (res2);
436 tabular_strategy_free (res3);
437 tabular_strategy_free (strat_r0);
438 tabular_strategy_free (strat_r1);
439
440
441 return EXIT_SUCCESS;
442 }
443