1 /*
2 Copyright (C) 2006, 2011, 2016 William Hart
3 Copyright (C) 2015 Nitin Kumar
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <https://www.gnu.org/licenses/>.
11 */
12
13 #ifndef QSIEVE_H
14 #define QSIEVE_H
15
16 #include <stdio.h>
17 #include <math.h>
18 #include <stdlib.h>
19 #include <string.h>
20 #include <stdint.h>
21
22 #include "flint.h"
23 #include "fmpz.h"
24 #include "ulong_extras.h"
25 #include "fmpz_vec.h"
26 #include "fmpz_factor.h"
27 #include "thread_support.h"
28
29 #ifdef __cplusplus
30 extern "C" {
31 #endif
32
33 #define QS_DEBUG 0 /* level of debug information printed, 0 = none */
34
35 #define BITS_ADJUST 25 /* add to sieve entries to compensate approximations */
36
37 #define BLOCK_SIZE (4*65536) /* size of sieving cache block */
38
39 typedef struct prime_t
40 {
41 mp_limb_t pinv; /* precomputed inverse */
42 int p; /* prime */
43 char size;
44 } prime_t;
45
46 typedef struct fac_t /* struct for factors of relations */
47 {
48 slong ind;
49 slong exp;
50 } fac_t;
51
52 typedef struct la_col_t /* matrix column */
53 {
54 slong * data; /* The list of occupied rows in this column */
55 slong weight; /* Number of nonzero entries in this column */
56 slong orig; /* Original relation number */
57 } la_col_t;
58
59
60 typedef struct hash_t /* entry in hash table */
61 {
62 mp_limb_t prime; /* value of prime */
63 mp_limb_t next; /* next prime which have same hash value as 'prime' */
64 mp_limb_t count; /* number of occurrence of 'prime' */
65 } hash_t;
66
67 typedef struct relation_t /* format for relation */
68 {
69 mp_limb_t lp; /* large prime, is 1, if relation is full */
70 slong num_factors; /* number of factors, excluding small factor */
71 slong small_primes; /* number of small factors */
72 slong * small; /* exponent of small factors */
73 fac_t * factor; /* factor of relation */
74 fmpz_t Y; /* square root of sieve value for relation */
75 } relation_t;
76
77 typedef struct qs_poly_s
78 {
79 fmpz_t B; /* current B coeff of poly */
80 int * soln1; /* first start position in sieve per prime */
81 int * soln2; /* second start position in sieve per prime */
82 int * posn1; /* temp space for sieving */
83 int * posn2; /* temp space for sieving */
84 slong * small; /* exponents of small prime factors in relations */
85 fac_t * factor; /* factors for a relation */
86 slong num_factors; /* number of factors found in a relation */
87 } qs_poly_s;
88
89 typedef qs_poly_s qs_poly_t[1];
90
91 typedef struct qs_s
92 {
93 volatile slong index_j;
94 #if FLINT_USES_PTHREAD
95 pthread_mutex_t mutex;
96 #endif
97 thread_pool_handle * handles;
98 slong num_handles;
99
100 fmpz_t n; /* Number to factor */
101
102 flint_bitcnt_t bits; /* Number of bits of n */
103
104 ulong ks_primes; /* number of Knuth-Schroeppel primes */
105
106 mp_limb_t k; /* Multiplier */
107 fmpz_t kn; /* kn as a multiprecision integer */
108
109 slong num_primes; /* number of factor base primes including k and 2 */
110
111 prime_t * factor_base; /* data about factor base primes */
112
113 int * sqrts; /* square roots of kn mod factor base primes */
114
115 slong small_primes; /* number of primes to not sieve with */
116 slong second_prime; /* index of first prime bigger than block size */
117 slong sieve_size; /* size of sieve to use */
118
119 unsigned char sieve_bits; /* sieve threshold */
120 unsigned char sieve_fill; /* for biasing sieve values */
121
122 /***************************************************************************
123 POLYNOMIAL DATA
124 **************************************************************************/
125
126 fmpz_t A; /* current value of coeff A of poly Ax^2 + Bx + C */
127 fmpz_t B; /* B value of poly */
128
129 mp_limb_t * A_ind; /* indices of factor base primes dividing A */
130
131 fmpz_t * A_divp; /* A_divp[i] = (A/p_i),
132 where the p_i are the prime factors of A */
133 mp_limb_t * B0_terms; /* B0_terms[i] = min(gamma_i, p - gamma_i) where
134 gamma_i = (sqrt(kn)*(A_divp[i])^(-1)) mod p_i,
135 where the p_i are the prime factors of A */
136
137 fmpz_t * B_terms; /* B_terms[i] = A_divp[i]*B0_terms[i] (multprec) */
138
139 mp_limb_t * A_inv; /* A_inv[k] = A^(-1) mod p_k, for FB prime p_k */
140 mp_limb_t ** A_inv2B; /* A_inv2B[i][k] = 2 * B_terms[i] * A^(-1) mod p_k
141 for FB prime p_k */
142
143 int * soln1; /* soln1[k] = first poly root mod FB prime p_k */
144 int * soln2; /* soln2[k] = second poly root mod FB prime p_k */
145
146 fmpz_t target_A; /* approximate target value for A coeff of poly */
147
148 fmpz_t upp_bound; /* upper bound on desired A = 2*target_A */
149 fmpz_t low_bound; /* lower bound on desired A = target_A/2 */
150
151 slong s; /* number of prime factors of A */
152 slong low; /* minimum offset in factor base,
153 for possible factors of A */
154 slong high; /* maximum offset in factor base,
155 for possible factors of A */
156 slong span; /* total number of possible factors for A */
157
158 /*
159 parameters for calculating next subset of possible factors of A, giving
160 a lexicographic ordering of all such tuples
161 */
162 slong h; /* tuple entry we just set, numbered from 1 at end of tuple */
163 slong m; /* last value we just set a tuple entry to */
164 slong A_ind_diff; /* diff. between indices of (s-1) and (s-2)-th A-factor */
165 mp_limb_t * curr_subset; /* current tuple */
166 mp_limb_t * first_subset; /* first tuple, in case of restart */
167 mp_limb_t j; /* index of s-th factor of first A, if s > 3 */
168
169 #if QS_DEBUG
170 slong poly_count; /* keep track of the number of polynomials used */
171 #endif
172
173 qs_poly_s * poly; /* poly data per thread */
174
175 /***************************************************************************
176 RELATION DATA
177 ***************************************************************************/
178
179 FILE * siqs; /* pointer to file for storing relations */
180 char * fname; /* name of file used for relations */
181
182 slong full_relation; /* number of full relations */
183 slong num_cycles; /* number of possible full relations from partials */
184
185 slong vertices; /* number of different primes in partials */
186 slong components; /* equal to 1 */
187 slong edges; /* total number of partials */
188
189 slong table_size; /* size of table */
190 hash_t * table; /* store 'prime' occurring in partial */
191 mp_limb_t * hash_table; /* to keep track of location of primes in 'table' */
192
193 slong extra_rels; /* number of extra relations beyond num_primes */
194 slong max_factors; /* maximum number of factors a relation can have */
195
196 fmpz * Y_arr; /* array of Y's corresponding to relations */
197 slong * curr_rel; /* current relation in array of relations */
198 slong * relation; /* relation array */
199
200 slong buffer_size; /* size of buffer of relations */
201 slong num_relations; /* number of relations so far */
202
203 ulong small_factor; /* small factor found when merging relations */
204
205 /***************************************************************************
206 LINEAR ALGEBRA DATA
207 ***************************************************************************/
208
209 la_col_t * matrix; /* the main matrix over GF(2) in sparse format */
210 la_col_t ** qsort_arr; /* array of columns ready to be sorted */
211
212 slong columns; /* number of columns in matrix so far */
213
214 /***************************************************************************
215 SQUARE ROOT DATA
216 ***************************************************************************/
217
218 slong * prime_count; /* counts of the exponents of primes appearing in the square */
219
220 } qs_s;
221
222 typedef qs_s qs_t[1];
223
224 /*
225 Tuning parameters { bits, ks_primes, fb_primes, small_primes, sieve_size}
226 for qsieve_factor_threaded where:
227 * bits is the number of bits of n
228 * ks_primes is the max number of primes to try in Knuth-Schroeppel function
229 * fb_primes is the number of factor base primes to use (including k and 2)
230 * small_primes is the number of small primes to not factor with (including k and 2)
231 * sieve_size is the size of the sieve to use
232 * sieve_bits - sieve_fill
233 */
234
235 #if 0 /* TODO have the tuning values taken from here if multithreaded */
236
237 static const mp_limb_t qsieve_tune[][6] =
238 {
239 {10, 50, 100, 5, 2 * 2000, 30}, /* */
240 {20, 50, 120, 6, 2 * 2500, 30}, /* */
241 {30, 50, 150, 6, 2 * 2000, 31}, /* */
242 {40, 50, 150, 8, 2 * 3000, 32}, /* 12 digits */
243 {50, 50, 150, 8, 2 * 3000, 34}, /* 15 digits */
244 {60, 50, 150, 9, 2 * 3500, 36}, /* 18 digits */
245 {70, 100, 200, 9, 2 * 4000, 42}, /* 21 digits */
246 {80, 100, 200, 9, 2 * 6000, 44}, /* 24 digits */
247 {90, 100, 200, 9, 2 * 6000, 50}, /* */
248 {100, 100, 300, 9, 2 * 7000, 54}, /* */
249 {110, 100, 500, 9, 2 * 25000, 62}, /* 31 digits */
250 {120, 100, 800, 9, 2 * 30000, 64}, /* */
251 {130, 100, 1000, 9, 2 * 30000, 64}, /* 41 digits */
252 {140, 100, 1200, 9, 2 * 30000, 66}, /* */
253 {150, 100, 1500, 10, 2 * 32000, 68}, /* 45 digit */
254 {160, 150, 1800, 11, 2 * 32000, 70}, /* */
255 {170, 150, 2000, 12, 2 * 32000, 72}, /* 50 digits */
256 {180, 150, 2500, 12, 2 * 32000, 73}, /* */
257 {190, 150, 2800, 12, 2 * 32000, 76}, /* */
258 {200, 200, 4000, 12, 2 * 32000, 80}, /* 60 digits */
259 {210, 100, 3600, 12, 2 * 32000, 83}, /* */
260 {220, 300, 6000, 15, 2 * 65536, 87}, /* */
261 {230, 350, 8500, 17, 3 * 65536, 90}, /* 70 digits */
262 {240, 400, 10000, 19, 4 * 65536, 93}, /* */
263 {250, 500, 15000, 19, 4 * 65536, 97}, /* 75 digits */
264 {260, 600, 25000, 25, 4 * 65536, 100}, /* 80 digits */
265 {270, 800, 35000, 27, 5 * 65536, 104} /* */
266 };
267
268 #else /* currently tuned for four threads */
269
270 static const mp_limb_t qsieve_tune[][6] =
271 {
272 {10, 50, 90, 5, 2 * 1500, 18}, /* */
273 {20, 50, 90, 6, 2 * 1600, 18}, /* */
274 {30, 50, 100, 6, 2 * 1800, 19}, /* */
275 {40, 50, 100, 8, 2 * 2000, 20}, /* 13 digits */
276 {50, 50, 100, 8, 2 * 2500, 22}, /* 16 digits */
277 {60, 50, 100, 9, 2 * 3000, 24}, /* 19 digits */
278 {70, 100, 250, 9, 2 * 6000, 25}, /* 22 digits */
279 {80, 100, 250, 9, 2 * 8000, 26}, /* 25 digits */
280 {90, 100, 250, 9, 2 * 9000, 30}, /* 28 digits */
281 {100, 100, 250, 9, 2 * 10000, 34}, /* 31 digits */
282 {110, 100, 250, 9, 2 * 30000, 38}, /* 34 digits */
283 {120, 100, 700, 9, 2 * 40000, 49}, /* 37 digits */
284 {130, 100, 800, 9, 2 * 50000, 59}, /* 40 digits */
285 {140, 100, 1200, 9, 2 * 65536, 66}, /* 43 digits */
286 {150, 100, 1500, 10, 2 * 65536, 70}, /* 46 digit */
287 {160, 150, 2000, 11, 4 * 65536, 73}, /* 49 digit */
288 {170, 150, 2000, 12, 4 * 65536, 75}, /* 52 digits */
289 {180, 150, 3000, 12, 4 * 65536, 76}, /* 55 digits */
290 {190, 150, 3000, 13, 4 * 65536, 78}, /* 58 digit */
291 {200, 200, 4500, 14, 4 * 65536, 81}, /* 61 digits */
292 {210, 100, 8000, 14, 12 * 65536, 84}, /* 64 digits */
293 {220, 300, 10000, 15, 12 * 65536, 88}, /* 67 digits */
294 {230, 400, 20000, 17, 20 * 65536, 90}, /* 70 digits */
295 {240, 450, 20000, 19, 20 * 65536, 93}, /* 73 digis */
296 {250, 500, 22000, 22, 24 * 65536, 97}, /* 76 digits */
297 {260, 600, 25000, 25, 24 * 65536, 100}, /* 79 digits */
298 {270, 800, 35000, 27, 28 * 65536, 102}, /* 82 digits */
299 {280, 900, 40000, 29, 28 * 65536, 104}, /* 85 digits */
300 {290, 1000, 60000, 29, 32 * 65536, 106}, /* 88 digits */
301 {300, 1100, 140000, 30, 32 * 65536, 108} /* 91 digits */
302 };
303
304 #endif
305
306 /* number of entries in the tuning table */
307 #define QS_TUNE_SIZE (sizeof(qsieve_tune)/(6*sizeof(mp_limb_t)))
308
309 FLINT_DLL void qsieve_init(qs_t qs_inf, const fmpz_t n);
310
311 FLINT_DLL mp_limb_t qsieve_knuth_schroeppel(qs_t qs_inf);
312
313 FLINT_DLL void qsieve_clear(qs_t qs_inf);
314
315 FLINT_DLL void qsieve_factor(fmpz_factor_t factors, const fmpz_t n);
316
317 FLINT_DLL prime_t * compute_factor_base(mp_limb_t * small_factor, qs_t qs_inf,
318 slong num_primes);
319
320 FLINT_DLL mp_limb_t qsieve_primes_init(qs_t qs_inf);
321
322 FLINT_DLL mp_limb_t qsieve_primes_increment(qs_t qs_inf, mp_limb_t delta);
323
324 FLINT_DLL mp_limb_t qsieve_poly_init(qs_t qs_inf);
325
326 FLINT_DLL int qsieve_init_A(qs_t qs_inf);
327
328 FLINT_DLL void qsieve_reinit_A(qs_t qs_inf);
329
330 FLINT_DLL int qsieve_next_A(qs_t qs_inf);
331
332 FLINT_DLL void qsieve_init_poly_first(qs_t qs_inf);
333
334 FLINT_DLL void qsieve_init_poly_next(qs_t qs_inf, slong i);
335
336 FLINT_DLL void qsieve_compute_C(fmpz_t C, qs_t qs_inf, qs_poly_t poly);
337
338 FLINT_DLL void qsieve_poly_copy(qs_poly_t poly, qs_t qs_inf);
339
340 FLINT_DLL void qsieve_poly_clear(qs_t qs_inf);
341
342 FLINT_DLL void qsieve_do_sieving(qs_t qs_inf, unsigned char * sieve, qs_poly_t poly);
343
344 FLINT_DLL void qsieve_do_sieving2(qs_t qs_inf, unsigned char * sieve, qs_poly_t poly);
345
346 FLINT_DLL slong qsieve_evaluate_candidate(qs_t qs_inf, ulong i, unsigned char * sieve, qs_poly_t poly);
347
348 FLINT_DLL slong qsieve_evaluate_sieve(qs_t qs_inf, unsigned char * sieve, qs_poly_t poly);
349
350 FLINT_DLL slong qsieve_collect_relations(qs_t qs_inf, unsigned char * sieve);
351
352 FLINT_DLL void qsieve_linalg_init(qs_t qs_inf);
353
354 FLINT_DLL void qsieve_linalg_realloc(qs_t qs_inf);
355
356 FLINT_DLL void qsieve_linalg_clear(qs_t qs_inf);
357
358 FLINT_DLL int qsieve_relations_cmp(const void * a, const void * b);
359
360 FLINT_DLL slong qsieve_merge_relations(qs_t qs_inf);
361
362 FLINT_DLL void qsieve_write_to_file(qs_t qs_inf, mp_limb_t prime,
363 fmpz_t Y, qs_poly_t poly);
364
365 FLINT_DLL hash_t * qsieve_get_table_entry(qs_t qs_inf, mp_limb_t prime);
366
367 FLINT_DLL void qsieve_add_to_hashtable(qs_t qs_inf, mp_limb_t prime);
368
369 FLINT_DLL relation_t qsieve_parse_relation(qs_t qs_inf, char * str);
370
371 FLINT_DLL relation_t qsieve_merge_relation(qs_t qs_inf, relation_t a, relation_t b);
372
373 FLINT_DLL int qsieve_compare_relation(const void * a, const void * b);
374
375 FLINT_DLL int qsieve_remove_duplicates(relation_t * rel_list, slong num_relations);
376
377 FLINT_DLL void qsieve_insert_relation(qs_t qs_inf, relation_t * rel_list,
378 slong num_relations);
379
380 FLINT_DLL int qsieve_process_relation(qs_t qs_inf);
381
insert_col_entry(la_col_t * col,slong entry)382 static __inline__ void insert_col_entry(la_col_t * col, slong entry)
383 {
384 if (((col->weight >> 4) << 4) == col->weight) /* need more space */
385 {
386 if (col->weight != 0) col->data =
387 (slong *) flint_realloc(col->data, (col->weight + 16)*sizeof(slong));
388 else col->data = (slong *) flint_malloc(16*sizeof(slong));
389 }
390
391 col->data[col->weight] = entry;
392 col->weight++;
393 }
394
swap_cols(la_col_t * col2,la_col_t * col1)395 static __inline__ void swap_cols(la_col_t * col2, la_col_t * col1)
396 {
397 la_col_t temp;
398
399 temp.weight = col1->weight;
400 temp.data = col1->data;
401 temp.orig = col1->orig;
402
403 col1->weight = col2->weight;
404 col1->data = col2->data;
405 col1->orig = col2->orig;
406
407 col2->weight = temp.weight;
408 col2->data = temp.data;
409 col2->orig = temp.orig;
410 }
411
clear_col(la_col_t * col)412 static __inline__ void clear_col(la_col_t * col)
413 {
414 col->weight = 0;
415 }
416
free_col(la_col_t * col)417 static __inline__ void free_col(la_col_t * col)
418 {
419 if (col->weight) flint_free(col->data);
420 }
421
422 FLINT_DLL uint64_t get_null_entry(uint64_t * nullrows, slong i, slong l);
423
424 FLINT_DLL void reduce_matrix(qs_t qs_inf, slong *nrows, slong *ncols, la_col_t *cols);
425
426 FLINT_DLL uint64_t * block_lanczos(flint_rand_t state, slong nrows,
427 slong dense_rows, slong ncols, la_col_t *B);
428
429 FLINT_DLL void qsieve_square_root(fmpz_t X, fmpz_t Y, qs_t qs_inf,
430 uint64_t * nullrows, slong ncols, slong l, fmpz_t N);
431
432 #ifdef __cplusplus
433 }
434 #endif
435
436 #endif
437