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 "flint/flint.h"
20 #include "flint/ulong_extras.h"
21 #include "flint/fmpz.h"
22 #include "qfb.h"
23
qfb_prime_form(qfb_t r,fmpz_t D,fmpz_t p)24 void qfb_prime_form(qfb_t r, fmpz_t D, fmpz_t p)
25 {
26 fmpz_t q, rem, s, t;
27
28 fmpz_init(s);
29
30 if (fmpz_cmp_ui(p, 2) == 0) /* special case, p = 2 */
31 {
32 ulong m8 = fmpz_fdiv_ui(D, 8);
33 if (m8 == 4)
34 fmpz_set_ui(r->b, 2);
35 else
36 fmpz_set_ui(r->b, m8);
37
38 fmpz_sub_ui(s, D, m8);
39 fmpz_neg(s, s);
40 fmpz_fdiv_q_2exp(r->c, s, 3);
41 fmpz_set(r->a, p);
42
43 fmpz_clear(s);
44
45 return;
46 }
47
48 fmpz_init(t);
49
50 fmpz_mod(t, D, p);
51 if (fmpz_is_zero(t)) /* special case, p | D */
52 {
53 fmpz_init(q);
54 fmpz_init(rem);
55
56 fmpz_fdiv_q(s, D, p); /* s = D/p */
57 if (fmpz_is_zero(s))
58 fmpz_set(t, s);
59 else
60 fmpz_sub(t, p, s); /* t = -s mod p */
61 while (fmpz_fdiv_ui(t, 4) != 0) /* ensure t is divisible by 4 */
62 fmpz_add(t, t, p);
63 fmpz_add(t, t, s); /* t = first possible value for d/p + 4c */
64 fmpz_fdiv_q(t, t, p);
65 fmpz_sqrtrem(q, rem, t);
66 if (!fmpz_is_zero(rem)) /* t/p is not a square */
67 {
68 if (fmpz_is_even(t)) /* b/p is even as t/p is */
69 fmpz_add_ui(q, q, 1 + fmpz_is_even(q));
70 else
71 fmpz_add_ui(q, q, 1 + fmpz_is_odd(q));
72 } /* q = b/p */
73
74 fmpz_mul(r->b, q, p);
75 fmpz_mul(q, q, q);
76 fmpz_mul(q, q, p);
77 fmpz_sub(q, q, s);
78 fmpz_fdiv_q_2exp(r->c, q, 2);
79 fmpz_set(r->a, p);
80
81 fmpz_clear(q);
82 fmpz_clear(rem);
83 } else
84 {
85 fmpz_sqrtmod(t, t, p);
86 fmpz_sub(s, D, t);
87
88 if (fmpz_is_odd(s))
89 fmpz_sub(t, p, t);
90
91 fmpz_set(r->a, p);
92 fmpz_set(r->b, t);
93 fmpz_mul(r->c, r->b, r->b);
94 fmpz_sub(r->c, r->c, D);
95 fmpz_divexact(r->c, r->c, r->a);
96 fmpz_fdiv_q_2exp(r->c, r->c, 2);
97 }
98
99 fmpz_clear(s);
100 fmpz_clear(t);
101 }
102