flint-sys 0.9.0

Bindings to the FLINT C library
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

/*
    Copyright (C) 2012 William Hart

    This file is part of FLINT.

    FLINT is free software: you can redistribute it and/or modify it under
    the terms of the GNU Lesser General Public License (LGPL) as published
    by the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
*/

#include <mpfr.h>
#include "ulong_extras.h"
#include "qfb.h"

int qfb_exponent_grh(fmpz_t exponent, fmpz_t n, ulong B1, ulong B2_sqrt)
{
   fmpz_t p, exp, n2;
   mpz_t mn;
   qfb_t f;
   ulong pr, nmodpr, s, grh_limit;
   mpfr_t lim;
   int ret = 1;
   n_primes_t iter;

   n_primes_init(iter);
   fmpz_init(p);
   fmpz_init(n2);
   fmpz_init(exp);
   qfb_init(f);

   flint_mpz_init_set_readonly(mn, n);
   mpfr_init_set_z(lim, mn, MPFR_RNDA);
   mpfr_abs(lim, lim, MPFR_RNDU);
   mpfr_log(lim, lim, MPFR_RNDU);
   mpfr_mul(lim, lim, lim, MPFR_RNDU);
   mpfr_mul_ui(lim, lim, 6, MPFR_RNDU);
   grh_limit = mpfr_get_ui(lim, MPFR_RNDU);

   fmpz_set_ui(exponent, 1);

   /* find odd prime such that n is a square mod p */
   pr = 0;
   for (pr = 1; pr < grh_limit; )
   {
      do
      {
         pr = n_primes_next(iter);
         nmodpr = fmpz_fdiv_ui(n, pr);
      } while ((pr == 2 && ((s = fmpz_fdiv_ui(n, 8)) == 2 || s == 3 || s == 5))
         || (pr != 2 && nmodpr != 0 && n_jacobi(nmodpr, pr) < 0));

      if (pr < grh_limit)
      {
         fmpz_set_ui(p, pr);

         /* find prime form of discriminant n */
         qfb_prime_form(f, n, p);
         fmpz_set(n2, n);

         /* deal with non-fundamental discriminants */
         if (nmodpr == 0 && fmpz_fdiv_ui(f->c, pr) == 0)
         {
            fmpz_fdiv_q_ui(f->a, f->a, pr);
            fmpz_fdiv_q_ui(f->b, f->b, pr);
            fmpz_fdiv_q_ui(f->c, f->c, pr);
            fmpz_fdiv_q_ui(n2, n2, pr*pr);
         }
         if (pr == 2 && fmpz_is_even(f->a)
                     && fmpz_is_even(f->b) && fmpz_is_even(f->c))
         {
            fmpz_fdiv_q_2exp(f->a, f->a, 1);
            fmpz_fdiv_q_2exp(f->b, f->b, 1);
            fmpz_fdiv_q_2exp(f->c, f->c, 1);
            fmpz_fdiv_q_2exp(n2, n2, 2);
         }

         qfb_reduce(f, f, n2);

         if (!fmpz_is_one(exponent))
            qfb_pow(f, f, n2, exponent);

         if (!qfb_exponent_element(exp, f, n2, B1, B2_sqrt))
         {
            ret = 0;
            goto cleanup;
         }

         if (!fmpz_is_one(exp))
            fmpz_mul(exponent, exponent, exp);
      }
   }

cleanup:
   qfb_clear(f);
   fmpz_clear(p);
   fmpz_clear(n2);
   fmpz_clear(exp);
   n_primes_clear(iter);
   mpfr_clear(lim);

   flint_mpz_clear_readonly(mn);

   return ret;
}