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 (C) 2012 William Hart
14
15 ******************************************************************************/
16
17 #include <stdlib.h>
18 #include <gmp.h>
19 #include "qfb.h"
20
21 /*
22 Iterate through all factors of a number given factorisation
23 into n prime powers whose maximum values are stored in exp,
24 storing the values at the current iteration in pows.
25 */
pow_incr(int * pows,int * exp,int n)26 int pow_incr(int * pows, int * exp, int n)
27 {
28 int i;
29
30 for (i = 0; i < n; i++)
31 {
32 pows[i]++;
33 if (pows[i] > exp[i])
34 pows[i] = 0;
35 else
36 return 1;
37 }
38
39 return 0;
40 }
41
qfb_reduced_forms_large(qfb ** forms,slong d)42 slong qfb_reduced_forms_large(qfb ** forms, slong d)
43 {
44 slong a, j, k, p, alim, alloc, num, roots, sqrt, i, prod, prime_i;
45 mp_srcptr primes;
46 const double * prime_inverses;
47 mp_limb_t a2;
48 n_factor_t * fac;
49
50 if (d >= 0)
51 {
52 printf("Exception: qfb_reduced_forms not implemented for positive discriminant.\n");
53 abort();
54 }
55
56 alim = n_sqrt(-d/3); /* maximum a value to check */
57
58 (*forms) = NULL;
59 alloc = 0;
60
61 if ((-d & 3) == 2 || (-d & 3) == 1) /* ensure d is 0, 1 mod 4 */
62 return 0;
63
64 fac = flint_calloc(alim + 1, sizeof(n_factor_t));
65
66 for (a = 2; a <= alim; a += 2) /* find powers of 2 dividing 4a values */
67 {
68 a2 = a;
69 fac[a].exp[0] = n_remove(&a2, 2) + 2;
70 fac[a].p[0] = 2;
71 fac[a].num = 1;
72 }
73
74 for (a = 1; a <= alim; a += 2)
75 {
76 fac[a].exp[0] = 2;
77 fac[a].p[0] = 2;
78 fac[a].num = 1;
79 }
80
81 sqrt = n_sqrt(alim);
82 primes = n_primes_arr_readonly(FLINT_MAX(sqrt, 10000));
83 prime_inverses = n_prime_inverses_arr_readonly(FLINT_MAX(sqrt, 10000));
84
85 prime_i = 1;
86 while ((p = primes[prime_i]) <= sqrt) /* sieve for factors of 4a values */
87 {
88 for (a = p; a <= alim; a+= p)
89 {
90 a2 = a;
91 num = fac[a].num;
92 fac[a].exp[num] = n_remove2_precomp(&a2, p, prime_inverses[prime_i]);
93 fac[a].p[num] = p;
94 fac[a].num++;
95 }
96 prime_i++;
97 }
98
99 for (a = 1; a <= alim; a++) /* write any remaining prime factor of each 4a value */
100 {
101 prod = 1;
102 for (i = 0; i < fac[a].num; i++)
103 prod *= n_pow(fac[a].p[i], fac[a].exp[i]);
104 p = (4*a)/prod;
105 if (p != 1)
106 {
107 num = fac[a].num;
108 fac[a].exp[num] = 1;
109 fac[a].p[num] = p;
110 fac[a].num++;
111 }
112 }
113
114 num = 0;
115
116 for (a = 1; a <= alim; a++) /* loop through possible a's */
117 {
118 mp_limb_t * s;
119
120 roots = n_sqrtmodn(&s, n_negmod((-d)%(4*a), 4*a), fac + a);
121
122 for (j = 0; j < roots; j++) /* loop through all square roots of d mod 4a */
123 {
124 mp_limb_signed_t b = s[j];
125
126 if (b > 2*a) b -= 4*a;
127
128 if (-a < b && b <= a) /* we may have a form */
129 {
130 /*
131 let B = FLINT_BITS
132 -sqrt(2^(B-1)/3) < b < sqrt(2^(B-1)/3)
133 0 < -d < 2^(B-1)
134 */
135 mp_limb_t c = ((mp_limb_t) (b*b) + (mp_limb_t) (-d))/(4*(mp_limb_t) a);
136
137 if (c >= (mp_limb_t) a && (b >= 0 || a != c)) /* we have a form */
138 {
139 mp_limb_t g;
140
141 if (b)
142 {
143 g = n_gcd(c, FLINT_ABS(b));
144 g = n_gcd(a, g);
145 } else
146 g = n_gcd(c, a);
147
148 if (g == 1) /* we have a primitive form */
149 {
150 if (num == alloc) /* realloc if necessary */
151 {
152 (*forms) = flint_realloc(*forms, (alloc + FLINT_MIN(alim, 100))*sizeof(qfb));
153 alloc += FLINT_MIN(alim, 100);
154 for (k = num; k < alloc; k++)
155 qfb_init((*forms) + k);
156 }
157
158 fmpz_set_si((*forms)[num].a, a); /* record form */
159 fmpz_set_si((*forms)[num].b, b);
160 fmpz_set_si((*forms)[num++].c, c);
161 }
162 }
163 }
164 }
165
166 flint_free(s);
167 }
168
169 flint_free(fac);
170
171 return num;
172 }
173
qfb_reduced_forms(qfb ** forms,slong d)174 slong qfb_reduced_forms(qfb ** forms, slong d)
175 {
176 slong a, b, k, c, p, blim, alloc, num, sqrt, i, prod, prime_i;
177 mp_srcptr primes;
178 const double * prime_inverses;
179 mp_limb_t b2, exp, primes_cutoff = 0;
180 n_factor_t * fac;
181 mp_limb_t * s;
182
183 if (d >= 0)
184 {
185 printf("Exception: qfb_reduced_forms not implemented for positive discriminant.\n");
186 abort();
187 }
188
189 blim = n_sqrt(-d/3); /* maximum a value to check */
190
191 (*forms) = NULL;
192 alloc = 0;
193
194 if ((-d & 3) == 2 || (-d & 3) == 1) /* ensure d is 0, 1 mod 4 */
195 return 0;
196
197 sqrt = n_sqrt(blim*blim - d);
198 n_nth_prime_bounds(&primes_cutoff, &primes_cutoff, sqrt);
199 if (primes_cutoff > FLINT_PRIMES_SMALL_CUTOFF*FLINT_PRIMES_SMALL_CUTOFF)
200 return qfb_reduced_forms_large(forms, d);
201
202 primes = n_primes_arr_readonly(FLINT_MAX(sqrt, 10000));
203 prime_inverses = n_prime_inverses_arr_readonly(FLINT_MAX(sqrt, 10000));
204
205 fac = flint_calloc(blim + 1, sizeof(n_factor_t));
206
207 prime_i = 1;
208 while ((p = primes[prime_i]) <= sqrt) /* sieve for factors of p^exp */
209 {
210 num = n_sqrtmod_primepow(&s, n_negmod((-d) % p, p), p, 1);
211
212 for (i = 0; i < num; i++) /* sieve with each sqrt mod p */
213 {
214 mp_limb_t off = s[i];
215 while (off <= blim)
216 {
217 b2 = (off*off - (mp_limb_t) d)/4;
218
219 fac[off].p[fac[off].num] = p;
220 fac[off].exp[fac[off].num] = n_remove2_precomp(&b2, p, prime_inverses[prime_i]);
221 fac[off].num++;
222
223 off += p;
224 }
225 }
226
227 prime_i++;
228
229 flint_free(s);
230 }
231
232 for (b = (d & 1); b <= blim; b += 2) /* write any remaining factors, including 2^exp */
233 {
234 b2 = ((mp_limb_t)(b*b - d))/4;
235
236 count_trailing_zeros(exp, b2); /* powers of 2 */
237 if (exp)
238 {
239 fac[b].p[fac[b].num] = 2;
240 fac[b].exp[fac[b].num] = exp;
241 fac[b].num++;
242 }
243
244 prod = 1;
245 for (i = 0; i < fac[b].num; i++)
246 prod *= n_pow(fac[b].p[i], fac[b].exp[i]);
247
248 b2 /= prod;
249
250 if (b2 != 1)
251 {
252 fac[b].p[fac[b].num] = b2;
253 fac[b].exp[fac[b].num] = 1;
254 fac[b].num++;
255 }
256 }
257
258 num = 0;
259 for (b = (d & 1); b <= blim; b += 2) /* compute forms for each b */
260 {
261 int pows[FLINT_MAX_FACTORS_IN_LIMB];
262 int n = fac[b].num;
263
264 b2 = ((mp_limb_t)(b*b - d))/4;
265
266 for (i = 0; i < n; i++)
267 pows[i] = 0;
268
269 do
270 {
271 a = 1;
272
273 for (i = 0; i < n; i++)
274 a *= n_pow(fac[b].p[i], pows[i]);
275
276 c = b2 / a;
277
278 if (a <= c && b <= a) /* we have a form */
279 {
280 mp_limb_t g;
281
282 if (b)
283 {
284 g = n_gcd(c, b);
285 g = n_gcd(a, g);
286 } else
287 g = n_gcd(c, a);
288
289 if (g == 1) /* primitive form */
290 {
291 if (num == alloc) /* realloc if necessary */
292 {
293 (*forms) = flint_realloc(*forms, (alloc + FLINT_MIN(blim, 100))*sizeof(qfb));
294 alloc += FLINT_MIN(blim, 100);
295 for (k = num; k < alloc; k++)
296 qfb_init((*forms) + k);
297 }
298
299 fmpz_set_si((*forms)[num].a, a); /* record form */
300 fmpz_set_si((*forms)[num].b, b);
301 fmpz_set_si((*forms)[num++].c, c);
302
303 if (b && a != c && b != a)
304 {
305 if (num == alloc) /* realloc if necessary */
306 {
307 (*forms) = flint_realloc(*forms, (alloc + FLINT_MIN(blim, 100))*sizeof(qfb));
308 alloc += FLINT_MIN(blim, 100);
309 for (k = num; k < alloc; k++)
310 qfb_init((*forms) + k);
311 }
312
313 fmpz_set_si((*forms)[num].a, a); /* record opposite form */
314 fmpz_set_si((*forms)[num].b, -b);
315 fmpz_set_si((*forms)[num++].c, c);
316 }
317 }
318 }
319 } while (pow_incr(pows, fac[b].exp, n));
320 }
321
322 flint_free(fac);
323
324 return num;
325 }
326