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 <mpfr.h>
20 #include "qfb.h"
21
qfb_exponent_grh(fmpz_t exponent,fmpz_t n,ulong B1,ulong B2_sqrt)22 int qfb_exponent_grh(fmpz_t exponent, fmpz_t n, ulong B1, ulong B2_sqrt)
23 {
24 fmpz_t p, exp, n2;
25 mpz_t mn;
26 qfb_t f;
27 ulong pr, nmodpr, s, grh_limit;
28 mpfr_t lim;
29 int ret = 1;
30 n_primes_t iter;
31
32 n_primes_init(iter);
33 fmpz_init(p);
34 fmpz_init(n2);
35 fmpz_init(exp);
36 qfb_init(f);
37
38 flint_mpz_init_set_readonly(mn, n);
39 mpfr_init_set_z(lim, mn, MPFR_RNDA);
40 mpfr_abs(lim, lim, MPFR_RNDU);
41 mpfr_log(lim, lim, MPFR_RNDU);
42 mpfr_mul(lim, lim, lim, MPFR_RNDU);
43 mpfr_mul_ui(lim, lim, 6, MPFR_RNDU);
44 grh_limit = mpfr_get_ui(lim, MPFR_RNDU);
45
46 fmpz_set_ui(exponent, 1);
47
48 /* find odd prime such that n is a square mod p */
49 pr = 0;
50 for (pr = 1; pr < grh_limit; )
51 {
52 do
53 {
54 pr = n_primes_next(iter);
55 nmodpr = fmpz_fdiv_ui(n, pr);
56 } while ((pr == 2 && ((s = fmpz_fdiv_ui(n, 8)) == 2 || s == 3 || s == 5))
57 || (pr != 2 && nmodpr != 0 && n_jacobi(nmodpr, pr) < 0));
58
59 if (pr < grh_limit)
60 {
61 fmpz_set_ui(p, pr);
62
63 /* find prime form of discriminant n */
64 qfb_prime_form(f, n, p);
65 fmpz_set(n2, n);
66
67 /* deal with non-fundamental discriminants */
68 if (nmodpr == 0 && fmpz_fdiv_ui(f->c, pr) == 0)
69 {
70 fmpz_fdiv_q_ui(f->a, f->a, pr);
71 fmpz_fdiv_q_ui(f->b, f->b, pr);
72 fmpz_fdiv_q_ui(f->c, f->c, pr);
73 fmpz_fdiv_q_ui(n2, n2, pr*pr);
74 }
75 if (pr == 2 && fmpz_is_even(f->a)
76 && fmpz_is_even(f->b) && fmpz_is_even(f->c))
77 {
78 fmpz_fdiv_q_2exp(f->a, f->a, 1);
79 fmpz_fdiv_q_2exp(f->b, f->b, 1);
80 fmpz_fdiv_q_2exp(f->c, f->c, 1);
81 fmpz_fdiv_q_2exp(n2, n2, 2);
82 }
83
84 qfb_reduce(f, f, n2);
85
86 if (!fmpz_is_one(exponent))
87 qfb_pow(f, f, n2, exponent);
88
89 if (!qfb_exponent_element(exp, f, n2, B1, B2_sqrt))
90 {
91 ret = 0;
92 goto cleanup;
93 }
94
95 if (!fmpz_is_one(exp))
96 fmpz_mul(exponent, exponent, exp);
97 }
98 }
99
100 cleanup:
101 qfb_clear(f);
102 fmpz_clear(p);
103 fmpz_clear(n2);
104 fmpz_clear(exp);
105 n_primes_clear(iter);
106
107 flint_mpz_clear_readonly(mn);
108
109 return ret;
110 }
111