1 /*=============================================================================
2
3 This file is part of Antic.
4
5 Antic is free software: you can redistribute it and/or modify it under
6 the terms of the GNU Lesser General Public License (LGPL) as published
7 by the Free Software Foundation; either version 2.1 of the License, or
8 (at your option) any later version. See <http://www.gnu.org/licenses/>.
9
10 =============================================================================*/
11 /******************************************************************************
12
13 Copyright 2012 William Hart
14
15 ******************************************************************************/
16
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include "profiler.h"
20 #include "flint.h"
21 #include "fmpz.h"
22 #include "qfb.h"
23
main(void)24 int main(void)
25 {
26 fmpz_t g, n, n0, p, L;
27 slong iters, i, j, depth, jmax = 10;
28 qfb_t pow, twopow;
29 ulong pr, nmodpr, mult, n0mod4;
30 qfb_hash_t * qhash;
31 int done;
32
33 fmpz_init(g);
34 fmpz_init(n);
35 fmpz_init(n0);
36 fmpz_init(p);
37 fmpz_init(L);
38
39 qfb_init(pow);
40 qfb_init(twopow);
41
42 printf("Enter number to be factored: "); fflush(stdout);
43 if (!fmpz_read(n0))
44 {
45 printf("Read failed\n");
46 abort();
47 }
48 fmpz_neg(n0, n0);
49
50 iters = 10;
51
52 /* find prime such that n is a square mod p (or p divides n) */
53 if (fmpz_is_even(n0))
54 {
55 printf("Factor: 2\n");
56 return 0;
57 }
58
59 mult = 1;
60
61 while (1) /* keep increasing iterations for each multiplier until done */
62 {
63 printf("iters = %ld, multipliers = %ld\n", iters, jmax);
64
65 for (j = 0; j < jmax; j++) /* loop over jmax different multipliers */
66 {
67
68 done = 1;
69
70 if (mult != 1 && fmpz_fdiv_ui(n0, mult) == 0)
71 {
72 printf("Factor: %ld\n", mult);
73 return 0;
74 }
75
76 fmpz_mul_ui(n, n0, mult);
77 if (fmpz_fdiv_ui(n, 4) == 3)
78 fmpz_mul_2exp(n, n, 2);
79
80 pr = 2;
81 fmpz_abs(L, n);
82 fmpz_root(L, L, 4);
83
84 do
85 {
86 pr = n_nextprime(pr, 0);
87 while (mult % pr == 0)
88 pr = n_nextprime(pr, 0);
89
90 nmodpr = fmpz_fdiv_ui(n, pr);
91
92 if (nmodpr == 0) /* pr is a factor */
93 {
94 printf("Factor: %lu\n", pr);
95 return 0;
96 }
97 } while (n_jacobi(nmodpr, pr) < 0);
98
99 fmpz_set_ui(p, pr);
100
101 /* find prime form of discriminant n */
102 qfb_prime_form(pow, n, p);
103
104 /* raise to various powers of small primes */
105 qfb_pow_ui(pow, pow, n, 59049); /* 3^10 */
106 qfb_pow_ui(twopow, pow, n, 4096);
107 if (qfb_is_principal_form(twopow, n))
108 goto done;
109
110 qfb_pow_ui(pow, pow, n, 390625); /* 5^8 */
111 qfb_pow_ui(twopow, pow, n, 4096);
112 if (qfb_is_principal_form(twopow, n))
113 goto done;
114
115 qfb_pow_ui(pow, pow, n, 117649); /* 7^6 */
116 qfb_pow_ui(twopow, pow, n, 4096);
117 if (qfb_is_principal_form(twopow, n))
118 goto done;
119
120 qfb_pow_ui(pow, pow, n, 14641); /* 11^4 */
121 qfb_pow_ui(twopow, pow, n, 4096);
122 if (qfb_is_principal_form(twopow, n))
123 goto done;
124
125 qfb_pow_ui(pow, pow, n, 169); /* 13^2 */
126 qfb_pow_ui(twopow, pow, n, 4096);
127 if (qfb_is_principal_form(twopow, n))
128 goto done;
129
130 depth = FLINT_BIT_COUNT(iters) + 1;
131 qhash = qfb_hash_init(depth);
132
133 pr = 13;
134 for (i = 0; i < iters; i++)
135 {
136 pr = n_nextprime(pr, 0);
137 qfb_pow_ui(pow, pow, n, pr*pr);
138 qfb_pow_ui(twopow, pow, n, 4096);
139 if (qfb_is_principal_form(twopow, n)) /* found factor */
140 break;
141 qfb_hash_insert(qhash, twopow, pow, i, depth);
142 }
143
144 if (i == iters) /* stage 2 */
145 {
146 ulong jump = iters*iters;
147 slong iters2;
148
149 for (i = 0; i < iters; i++)
150 {
151 qfb_pow_ui(pow, pow, n, jump);
152 qfb_pow_ui(twopow, pow, n, 4096);
153 if (qfb_is_principal_form(twopow, n)) /* found factor */
154 break;
155 iters2 = qfb_hash_find(qhash, twopow, depth);
156 if (iters2 != -1) /* found factor */
157 {
158 if (fmpz_sgn(qhash[iters2].q->b) == fmpz_sgn(twopow->b))
159 qfb_inverse(qhash[iters2].q2, qhash[iters2].q2);
160
161 qfb_nucomp(pow, pow, qhash[iters2].q2, n, L);
162 qfb_reduce(pow, pow, n);
163
164 break;
165 }
166 }
167
168 if (i == iters)
169 done = 0;
170 }
171
172 done:
173 if (done)
174 {
175 for (i = 0; i < 12; i++)
176 {
177 qfb_pow_ui(twopow, pow, n, 2);
178 if (qfb_is_principal_form(twopow, n))
179 {
180 if (fmpz_cmpabs(pow->a, pow->b) != 0)
181 {
182 fmpz_abs(pow->b, pow->b);
183 fmpz_sub(g, pow->b, pow->a);
184 fmpz_sub(pow->a, g, pow->a);
185 }
186
187 fmpz_gcd(g, pow->a, n0);
188
189 if (!fmpz_is_one(g)) /* Success! */
190 {
191 printf("Factor: ");
192 fmpz_print(g);
193 printf("\n");
194 return 0;
195 }
196
197 done = 0;
198 break;
199 }
200 qfb_set(pow, twopow);
201 }
202 }
203
204 mult += 2;
205
206 qfb_hash_clear(qhash, depth);
207
208 }
209
210 iters *= 2;
211 jmax = (slong) (jmax * 1.15);
212
213 mult = iters/10;
214 mult |= 1;
215 }
216
217 qfb_clear(pow);
218 qfb_clear(twopow);
219
220 fmpz_clear(g);
221 fmpz_clear(n);
222 fmpz_clear(n0);
223 fmpz_clear(p);
224 fmpz_clear(L);
225
226 return 0;
227 }
228