1 /*
2 Copyright (C) 2010 Fredrik Johansson
3 Copyright (C) 2020 William Hart
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 #include <gmp.h>
14 #include "flint.h"
15 #include "fmpz.h"
16 #include "fmpz_vec.h"
17 #include "mpn_extras.h"
18 #include "ulong_extras.h"
19
20 int
fmpz_factor_trial(fmpz_factor_t factor,const fmpz_t n,slong num_primes)21 fmpz_factor_trial(fmpz_factor_t factor, const fmpz_t n, slong num_primes)
22 {
23 ulong exp;
24 mp_limb_t p;
25 mpz_t x;
26 mp_ptr xd;
27 mp_size_t xsize;
28 slong found;
29 int ret = 1;
30 slong * idx;
31 slong bits, i;
32 const mp_limb_t * primes;
33
34 if (num_primes > 3512 || num_primes < 0)
35 {
36 flint_printf("(fmpz_factor_trial) Number of primes must be in 0..3512\n");
37 flint_abort();
38 }
39
40 if (!COEFF_IS_MPZ(*n))
41 {
42 fmpz_factor_si(factor, *n);
43
44 return ret;
45 }
46
47 _fmpz_factor_set_length(factor, 0);
48
49 /* Make an mpz_t copy whose limbs will be mutated */
50 mpz_init(x);
51 fmpz_get_mpz(x, n);
52 if (x->_mp_size < 0)
53 {
54 x->_mp_size = -(x->_mp_size);
55 factor->sign = -1;
56 }
57 else
58 {
59 factor->sign = 1;
60 }
61
62 xd = x->_mp_d;
63 xsize = x->_mp_size;
64
65 /* Factor out powers of two */
66 xsize = flint_mpn_remove_2exp(xd, xsize, &exp);
67 if (exp != 0)
68 _fmpz_factor_append_ui(factor, UWORD(2), exp);
69
70 bits = fmpz_sizeinbase(n, 2) - exp;
71 idx = (slong *) flint_malloc((5 + bits/4)*sizeof(slong));
72
73 found = flint_mpn_factor_trial_tree(idx, xd, xsize, num_primes);
74
75 if (found)
76 {
77 primes = n_primes_arr_readonly(3512);
78
79 for (i = 0; i < found; i++)
80 {
81 p = primes[idx[i]];
82
83 exp = 1;
84 xsize = flint_mpn_divexact_1(xd, xsize, p);
85
86 /* Check if p^2 divides n */
87 if (flint_mpn_divisible_1_p(xd, xsize, p))
88 {
89 /* TODO: when searching for squarefree numbers
90 (Moebius function, etc), we can abort here. */
91 xsize = flint_mpn_divexact_1(xd, xsize, p);
92 exp = 2;
93 }
94
95 /* If we're up to cubes, then maybe there are higher powers */
96 if (exp == 2 && flint_mpn_divisible_1_p(xd, xsize, p))
97 {
98 xsize = flint_mpn_divexact_1(xd, xsize, p);
99 xsize = flint_mpn_remove_power_ascending(xd, xsize, &p, 1, &exp);
100 exp += 3;
101 }
102
103 _fmpz_factor_append_ui(factor, p, exp);
104 }
105 }
106
107 /* Any factor left? */
108 if (xsize > 1 || xd[0] != 1)
109 {
110 fmpz_t cofactor;
111 mpz_t mockx; /* do not free */
112
113 fmpz_init(cofactor);
114 mockx->_mp_d = xd;
115 mockx->_mp_size = xsize;
116 mockx->_mp_alloc = x->_mp_alloc;
117
118 fmpz_set_mpz(cofactor, mockx);
119 _fmpz_factor_append(factor, cofactor, 1);
120
121 fmpz_clear(cofactor);
122
123 ret = 0;
124 }
125
126 mpz_clear(x);
127
128 flint_free(idx);
129
130 return ret;
131 }
132