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 "qfb.h"
21
22 /*
23 This is based (with permission) on an implementation of NUCOMP
24 by Michael Jacobson which is an adaption to binary quadratic
25 forms of the algorithm presented in Appendix A of "Solving the
26 Pell Equation", Michael Jacobson and High Williams, CMS Books in
27 Mathematics, Springer 2009.
28 */
qfb_nucomp(qfb_t r,const qfb_t f,const qfb_t g,fmpz_t D,fmpz_t L)29 void qfb_nucomp(qfb_t r, const qfb_t f, const qfb_t g, fmpz_t D, fmpz_t L)
30 {
31 fmpz_t a1, a2, c2, ca, cb, cc, k, s, sp, ss, m, t, u2, v1, v2;
32
33 if (fmpz_cmp(f->a, g->a) > 0)
34 {
35 qfb_nucomp(r, g, f, D, L);
36 return;
37 }
38
39 fmpz_init(a1); fmpz_init(a2); fmpz_init(c2);
40 fmpz_init(ca); fmpz_init(cb); fmpz_init(cc);
41 fmpz_init(k); fmpz_init(m);
42 fmpz_init(s); fmpz_init(sp); fmpz_init(ss);
43 fmpz_init(t); fmpz_init(u2); fmpz_init(v1); fmpz_init(v2);
44
45 /* nucomp calculation */
46
47 fmpz_set(a1, f->a);
48 fmpz_set(a2, g->a);
49 fmpz_set(c2, g->c);
50
51 fmpz_add(ss, f->b, g->b);
52 fmpz_fdiv_q_2exp(ss, ss, 1);
53
54 fmpz_sub(m, f->b, g->b);
55 fmpz_fdiv_q_2exp(m, m, 1);
56
57 fmpz_fdiv_r(t, a2, a1);
58 if (fmpz_is_zero(t))
59 {
60 fmpz_set_ui(v1, 0);
61 fmpz_set(sp, a1);
62 } else
63 fmpz_gcdinv(sp, v1, t, a1);
64
65 fmpz_mul(k, m, v1);
66 fmpz_fdiv_r(k, k, a1);
67
68 if (!fmpz_is_one(sp))
69 {
70 fmpz_xgcd(s, v2, u2, ss, sp);
71
72 fmpz_mul(k, k, u2);
73 fmpz_mul(t, v2, c2);
74 fmpz_sub(k, k, t);
75
76 if (!fmpz_is_one(s))
77 {
78 fmpz_divexact(a1, a1, s);
79 fmpz_divexact(a2, a2, s);
80 fmpz_mul(c2, c2, s);
81 }
82
83 fmpz_fdiv_r(k, k, a1);
84 }
85
86 if (fmpz_cmp(a1, L) < 0)
87 {
88 fmpz_mul(t, a2, k);
89
90 fmpz_mul(ca, a2, a1);
91
92 fmpz_mul_2exp(cb, t, 1);
93 fmpz_add(cb, cb, g->b);
94
95 fmpz_add(cc, g->b, t);
96 fmpz_mul(cc, cc, k);
97 fmpz_add(cc, cc, c2);
98
99 fmpz_divexact(cc, cc, a1);
100 } else
101 {
102 fmpz_t m1, m2, r1, r2, co1, co2, temp;
103
104 fmpz_init(m1); fmpz_init(m2); fmpz_init(r1); fmpz_init(r2);
105 fmpz_init(co1); fmpz_init(co2); fmpz_init(temp);
106
107 fmpz_set(r2, a1);
108 fmpz_set(r1, k);
109
110 fmpz_xgcd_partial(co2, co1, r2, r1, L);
111
112 fmpz_mul(t, a2, r1);
113 fmpz_mul(m1, m, co1);
114 fmpz_add(m1, m1, t);
115 fmpz_divexact(m1, m1, a1);
116
117 fmpz_mul(m2, ss, r1);
118 fmpz_mul(temp, c2, co1);
119 fmpz_sub(m2, m2, temp);
120 fmpz_divexact(m2, m2, a1);
121
122 fmpz_mul(ca, r1, m1);
123 fmpz_mul(temp, co1, m2);
124 if (fmpz_sgn(co1) < 0)
125 fmpz_sub(ca, ca, temp);
126 else
127 fmpz_sub(ca, temp, ca);
128
129 fmpz_mul(cb, ca, co2);
130 fmpz_sub(cb, t, cb);
131 fmpz_mul_2exp(cb, cb, 1);
132 fmpz_divexact(cb, cb, co1);
133 fmpz_sub(cb, cb, g->b);
134 fmpz_mul_2exp(temp, ca, 1);
135 fmpz_fdiv_r(cb, cb, temp);
136
137 fmpz_mul(cc, cb, cb);
138 fmpz_sub(cc, cc, D);
139 fmpz_divexact(cc, cc, ca);
140 fmpz_fdiv_q_2exp(cc, cc, 2);
141
142 if (fmpz_sgn(ca) < 0)
143 {
144 fmpz_neg(ca, ca);
145 fmpz_neg(cc, cc);
146 }
147
148 fmpz_clear(m1); fmpz_clear(m2); fmpz_clear(r1); fmpz_clear(r2);
149 fmpz_clear(co1); fmpz_clear(co2); fmpz_clear(temp);
150 }
151
152 fmpz_set(r->a, ca);
153 fmpz_set(r->b, cb);
154 fmpz_set(r->c, cc);
155
156 fmpz_clear(ca); fmpz_clear(cb); fmpz_clear(cc);
157 fmpz_clear(k); fmpz_clear(m);
158 fmpz_clear(s); fmpz_clear(sp); fmpz_clear(ss);
159 fmpz_clear(t); fmpz_clear(u2); fmpz_clear(v1); fmpz_clear(v2);
160 fmpz_clear(a1); fmpz_clear(a2); fmpz_clear(c2);
161 }
162