/*-------------------------------------------------------------------- This source distribution is placed in the public domain by its author, Jason Papadopoulos. You may use it for any purpose, free of charge, without having to notify anyone. I disclaim any responsibility for any errors. Optionally, please be nice and tell me if you find this source to be useful. Again optionally, if you add to the functionality present here please consider making those additions public too, so that others may benefit from your work. --jasonp@boo.net 6/3/07 --------------------------------------------------------------------*/ #ifndef _MPQS_H_ #define _MPQS_H_ /* An implementation of the Self-Initializing Multiple Polynomial Quadratic Sieve algorithm for integer factorization. */ #include #ifdef __cplusplus extern "C" { #endif /*------------------- SIEVE RELATED DECLARATIONS ---------------------*/ /* There's a limit to how many factors can contribute to SIQS polynomials */ #define MAX_POLY_FACTORS 20 /* The size of a factor base prime is limited too */ #define MAX_FB_PRIME 0x3ffffff /* representing one prime in the factor base */ #define INVALID_ROOT ((uint32)(-1)) typedef struct { uint32 prime : 26; /* the factor base prime */ uint32 logprime : 5; /* log2(prime)rounded to nearest */ uint32 rcorrect : 1; /* used to correct multiples by recip */ uint32 root1; /* the two square roots of an MPQS */ uint32 root2; /* polynomial (INVALID_ROOT means don't use) */ uint32 recip; /* integer reciprocal of prime */ } fb_t; /* A packed representation of the above structure, used during the sieving phase only. Careful: the room allowed here for the factor base prime is less than is allowed in the ordinary factor base. That's okay, since this structure is only used for factor base primes less than the sieve block size */ typedef struct { uint32 prime : 24; /* the factor base prime */ uint8 logprime; /* log2(prime)rounded to nearest */ uint32 next_loc1; /* the two next sieve locations for */ uint32 next_loc2; /* the above two roots */ } packed_fb_t; /* Offset 0 in the array of factor base primes is reserved for the "prime" -1. This allows the sieve to be over both positive and negative values. When actually dividing by factor base primes, start with the fb_t at this position in the factor base array: */ #define MIN_FB_OFFSET 1 /* structure representing one MPQS polynomial. Every relation found by the sieve must be associated with (i.e. must refer to) one of these. For sieve offset x, the MPQS polynomial is a * x + b */ typedef struct poly_t { uint32 a_idx; /* offset into a list of 'a' values */ mp_t b; /* the MPQS 'b' value */ } poly_t; /* If a sieve value turns out to be smooth, the following is used to record information for use in later stages of the algorithm. One sieve value is packed into a relation_t, and an la_col_t contains a collection of sieve values */ #define POSITIVE 0 #define NEGATIVE 1 typedef struct relation_t { uint32 sieve_offset; /* value of x in a*x+b */ uint32 poly_idx; /* pointer to this relation's MPQS poly */ uint32 num_factors; /* size of list of factors */ uint32 large_prime[2]; /* for partial relations, the leftover factor(s). Either factor may equal 1 */ uint32 *fb_offsets; /* The array of offsets into the factor base which contain primes dividing this sieve value. Duplicate offsets are possible. Sorted in ascending order */ } relation_t; /* Once the sieving stage has collected enough relations, we'll need a few extra relations to make sure that the linear algebra stage finds linear dependencies */ #define NUM_EXTRA_RELATIONS 64 /* The sieving phase uses a hashtable to handle the large factor base primes. The hashtable is an array of type hashtable_t, each entry of which contains an array of type hash_entry_t */ typedef struct { uint32 logprime : 8; /* the logprime field from the factor base */ uint32 prime_index : 24; /* factor base array offset for entry */ uint32 sieve_offset; /* offset within one hash bin */ } hash_entry_t; typedef struct { uint32 num_alloc; /* total list size this hashtable position */ uint32 num_used; /* number occupied at this hashtable position */ hash_entry_t *list; /* list of entries at this hashtable position */ } hashtable_t; /* Configuration of the sieving code requires passing the following structure */ typedef struct { uint32 bits; /* size of integer this config info applies to */ uint32 fb_size; /* number of factor base primes */ uint32 large_mult; /* the large prime multiplier */ uint32 sieve_size; /* the size of the sieve (actual sieve is 2x this) */ } sieve_param_t; /* The sieving code needs a special structure for determining the number of cycles in a collection of partial relations */ #define LOG2_CYCLE_HASH 22 typedef struct { uint32 next; uint32 prime; uint32 data; uint32 count; } cycle_t; /* To avoid huge parameter lists passed between sieving routines, all of the relevant data used in the sieving phase is packed into a single structure. Routines take the information they need out of this. */ typedef struct { msieve_obj *obj; /* object controlling entire factorization */ mp_t *n; /* the number to factor (scaled by multiplier)*/ uint32 multiplier; /* small multiplier for n (may be composite) */ /* bookkeeping information for sieving */ uint32 sieve_block_size; /* bytes in one sieve block */ uint8 *sieve_array; /* scratch space used for one sieve block */ fb_t *factor_base; /* the factor base to use */ packed_fb_t *packed_fb; /* scratch space for a packed version of it */ uint32 fb_size; /* number of factor base primes (includes -1)*/ uint32 num_sieve_blocks; /* number of sieve blocks in sieving interval*/ uint32 poly_block; /* number of FB primes in a block of work */ uint32 fb_block; /* number of polynomials simultaneously sieved */ uint32 sieve_small_fb_start; /* starting FB offset for sieving */ uint32 sieve_large_fb_start; uint32 tf_small_recip1_cutoff; uint32 tf_small_recip2_cutoff; uint32 tf_med_recip1_cutoff; uint32 tf_med_recip2_cutoff; uint32 tf_large_cutoff; hashtable_t *hashtable; /* hash bins for sieve values */ uint32 cutoff1; /* if log2(sieve value) exceeds this number, a little trial division is performed */ uint32 cutoff2; /* if log2(sieve value) exceeds this number, full trial division is performed */ /* bookkeeping information for creating polynomials */ uint32 *modsqrt_array; /* a square root of n mod each FB prime */ mp_t target_a; /* optimal value of 'a' */ uint32 a_bits; /* required number of bits in the 'a' value of all MPQS polynomials */ uint32 total_poly_a; /* total number of polynomial 'a' values */ mp_t *poly_a_list; /* list of 'a' values for MPQS polys */ poly_t *poly_list; /* list of MPQS polynomials */ mp_t curr_a; /* the current 'a' value */ mp_t *curr_b; /* list of all the 'b' values for that 'a' */ uint32 *root_correction; /* see poly.c */ uint8 *next_poly_action; /* see poly.c */ uint32 num_poly_factors; /* number of factors in poly 'a' value */ uint32 factor_bounds[25]; /* bounds on FB offsets for primes that */ /* can be factors of MPQS polys */ uint32 poly_factors[MAX_POLY_FACTORS]; /* factorization of curr. 'a' */ uint8 factor_bits[MAX_POLY_FACTORS]; /* size of each factor of 'a' */ mp_t poly_tmp_b[MAX_POLY_FACTORS]; /* temporary quantities */ uint32 *poly_b_array; /* precomputed values for all factor base primes, used to compute new polynomials */ uint32 *poly_b_small[MAX_POLY_FACTORS]; /* bookkeeping information for double large primes */ uint32 large_prime_max; /* the cutoff value for keeping a partial relation; actual value, not a multiplier */ mp_t max_fb2; /* the square of the largest factor base prime */ mp_t large_prime_max2; /* the cutoff value for factoring partials */ relation_t *relation_list; /* list of full/partial relations */ uint32 num_relations; /* number of relations in list */ la_col_t *cycle_list; /* cycles derived from relations */ uint32 num_cycles; /* number of cycles in list */ cycle_t *cycle_table; /* list of all the vertices in the graph */ uint32 cycle_table_size; /* number of vertices filled in the table */ uint32 cycle_table_alloc; /* number of cycle_t structures allocated */ uint32 *cycle_hashtable; /* hashtable to index into cycle_table */ uint32 components; /* connected components (see relation.c) */ uint32 vertices; /* vertices in graph (see relation.c) */ } sieve_conf_t; /* attempt to trial factor one sieve value */ uint32 check_sieve_val(sieve_conf_t *conf, int32 sieve_offset, uint32 bits, mp_t *a, mp_t *b, mp_t *c, uint32 poly_index, hashtable_t *hash_bucket); /* the core of the sieving code */ typedef uint32 (*qs_core_sieve_fcn)(sieve_conf_t *conf, uint32 poly_start, uint32 num_poly); #define DECLARE_SIEVE_FCN(name) \ uint32 name(sieve_conf_t *conf, \ uint32 poly_start, \ uint32 num_poly); DECLARE_SIEVE_FCN(qs_core_sieve_generic_32k); DECLARE_SIEVE_FCN(qs_core_sieve_generic_64k); #if defined(__GNUC__) #if defined(__i386__) DECLARE_SIEVE_FCN(qs_core_sieve_p2_64k); DECLARE_SIEVE_FCN(qs_core_sieve_p3_64k); DECLARE_SIEVE_FCN(qs_core_sieve_p4_64k); DECLARE_SIEVE_FCN(qs_core_sieve_pm_32k); DECLARE_SIEVE_FCN(qs_core_sieve_core_32k); DECLARE_SIEVE_FCN(qs_core_sieve_k7_64k); DECLARE_SIEVE_FCN(qs_core_sieve_k7xp_64k); DECLARE_SIEVE_FCN(qs_core_sieve_k8_64k); #elif defined(__x86_64__) DECLARE_SIEVE_FCN(qs_core_sieve_p4_64k); DECLARE_SIEVE_FCN(qs_core_sieve_core_32k); DECLARE_SIEVE_FCN(qs_core_sieve_k8_64k); #endif #elif defined(_MSC_VER) && (_MSC_VER >= 1400) DECLARE_SIEVE_FCN(qs_core_sieve_vc8_32k); DECLARE_SIEVE_FCN(qs_core_sieve_vc8_64k); #endif /* lightweight profiling of the sieve routines */ #ifdef SIEVE_TIMING #define TIME1(var) { uint64 tmp_##var = read_clock(); #define TIME2(var) (var) += read_clock() - tmp_##var; } extern uint64 next_poly_small_time; extern uint64 next_poly_large_time; extern uint64 plus_init_time; extern uint64 bucket_time; extern uint64 sieve_small_time; extern uint64 sieve_large_time; extern uint64 tf_plus_scan_time; extern uint64 tf_total_time; #else #define TIME1(var) /* nothing */ #define TIME2(var) /* nothing */ #endif /* pull out the large primes from a relation read from the savefile */ void read_large_primes(char *buf, uint32 *prime1, uint32 *prime2); /* given the primes from a sieve relation, add that relation to the graph used for tracking cycles */ void add_to_cycles(sieve_conf_t *conf, uint32 large_prime1, uint32 large_prime2); /* encapsulate all of the information conerning a sieve relation and dump it to the savefile */ void save_relation(sieve_conf_t *conf, uint32 sieve_offset, uint32 *fb_offsets, uint32 num_factors, uint32 poly_index, uint32 large_prime1, uint32 large_prime2); void qs_print_to_savefile(msieve_obj *obj, char *buf); void qs_flush_savefile(msieve_obj *obj); /* perform postprocessing on a list of relations */ void qs_filter_relations(sieve_conf_t *conf); /* Initialize/free polynomial scratch data. Note that poly_free can only be called after the relation filtering phase */ void poly_init(sieve_conf_t *conf, uint32 sieve_size); void poly_free(sieve_conf_t *conf); /* compute a random polynomial 'a' value, and also compute all of the 'b' values, all of the precomputed quantities for the 'b' values, and all of the initial roots in the factor base */ void build_base_poly(sieve_conf_t *conf); /* As above, except it's assumed that conf->poly_factors and conf->curr_a are already filled in. Factor base roots are not affected */ void build_derived_poly(sieve_conf_t *conf); /* The main function to perform sieving. obj is the object controlling this factorization n is the number to be factored poly_list is a linked list of all the polynomials created factor_base contains the factor base (duh) modsqrt_array contains modular square roots of FB primes params is used to initialize and configure the sieve code multiplier is a small integer by which n is scaled relation_list is the list of relations from the sieving stage num_relations is the size of relation_list cycle_list is the list of cycles that the QS filtering code builds num_cycles is the number of cycles */ void do_sieving(msieve_obj *obj, mp_t *n, mp_t **poly_a_list, poly_t **poly_list, fb_t *factor_base, uint32 *modsqrt_array, sieve_param_t *params, uint32 multiplier, relation_t **relation_list, uint32 *num_relations, la_col_t **cycle_list, uint32 *num_cycles); void qs_free_relation_list(relation_t *list, uint32 num_relations); /*--------------LINEAR ALGEBRA RELATED DECLARATIONS ---------------------*/ /* Find linear dependencies. The number of nontrivial dependencies found is returned obj is the object controlling this factorization fb_size is the size of the factor base bitfield is an array of num_cycles numbers. Bit i of word j tells whether cycle_list[j] is used in nullspace vector i. Essentially, bitfield[] is a collection of 64 nullspace vectors packed together. relation_list is the list of relations from the sieving stage num_relations is the size of relation_list cycle_list is the list of cycles that the linear algebra code constructs (on input it is the list of cycles that the QS filtering code has found) num_cycles is the number of cycles. On input this is the size of cycle_list. The linear algebra code can change this number; the only guarantee is that its final value is at least fb_size + NUM_EXTRA_RELATIONS */ void solve_linear_system(msieve_obj *obj, uint32 fb_size, uint64 **bitfield, relation_t *relation_list, la_col_t *cycle_list, uint32 *num_cycles); /*-------------- MPQS SQUARE ROOT RELATED DECLARATIONS ---------------------*/ uint32 find_factors(msieve_obj *obj, mp_t *n, fb_t *factor_base, uint32 fb_size, la_col_t *vectors, uint32 vsize, relation_t *relation_list, uint64 *null_vectors, uint32 multiplier, mp_t *a_list, poly_t *poly_list, factor_list_t *factor_list); #ifdef __cplusplus } #endif #endif /* _MPQS_H_ */