include/flint/qsieve.h

0001 /*
0002     Copyright (C) 2006, 2011, 2016 William Hart
0003     Copyright (C) 2015 Nitin Kumar
0004
0005     This file is part of FLINT.
0006
0007     FLINT is free software: you can redistribute it and/or modify it under
0008     the terms of the GNU Lesser General Public License (LGPL) as published
0009     by the Free Software Foundation; either version 3 of the License, or
0010     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0011 */
0012
0013 #ifndef QSIEVE_H
0014 #define QSIEVE_H
0015
0016 #include <stdint.h>
0017 #include "fmpz_types.h"
0018
0019 #if FLINT_USES_PTHREAD
0020 # include <pthread.h>
0021 #endif
0022
0023 #ifdef __cplusplus
0024 extern "C" {
0025 #endif
0026
0027 /* Windows systems may define `small` macro, which leads to conflicts */
0028 #undef small
0029
0030 #define QS_DEBUG 0 /* level of debug information printed, 0 = none */
0031
0032 #define BITS_ADJUST 25 /* add to sieve entries to compensate approximations */
0033
0034 #define BLOCK_SIZE (4*65536) /* size of sieving cache block */
0035
0036 typedef struct
0037 {
0038    ulong pinv;     /* precomputed inverse */
0039    ulong pinv2;    /* alternative precomputed inverse */
0040    int p;              /* prime */
0041    char size;
0042 } prime_t;
0043
0044 typedef struct          /* struct for factors of relations */
0045 {
0046    slong ind;
0047    slong exp;
0048 } fac_t;
0049
0050 typedef struct           /* matrix column */
0051 {
0052    slong * data;        /* The list of occupied rows in this column */
0053    slong weight;        /* Number of nonzero entries in this column */
0054    slong orig;         /* Original relation number */
0055 } la_col_t;
0056
0057
0058 typedef struct          /* entry in hash table */
0059 {
0060    ulong prime;    /* value of prime */
0061    ulong next;     /* next prime which have same hash value as 'prime' */
0062    ulong count;    /* number of occurrence of 'prime' */
0063 } hash_t;
0064
0065 typedef struct             /* format for relation */
0066 {
0067    ulong lp;          /* large prime, is 1, if relation is full */
0068    slong num_factors;     /* number of factors, excluding small factor */
0069    slong small_primes;   /* number of small factors */
0070    slong * small;         /* exponent of small factors */
0071    fac_t * factor;        /* factor of relation */
0072    fmpz_t Y;              /* square root of sieve value for relation */
0073 } relation_t;
0074
0075 typedef struct
0076 {
0077    fmpz_t B;          /* current B coeff of poly */
0078    int * soln1;       /* first start position in sieve per prime */
0079    int * soln2;       /* second start position in sieve per prime */
0080    int * posn1;       /* temp space for sieving */
0081    int * posn2;       /* temp space for sieving */
0082    slong * small;     /* exponents of small prime factors in relations */
0083    fac_t * factor;    /* factors for a relation */
0084    slong num_factors; /* number of factors found in a relation */
0085 } qs_poly_s;
0086
0087 typedef qs_poly_s qs_poly_t[1];
0088
0089 typedef struct
0090 {
0091    volatile slong index_j;
0092 #if FLINT_USES_PTHREAD
0093    pthread_mutex_t mutex;
0094 #endif
0095    thread_pool_handle * handles;
0096    slong num_handles;
0097
0098    fmpz_t n;               /* Number to factor */
0099
0100    flint_bitcnt_t bits;    /* Number of bits of n */
0101
0102    ulong ks_primes;        /* number of Knuth-Schroeppel primes */
0103
0104    slong fb_primes;        /* number of factor base primes to use (incl. k and 2) */
0105
0106    ulong k;            /* Multiplier */
0107    fmpz_t kn;              /* kn as a multiprecision integer */
0108
0109    slong num_primes;       /* number of factor base primes including k and 2 */
0110
0111    prime_t * factor_base;  /* data about factor base primes */
0112
0113    int * sqrts;            /* square roots of kn mod factor base primes */
0114
0115    slong small_primes;     /* number of primes to not sieve with */
0116    slong second_prime;     /* index of first prime bigger than block size */
0117    slong sieve_size;       /* size of sieve to use */
0118
0119    unsigned char sieve_bits;  /* sieve threshold */
0120    unsigned char sieve_fill;  /* for biasing sieve values */
0121
0122    /***************************************************************************
0123                        POLYNOMIAL DATA
0124     **************************************************************************/
0125
0126    fmpz_t A;                /* current value of coeff A of poly Ax^2 + Bx + C */
0127    fmpz_t B;                /* B value of poly */
0128
0129    ulong * A_ind;       /* indices of factor base primes dividing A */
0130
0131    fmpz_t * A_divp;         /* A_divp[i] = (A/p_i),
0132                                where the p_i are the prime factors of A */
0133    ulong * B0_terms;    /* B0_terms[i] = min(gamma_i, p - gamma_i) where
0134                                gamma_i = (sqrt(kn)*(A_divp[i])^(-1)) mod p_i,
0135                                where the p_i are the prime factors of A */
0136
0137    fmpz_t * B_terms;        /* B_terms[i] = A_divp[i]*B0_terms[i] (multprec) */
0138
0139    ulong * A_inv;       /* A_inv[k] = A^(-1) mod p_k, for FB prime p_k */
0140    ulong ** A_inv2B;    /* A_inv2B[i][k] = 2 * B_terms[i] * A^(-1) mod p_k
0141                                for FB prime p_k */
0142
0143    int * soln1;             /* soln1[k] = first poly root mod FB prime p_k */
0144    int * soln2;             /* soln2[k] = second poly root mod FB prime p_k */
0145
0146    fmpz_t target_A;         /* approximate target value for A coeff of poly */
0147
0148    fmpz_t upp_bound;        /* upper bound on desired A = 2*target_A */
0149    fmpz_t low_bound;        /* lower bound on desired A = target_A/2 */
0150
0151    slong s;                 /* number of prime factors of A */
0152    slong low;               /* minimum offset in factor base,
0153                                for possible factors of A */
0154    slong high;              /* maximum offset in factor base,
0155                                for possible factors of A */
0156    slong span;              /* total number of possible factors for A */
0157
0158    /*
0159       parameters for calculating next subset of possible factors of A, giving
0160       a lexicographic ordering of all such tuples
0161    */
0162    slong h; /* tuple entry we just set, numbered from 1 at end of tuple */
0163    slong m; /* last value we just set a tuple entry to */
0164    slong A_ind_diff; /* diff. between indices of (s-1) and (s-2)-th A-factor */
0165    ulong * curr_subset; /* current tuple */
0166    ulong * first_subset; /* first tuple, in case of restart */
0167    ulong j; /* index of s-th factor of first A, if s > 3 */
0168
0169 #if QS_DEBUG
0170    slong poly_count;         /* keep track of the number of polynomials used */
0171 #endif
0172
0173    qs_poly_s * poly;         /* poly data per thread */
0174
0175    /***************************************************************************
0176                        RELATION DATA
0177    ***************************************************************************/
0178
0179    FLINT_FILE * siqs;           /* pointer to file for storing relations */
0180    char * fname;          /* name of file used for relations */
0181
0182    slong full_relation;   /* number of full relations */
0183    slong num_cycles;      /* number of possible full relations from partials */
0184
0185    slong vertices;        /* number of different primes in partials */
0186    slong components;      /* equal to 1 */
0187    slong edges;           /* total number of partials */
0188
0189    slong table_size;      /* size of table */
0190    hash_t * table;        /* store 'prime' occurring in partial */
0191    ulong * hash_table;  /* to keep track of location of primes in 'table' */
0192
0193    slong extra_rels;      /* number of extra relations beyond num_primes */
0194    slong max_factors;     /* maximum number of factors a relation can have */
0195
0196    fmpz * Y_arr;          /* array of Y's corresponding to relations */
0197    slong * curr_rel;      /* current relation in array of relations */
0198    slong * relation;      /* relation array */
0199
0200    slong buffer_size;     /* size of buffer of relations */
0201    slong num_relations;   /* number of relations so far */
0202
0203    ulong small_factor;    /* small factor found when merging relations */
0204
0205    /***************************************************************************
0206                        LINEAR ALGEBRA DATA
0207    ***************************************************************************/
0208
0209    la_col_t * matrix; /* the main matrix over GF(2) in sparse format */
0210    la_col_t ** qsort_arr; /* array of columns ready to be sorted */
0211
0212    slong columns; /* number of columns in matrix so far */
0213
0214    /***************************************************************************
0215                        SQUARE ROOT DATA
0216    ***************************************************************************/
0217
0218    slong * prime_count; /* counts of the exponents of primes appearing in the square */
0219
0220 } qs_s;
0221
0222 typedef qs_s qs_t[1];
0223
0224 /*
0225    Tuning parameters { bits, ks_primes, fb_primes, small_primes, sieve_size}
0226    for qsieve_factor_threaded where:
0227      * bits is the number of bits of n
0228      * ks_primes is the max number of primes to try in Knuth-Schroeppel function
0229      * fb_primes is the number of factor base primes to use (including k and 2)
0230      * small_primes is the number of small primes to not factor with (including k and 2)
0231      * sieve_size is the size of the sieve to use
0232      * sieve_bits - sieve_fill
0233 */
0234
0235 static const ulong qsieve_tune[][6] =
0236 {
0237     /* Generated with qsieve-tune */
0238     {10,  49,     84,  4,    4300,  25},  /* 3 digits, 0.00 ms/semiprime */
0239     {20,  10,    120,  5,    4800,  26},  /* 6 digits, 0.01 ms/semiprime */
0240     {30,  14,     25,  3,    4500,  19},  /* 9 digits, 0.72 ms/semiprime */
0241     {40,  28,     73,  9,    5700,  24},  /* 12 digits, 1.05 ms/semiprime */
0242     {50,  25,     85,  7,    5400,  28},  /* 15 digits, 1.25 ms/semiprime */
0243     {60,  27,    100, 14,   11000,  28},  /* 18 digits, 1.45 ms/semiprime */
0244     {65,  26,    110, 19,   13000,  28},  /* 19 digits, 1.86 ms/semiprime */
0245     {70,  26,    180,  3,   12000,  43},  /* 21 digits, 2.26 ms/semiprime */
0246     {75,  26,    190,  4,   12000,  45},  /* 22 digits, 2.61 ms/semiprime */
0247     {80,  27,    190,  5,   12000,  45},  /* 24 digits, 3.16 ms/semiprime */
0248     {85,  28,    150,  5,   10000,  45},  /* 25 digits, 2.34 ms/semiprime */
0249     {90,  34,    150,  5,   10000,  45},  /* 27 digits, 2.68 ms/semiprime */
0250     {100, 40,    210,  6,   14000,  48},  /* 30 digits, 4.20 ms/semiprime */
0251     {110, 46,    220,  7,   13000,  48},  /* 33 digits, 6.50 ms/semiprime */
0252     {120, 50,    250,  7,   13000,  50},  /* 36 digits, 11.58 ms/semiprime */
0253     {130, 47,    310,  9,   16000,  50},  /* 39 digits, 19.37 ms/semiprime */
0254     {135, 39,    450, 10,   20000,  53},  /* 40 digits, 25.39 ms/semiprime */
0255     {140, 47,    580,  8,   17000,  56},  /* 41 digits, 33.47 ms/semiprime */
0256     {145, 47,    610,  9,   16000,  56},  /* 43 digits, 46.47 ms/semiprime */
0257     {150, 45,    610,  8,   17000,  59},  /* 45 digits, 59.05 ms/semiprime */
0258     {155, 43,    860, 10,   21000,  59},  /* 46 digits, 92.20 ms/semiprime */
0259     /* Legacy tuning values that seem reasonable */
0260     {160, 150,  2000, 11,   4 *  65536, 73}, /* 49 digit */
0261     {170, 150,  2000, 12,   4 *  65536, 75}, /* 52 digits */
0262     {180, 150,  3000, 12,   4 *  65536, 76}, /* 55 digits */
0263     {190, 150,  3000, 13,   4 *  65536, 78}, /* 58 digit */
0264     {200, 200,  4500, 14,   4 *  65536, 81}, /* 61 digits */
0265     {210, 100,  8000, 14,   12 *  65536, 84}, /* 64 digits */
0266     {220, 300, 10000, 15,   12 *  65536, 88}, /* 67 digits */
0267     {230, 400, 20000, 17,   20 *  65536, 90}, /* 70 digits */
0268     {240, 450, 20000, 19,   20 *  65536, 93}, /* 73 digis */
0269     {250, 500, 22000, 22,   24 *  65536, 97}, /* 76 digits */
0270     {260, 600, 25000, 25,   24 *  65536, 100}, /* 79 digits */
0271     {270, 800, 35000, 27,   28 *  65536, 102}, /* 82 digits */
0272     {280, 900, 40000, 29,   28 *  65536, 104}, /* 85 digits */
0273     {290, 1000, 60000, 29,  32 *  65536, 106}, /* 88 digits */
0274     {300, 1100, 140000, 30,  32 * 65536, 108} /* 91 digits */
0275 };
0276
0277 /* number of entries in the tuning table */
0278 #define QS_TUNE_SIZE (sizeof(qsieve_tune)/(6*sizeof(ulong)))
0279
0280 void qsieve_init(qs_t qs_inf, const fmpz_t n);
0281
0282 /* As qsieve_init, but with the tuning parameters supplied explicitly rather
0283    than read from the qsieve_tune table.  qsieve_init calls this with the
0284    table defaults for the bit-size of n. */
0285 void qsieve_init_with_tune(qs_t qs_inf, const fmpz_t n, ulong ks_primes,
0286         slong fb_primes, slong small_primes, slong sieve_size, ulong sieve_bits);
0287
0288 ulong qsieve_knuth_schroeppel(qs_t qs_inf);
0289
0290 void qsieve_clear(qs_t qs_inf);
0291
0292 void qsieve_factor(fmpz_factor_t factors, const fmpz_t n);
0293
0294 /* As qsieve_factor, but with the tuning parameters supplied explicitly. */
0295 void qsieve_factor_with_tune(fmpz_factor_t factors, const fmpz_t n,
0296         ulong ks_primes, slong fb_primes, slong small_primes,
0297         slong sieve_size, ulong sieve_bits);
0298
0299 prime_t * compute_factor_base(ulong * small_factor, qs_t qs_inf,
0300                                                              slong num_primes);
0301
0302 ulong qsieve_primes_init(qs_t qs_inf);
0303
0304 ulong qsieve_primes_increment(qs_t qs_inf, ulong delta);
0305
0306 ulong qsieve_poly_init(qs_t qs_inf);
0307
0308 int qsieve_init_A(qs_t qs_inf);
0309
0310 void qsieve_reinit_A(qs_t qs_inf);
0311
0312 int qsieve_next_A(qs_t qs_inf);
0313
0314 void qsieve_init_poly_first(qs_t qs_inf);
0315
0316 void qsieve_init_poly_next(qs_t qs_inf, slong i);
0317
0318 void qsieve_compute_C(fmpz_t C, qs_t qs_inf, qs_poly_t poly);
0319
0320 void qsieve_poly_copy(qs_poly_t poly, qs_t qs_inf);
0321
0322 void qsieve_poly_clear(qs_t qs_inf);
0323
0324 void qsieve_do_sieving(qs_t qs_inf, unsigned char * sieve, qs_poly_t poly);
0325
0326 void qsieve_do_sieving2(qs_t qs_inf, unsigned char * sieve, qs_poly_t poly);
0327
0328 slong qsieve_evaluate_candidate(qs_t qs_inf, ulong i, unsigned char * sieve, qs_poly_t poly);
0329
0330 slong qsieve_evaluate_sieve(qs_t qs_inf, unsigned char * sieve, qs_poly_t poly);
0331
0332 slong qsieve_collect_relations(qs_t qs_inf, unsigned char * sieve);
0333
0334 void qsieve_linalg_init(qs_t qs_inf);
0335
0336 void qsieve_linalg_realloc(qs_t qs_inf);
0337
0338 void qsieve_linalg_clear(qs_t qs_inf);
0339
0340 int qsieve_relations_cmp(const void * a, const void * b);
0341
0342 slong qsieve_merge_relations(qs_t qs_inf);
0343
0344 void qsieve_write_to_file(qs_t qs_inf, ulong prime,
0345                                                      const fmpz_t Y, const qs_poly_t poly);
0346
0347 hash_t * qsieve_get_table_entry(qs_t qs_inf, ulong prime);
0348
0349 void qsieve_add_to_hashtable(qs_t qs_inf, ulong prime);
0350
0351 relation_t qsieve_parse_relation(qs_t qs_inf);
0352
0353 relation_t qsieve_merge_relation(qs_t qs_inf, relation_t  a, relation_t  b);
0354
0355 int qsieve_compare_relation(const void * a, const void * b);
0356
0357 int qsieve_remove_duplicates(relation_t * rel_list, slong num_relations);
0358
0359 void qsieve_insert_relation(qs_t qs_inf, relation_t * rel_list,
0360                                                           slong num_relations);
0361
0362 int qsieve_process_relation(qs_t qs_inf);
0363
0364 static inline void insert_col_entry(la_col_t * col, slong entry)
0365 {
0366    if (((col->weight >> 4) << 4) == col->weight) /* need more space */
0367    {
0368        if (col->weight != 0) col->data =
0369            (slong *) flint_realloc(col->data, (col->weight + 16)*sizeof(slong));
0370        else col->data = (slong *) flint_malloc(16*sizeof(slong));
0371    }
0372
0373    col->data[col->weight] = entry;
0374    col->weight++;
0375 }
0376
0377 static inline void swap_cols(la_col_t * col2, la_col_t * col1)
0378 {
0379    la_col_t temp;
0380
0381    temp.weight = col1->weight;
0382    temp.data = col1->data;
0383    temp.orig = col1->orig;
0384
0385    col1->weight = col2->weight;
0386    col1->data = col2->data;
0387    col1->orig = col2->orig;
0388
0389    col2->weight = temp.weight;
0390    col2->data = temp.data;
0391    col2->orig = temp.orig;
0392 }
0393
0394 static inline void clear_col(la_col_t * col)
0395 {
0396    col->weight = 0;
0397 }
0398
0399 static inline void free_col(la_col_t * col)
0400 {
0401    if (col->weight) flint_free(col->data);
0402 }
0403
0404 uint64_t get_null_entry(uint64_t * nullrows, slong i, slong l);
0405
0406 void reduce_matrix(qs_t qs_inf, slong *nrows, slong *ncols, la_col_t *cols);
0407
0408 uint64_t * block_lanczos(flint_rand_t state, slong nrows,
0409             slong dense_rows, slong ncols, la_col_t *B);
0410
0411 void qsieve_square_root(fmpz_t X, fmpz_t Y, qs_t qs_inf,
0412    uint64_t * nullrows, slong ncols, slong l, fmpz_t N);
0413
0414 #ifdef __cplusplus
0415 }
0416 #endif
0417
0418 #endif