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 Copyright (C) 2020 Chia Network Inc
15
16 ******************************************************************************/
17
18 #include <stdlib.h>
19 #include <gmp.h>
20 #include "flint/flint.h"
21 #include "flint/ulong_extras.h"
22 #include "flint/fmpz.h"
23 #include "qfb.h"
24
25 static void
qfb_nudupl_gcdinv(fmpz_t d,fmpz_t a,fmpz_t t,const fmpz_t f,const fmpz_t g)26 qfb_nudupl_gcdinv(fmpz_t d, fmpz_t a, fmpz_t t, const fmpz_t f, const fmpz_t g)
27 {
28 if (fmpz_cmp(f, g) < 0)
29 fmpz_gcdinv(d, a, f, g);
30 else
31 {
32 fmpz_fdiv_r(t, f, g);
33 fmpz_gcdinv(d, a, t, g);
34 }
35 }
36
qfb_nudupl(qfb_t r,const qfb_t f,fmpz_t D,fmpz_t L)37 void qfb_nudupl(qfb_t r, const qfb_t f, fmpz_t D, fmpz_t L)
38 {
39 fmpz_t a1, b1, c1, ca, cb, cc, k, s, t, u2, v1, v2;
40
41 fmpz_init(a1); fmpz_init(b1); fmpz_init(c1);
42 fmpz_init(ca); fmpz_init(cb); fmpz_init(cc);
43 fmpz_init(k);
44 fmpz_init(s);
45 fmpz_init(t); fmpz_init(u2); fmpz_init(v1); fmpz_init(v2);
46
47 /* nucomp calculation */
48
49 fmpz_set(a1, f->a);
50 fmpz_set(c1, f->c);
51
52 fmpz_zero(k);
53
54 if (fmpz_cmpabs(b1, a1) == 0)
55 {
56 fmpz_set(s, a1);
57 fmpz_zero(v2);
58 } else if (fmpz_sgn(f->b) < 0)
59 {
60 fmpz_neg(b1, f->b);
61 qfb_nudupl_gcdinv(s, v2, t, b1, a1);
62 fmpz_neg(v2, v2);
63 } else
64 qfb_nudupl_gcdinv(s, v2, t, f->b, a1);
65
66 fmpz_mul(t, v2, c1);
67 fmpz_neg(k, t);
68
69 if (!fmpz_is_one(s))
70 {
71 fmpz_divexact(a1, a1, s);
72 fmpz_mul(c1, c1, s);
73 }
74
75 fmpz_fdiv_r(k, k, a1);
76
77 if (fmpz_cmp(a1, L) < 0)
78 {
79 fmpz_mul(t, a1, k);
80
81 fmpz_mul(ca, a1, a1);
82
83 fmpz_mul_2exp(cb, t, 1);
84 fmpz_add(cb, cb, f->b);
85
86 fmpz_add(cc, f->b, t);
87 fmpz_mul(cc, cc, k);
88 fmpz_add(cc, cc, c1);
89
90 fmpz_divexact(cc, cc, a1);
91 } else
92 {
93 fmpz_t m2, r1, r2, co1, co2, temp;
94
95 fmpz_init(m2); fmpz_init(r1); fmpz_init(r2);
96 fmpz_init(co1); fmpz_init(co2); fmpz_init(temp);
97
98 fmpz_set(r2, a1);
99 fmpz_set(r1, k);
100
101 fmpz_xgcd_partial(co2, co1, r2, r1, L);
102
103 fmpz_mul(t, a1, r1);
104
105 fmpz_mul(m2, f->b, r1);
106 fmpz_mul(temp, c1, co1);
107 fmpz_sub(m2, m2, temp);
108 fmpz_divexact(m2, m2, a1);
109
110 fmpz_mul(ca, r1, r1);
111 fmpz_mul(temp, co1, m2);
112 if (fmpz_sgn(co1) < 0)
113 fmpz_sub(ca, ca, temp);
114 else
115 fmpz_sub(ca, temp, ca);
116
117 fmpz_mul(cb, ca, co2);
118 fmpz_sub(cb, t, cb);
119 fmpz_mul_2exp(cb, cb, 1);
120 fmpz_divexact(cb, cb, co1);
121 fmpz_sub(cb, cb, f->b);
122 fmpz_mul_2exp(temp, ca, 1);
123 fmpz_fdiv_r(cb, cb, temp);
124
125 fmpz_mul(cc, cb, cb);
126 fmpz_sub(cc, cc, D);
127 fmpz_divexact(cc, cc, ca);
128 fmpz_fdiv_q_2exp(cc, cc, 2);
129
130 if (fmpz_sgn(ca) < 0)
131 {
132 fmpz_neg(ca, ca);
133 fmpz_neg(cc, cc);
134 }
135
136 fmpz_clear(m2); fmpz_clear(r1); fmpz_clear(r2);
137 fmpz_clear(co1); fmpz_clear(co2); fmpz_clear(temp);
138 }
139
140 fmpz_set(r->a, ca);
141 fmpz_set(r->b, cb);
142 fmpz_set(r->c, cc);
143
144 fmpz_clear(ca); fmpz_clear(cb); fmpz_clear(cc);
145 fmpz_clear(k);
146 fmpz_clear(s);
147 fmpz_clear(t); fmpz_clear(u2); fmpz_clear(v2);
148 fmpz_clear(a1); fmpz_clear(b1); fmpz_clear(c1);
149 }
150