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) 2010 Sebastian Pancratz
14 Copyright (C) 2011 Fredrik Johansson
15 Copyright (C) 2014 William Hart
16
17 ******************************************************************************/
18
19 #include "nf_elem.h"
20
_nf_elem_sub_lf(nf_elem_t a,const nf_elem_t b,const nf_elem_t c,const nf_t nf,int can)21 void _nf_elem_sub_lf(nf_elem_t a, const nf_elem_t b,
22 const nf_elem_t c, const nf_t nf, int can)
23 {
24 const fmpz * const p = LNF_ELEM_NUMREF(b);
25 const fmpz * const q = LNF_ELEM_DENREF(b);
26 const fmpz * const r = LNF_ELEM_NUMREF(c);
27 const fmpz * const s = LNF_ELEM_DENREF(c);
28 fmpz * const rnum = LNF_ELEM_NUMREF(a);
29 fmpz * const rden = LNF_ELEM_DENREF(a);
30 fmpz_t t;
31
32 if (can)
33 _fmpq_sub(rnum, rden, p, q, r, s);
34 else
35 {
36 /* Same denominator */
37 if (fmpz_equal(q, s))
38 {
39 fmpz_sub(rnum, p, r);
40 fmpz_set(rden, q);
41
42 return;
43 }
44
45 /* p/q is an integer */
46 if (fmpz_is_one(q))
47 {
48 fmpz_init(t);
49
50 fmpz_mul(t, p, s);
51 fmpz_sub(rnum, t, r);
52 fmpz_set(rden, s);
53
54 fmpz_clear(t);
55
56 return;
57 }
58
59 /* r/s is an integer */
60 if (fmpz_is_one(s))
61 {
62 fmpz_init(t);
63
64 fmpz_mul(t, r, q);
65 fmpz_sub(rnum, t, p);
66 fmpz_set(rden, q);
67
68 fmpz_clear(t);
69
70 return;
71 }
72
73 /*
74 We want to compute p/q - r/s which is (p*s - q*r, q*s).
75 */
76
77 fmpz_init(t);
78
79 fmpz_mul(t, q, r);
80 fmpz_mul(rnum, p, s);
81 fmpz_sub(rnum, rnum, t);
82 fmpz_mul(rden, q, s);
83
84 fmpz_clear(t);
85 }
86 }
87
_nf_elem_sub_qf(nf_elem_t a,const nf_elem_t b,const nf_elem_t c,const nf_t nf,int can)88 void _nf_elem_sub_qf(nf_elem_t a, const nf_elem_t b,
89 const nf_elem_t c, const nf_t nf, int can)
90 {
91 fmpz_t d;
92
93 const fmpz * const bnum = QNF_ELEM_NUMREF(b);
94 const fmpz * const bden = QNF_ELEM_DENREF(b);
95
96 const fmpz * const cnum = QNF_ELEM_NUMREF(c);
97 const fmpz * const cden = QNF_ELEM_DENREF(c);
98
99 fmpz * const anum = QNF_ELEM_NUMREF(a);
100 fmpz * const aden = QNF_ELEM_DENREF(a);
101
102 fmpz_init(d);
103 fmpz_one(d);
104
105 if (fmpz_equal(bden, cden))
106 {
107 fmpz_sub(anum, bnum, cnum);
108 fmpz_sub(anum + 1, bnum + 1, cnum + 1);
109 fmpz_sub(anum + 2, bnum + 2, cnum + 2);
110 fmpz_set(aden, bden);
111
112 if (can && !fmpz_is_one(aden))
113 {
114 fmpz_gcd(d, anum, anum + 1);
115 fmpz_gcd(d, d, anum + 2);
116 if (!fmpz_is_one(d))
117 {
118 fmpz_gcd(d, d, aden);
119
120 if (!fmpz_is_one(d))
121 {
122 fmpz_divexact(anum, anum, d);
123 fmpz_divexact(anum + 1, anum + 1, d);
124 fmpz_divexact(anum + 2, anum + 2, d);
125 fmpz_divexact(aden, aden, d);
126 }
127 }
128 }
129
130 fmpz_clear(d);
131
132 return;
133 }
134
135 if (!fmpz_is_one(bden) && !fmpz_is_one(cden))
136 fmpz_gcd(d, bden, cden);
137
138 if (fmpz_is_one(d))
139 {
140 fmpz_mul(anum, bnum, cden);
141 fmpz_mul(anum + 1, bnum + 1, cden);
142 fmpz_mul(anum + 2, bnum + 2, cden);
143 fmpz_submul(anum, cnum, bden);
144 fmpz_submul(anum + 1, cnum + 1, bden);
145 fmpz_submul(anum + 2, cnum + 2, bden);
146 fmpz_mul(aden, bden, cden);
147 } else
148 {
149 fmpz_t bden1;
150 fmpz_t cden1;
151
152 fmpz_init(bden1);
153 fmpz_init(cden1);
154
155 fmpz_divexact(bden1, bden, d);
156 fmpz_divexact(cden1, cden, d);
157
158 fmpz_mul(anum, bnum, cden1);
159 fmpz_mul(anum + 1, bnum + 1, cden1);
160 fmpz_mul(anum + 2, bnum + 2, cden1);
161 fmpz_submul(anum, cnum, bden1);
162 fmpz_submul(anum + 1, cnum + 1, bden1);
163 fmpz_submul(anum + 2, cnum + 2, bden1);
164
165 if (fmpz_is_zero(anum) && fmpz_is_zero(anum + 1) && fmpz_is_zero(anum + 2))
166 fmpz_one(aden);
167 else
168 {
169 if (can)
170 {
171 fmpz_t e;
172
173 fmpz_init(e);
174
175 fmpz_gcd(e, anum, anum + 1);
176 fmpz_gcd(e, e, anum + 2);
177 if (!fmpz_is_one(e))
178 fmpz_gcd(e, e, d);
179
180 if (fmpz_is_one(e))
181 fmpz_mul(aden, bden, cden1);
182 else
183 {
184 fmpz_divexact(anum, anum, e);
185 fmpz_divexact(anum + 1, anum + 1, e);
186 fmpz_divexact(anum + 2, anum + 2, e);
187 fmpz_divexact(bden1, bden, e);
188 fmpz_mul(aden, bden1, cden1);
189 }
190
191 fmpz_clear(e);
192 } else
193 fmpz_mul(aden, bden, cden1);
194 }
195
196 fmpz_clear(bden1);
197 fmpz_clear(cden1);
198 }
199
200 fmpz_clear(d);
201 }
202
nf_elem_sub_qf(nf_elem_t a,const nf_elem_t b,const nf_elem_t c,const nf_t nf)203 void nf_elem_sub_qf(nf_elem_t a, const nf_elem_t b,
204 const nf_elem_t c, const nf_t nf)
205 {
206 if (a == c)
207 {
208 nf_elem_t t;
209
210 nf_elem_init(t, nf);
211
212 _nf_elem_sub_qf(t, b, c, nf, 1);
213 nf_elem_swap(t, a, nf);
214
215 nf_elem_clear(t, nf);
216 } else
217 _nf_elem_sub_qf(a, b, c, nf, 1);
218 }
219