1 /*
2 Copyright (C) 2010 Fredrik Johansson
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <https://www.gnu.org/licenses/>.
10 */
11
12 #include <gmp.h>
13 #include "flint.h"
14 #include "fmpz.h"
15 #include "fmpz_vec.h"
16 #include "mpn_extras.h"
17 #include "ulong_extras.h"
18
19 int
fmpz_factor_trial_range(fmpz_factor_t factor,const fmpz_t n,ulong start,ulong num_primes)20 fmpz_factor_trial_range(fmpz_factor_t factor, const fmpz_t n, ulong start, ulong num_primes)
21 {
22 ulong exp;
23 mp_limb_t p;
24 mpz_t x;
25 mp_ptr xd;
26 mp_size_t xsize;
27 slong found;
28 slong trial_start, trial_stop;
29 int ret = 1;
30
31 if (!COEFF_IS_MPZ(*n))
32 {
33 fmpz_factor_si(factor, *n);
34
35 return ret;
36 }
37
38 _fmpz_factor_set_length(factor, 0);
39
40 /* Make an mpz_t copy whose limbs will be mutated */
41 mpz_init(x);
42 fmpz_get_mpz(x, n);
43 if (x->_mp_size < 0)
44 {
45 x->_mp_size = -(x->_mp_size);
46 factor->sign = -1;
47 }
48 else
49 {
50 factor->sign = 1;
51 }
52
53 xd = x->_mp_d;
54 xsize = x->_mp_size;
55
56 /* Factor out powers of two */
57 if (start == 0)
58 {
59 xsize = flint_mpn_remove_2exp(xd, xsize, &exp);
60 if (exp != 0)
61 _fmpz_factor_append_ui(factor, UWORD(2), exp);
62 }
63
64 trial_start = FLINT_MAX(1, start);
65 trial_stop = FLINT_MIN(start + 1000, start + num_primes);
66
67 do
68 {
69 found = flint_mpn_factor_trial(xd, xsize, trial_start, trial_stop);
70
71 if (found)
72 {
73 p = n_primes_arr_readonly(found+1)[found];
74 exp = 1;
75 xsize = flint_mpn_divexact_1(xd, xsize, p);
76
77 /* Check if p^2 divides n */
78 if (flint_mpn_divisible_1_p(xd, xsize, p))
79 {
80 /* TODO: when searching for squarefree numbers
81 (Moebius function, etc), we can abort here. */
82 xsize = flint_mpn_divexact_1(xd, xsize, p);
83 exp = 2;
84 }
85
86 /* If we're up to cubes, then maybe there are higher powers */
87 if (exp == 2 && flint_mpn_divisible_1_p(xd, xsize, p))
88 {
89 xsize = flint_mpn_divexact_1(xd, xsize, p);
90 xsize = flint_mpn_remove_power_ascending(xd, xsize, &p, 1, &exp);
91 exp += 3;
92 }
93
94 _fmpz_factor_append_ui(factor, p, exp);
95 /* flint_printf("added %wu %wu\n", p, exp); */
96
97 /* Continue using only trial division whilst it is successful.
98 This allows quickly factoring huge highly composite numbers
99 such as factorials, which can arise in some applications. */
100 trial_start = found + 1;
101 trial_stop = FLINT_MIN(trial_start + 1000, start + num_primes);
102 continue;
103 }
104 else
105 {
106 /* Insert primality test, perfect power test, other factoring
107 algorithms here... */
108 trial_start = trial_stop;
109 trial_stop = FLINT_MIN(trial_start + 1000, start + num_primes);
110 }
111 } while ((xsize > 1 || xd[0] != 1) && trial_start != trial_stop);
112
113 /* Any factor left? */
114 if (xsize > 1 || xd[0] != 1)
115 ret = 0;
116
117 mpz_clear(x);
118 return ret;
119 }
120