1 //==============================================================================================
2 //
3 // This file is part of LiDIA --- a library for computational number theory
4 //
5 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
6 //
7 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
8 //
9 //----------------------------------------------------------------------------------------------
10 //
11 // $Id$
12 //
13 // Author : Stefan Neis (SN)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifdef HAVE_CONFIG_H
20 # include "config.h"
21 #endif
22 #include "LiDIA/prime_ideal.h"
23 #include "LiDIA/debug.h"
24 #include <cctype>
25
26
27
28 #ifdef LIDIA_NAMESPACE
29 namespace LiDIA {
30 #endif
31
32
33
34 // Constructors & destructor
prime_ideal()35 prime_ideal::prime_ideal ()
36 : gen1(0),
37 gen2(),
38 e(0),
39 f(0),
40 valu(bigint(0))
41 {
42 debug_handler_l("prime_ideal", "constructor (void)", 1);
43 }
44
45
46
prime_ideal(const bigint & prim,const alg_number & a,lidia_size_t ram,lidia_size_t inert)47 prime_ideal::prime_ideal (const bigint & prim, const alg_number & a,
48 lidia_size_t ram, lidia_size_t inert)
49 : gen1(prim),
50 gen2(a),
51 e(ram),
52 f(inert),
53 valu(0, a.which_base())
54 {
55 debug_handler_l("prime_ideal", "constructor (const bigint &, "
56 "const alg_number &, lidia_size_t, lidia_size_t)", 1);
57 if (!is_prime(gen1, 10) || !a.denominator().is_one())
58 lidia_error_handler("prime_ideal", "constructor can only be called with\n"
59 "\t a prime number and\n"
60 "\t an integral algebraic number\n as arguments");
61
62 if (f == 0 || e == 0) {
63 module M = alg_ideal(gen1, gen2);
64 if (f == 0) {
65 f = a.degree() - M.coeff_matrix().get_no_of_columns();
66 }
67 if (e == 0) {
68 M %= gen1;
69 module P1;
70 module P2 = M;
71
72 do {
73 P1 = P2;
74 e++;
75 multiply(P2, P1, M);
76 P2 %= gen1;
77 } while (P2.coeff_matrix().get_no_of_columns() <
78 P1.coeff_matrix().get_no_of_columns() ||
79 (P2.coeff_matrix().is_column_zero(0) &&
80 !P1.coeff_matrix().is_column_zero(0)));
81 }
82 }
83 }
84
85
86
prime_ideal(const bigint & prim,const nf_base * O)87 prime_ideal::prime_ideal (const bigint & prim, const nf_base * O)
88 : gen1(prim),
89 gen2(O),
90 e(1),
91 f(O->degree()),
92 valu(0, O)
93 {
94 debug_handler_l("prime_ideal", "constructor (const bigint &, "
95 "const nf_base &)", 1);
96 if (!is_prime(gen1, 10))
97 lidia_error_handler("prime_ideal", "constructor can only be called with\n"
98 "\t a prime number and\n"
99 "\t (optionally) an nf_base/order/number_field");
100 }
101
102
103
104 //Computing the help variable for computing valuations(Cohen, p. 199-201)
105 void
compute_valu() const106 prime_ideal::compute_valu () const
107 {
108 debug_handler("prime_ideal", "in member-function compute_valu()");
109
110 // Set a to the second generator of the prime_ideal
111 lidia_size_t n = gen2.degree();
112 bigint fact;
113
114 bigmod_matrix LGS(n, n, gen1);
115
116 bigint *tmp = new bigint[n];
117
118 for (register lidia_size_t i = 0; i < n; tmp[i++] = 0) {
119 tmp[i] = 1; //tmp = e_i;
120 alg_number b (tmp, 1, gen2.which_base()); //b = w_i
121 multiply(b, b, gen2);
122 LGS.sto_column_vector(b.coeff_vector(), n, i);
123 }
124 delete[] tmp;
125
126 debug_handler_c("prime_ideal", "in member-function compute_valu()", 2,
127 std::cout << "Computing kernel of " << LGS << " mod " << gen1);
128 LGS.kernel(LGS, fact);
129 debug_handler_c("prime_ideal", "in member-function compute_valu()", 2,
130 std::cout << "kernel is " << LGS);
131
132 if (!fact.is_one()) {
133 lidia_error_handler("prime_ideal", "compute_valu()::"
134 "internal error 1 while computing kernel");
135 }
136 valu = alg_number(tmp = LGS.get_column(0), 1, gen2.which_base());
137 delete[] tmp;
138 }
139
140
141
142 // swap
143 void
swap(prime_ideal & a,prime_ideal & b)144 swap (prime_ideal &a, prime_ideal &b) {
145 swap(a.gen1, b.gen1);
146 swap(a.gen2, b.gen2);
147 swap(a.valu, b.valu);
148
149 lidia_size_t tmp;
150
151 tmp = a.e;
152 a.e = b.e;
153 b.e = tmp;
154
155 tmp = a.f;
156 a.f = b.f;
157 b.f = tmp;
158 }
159
160
161
162 // High-level functions:
163 long
ord(const prime_ideal & prim,const alg_ideal & MM)164 ord (const prime_ideal & prim, const alg_ideal & MM)
165 {
166 debug_handler_c("prime_ideal", "in ord(prime_ideal & alg_ideal &)", 3,
167 std::cout << "Computing valuation of " << MM << " at " << prim);
168 // const bigmod_matrix & Mcoeff = M.coeff_matrix();
169 // const bigint & Mmod = Mcoeff.get_modulus();
170 long v = 0;
171 if (prim.gen2.is_zero()) {
172 bigint res, tmp(MM.denominator());
173 div_rem(tmp, res, tmp, prim.gen1);
174 while (res.is_zero()) {
175 v --;
176 div_rem(tmp, res, tmp, prim.gen1);
177 }
178 bigmod_matrix M = MM.coeff_matrix();
179 if (v == 0)
180 // special algorithm for inert primes
181 while (divide(M, M, prim.gen1))
182 v++;
183 return v;
184 }
185 if (prim.valu.is_zero()) prim.compute_valu();
186
187 bigint res, tmp2(MM.denominator());
188
189 // Initialize with the valuation obtained from the denominator;
190 div_rem(tmp2, res, tmp2, prim.gen1);
191 while (res.is_zero()) {
192 v -= prim.e;
193 div_rem(tmp2, res, tmp2, prim.gen1);
194 }
195 debug_handler_c("prime_ideal", "in ord(prime_ideal &, alg_ideal &)", 3,
196 std::cout << "now v is " << v);
197
198 // check whether prim.gen1 divides all pseudo--diagonal elements.
199 remainder(res, MM.coeff_matrix().get_modulus(), prim.gen1);
200 if (!res.is_zero()) {
201 debug_handler_c("prime_ideal", "in ord(prime_ideal &, alg_ideal &)", 3,
202 std::cout << "returning(1) " << v << std::endl);
203 return v;
204 }
205
206 // Now add valuation of the numerator:
207 alg_number b;
208 alg_ideal M (MM);
209
210 debug_handler_c("prime_ideal", "in ord(prime_ideal &, alg_ideal &)", 3,
211 std::cout << "valu is " << prim.valu << std::endl);
212
213 do {
214 debug_handler_c("prime_ideal", "in ord(prime_ideal &, alg_ideal &)", 3,
215 std::cout << "At start of loop M is " << M << std::endl);
216 multiply(M, M, prim.valu);
217
218 debug_handler_c("prime_ideal", "in ord(prime_ideal &, alg_ideal &)", 3,
219 std::cout << "Now(1) M is " << M.base << std::endl);
220
221 // remainder(res, prim.O->base_denominator(), prim.p);
222 // if (!res.is_zero()){
223 // v++;
224 // std::cout << "now(1) v is "<<v<<","<<std::endl;
225 // std::cout << " since O is "<<(*prim.O)<<" i.e. field is ";
226 // std::cout << (*(prim.O->which_field())) << " with TRAFO";
227 // std::cout << prim.O->base_numerator()<<"/"<<prim.O->base_denominator();
228 // std::cout <<std::endl<<std::flush;
229 // divide(A, A, prim.p);
230 // std::cout << "Now(2) A is "<<A<<std::endl;
231 // continue;
232 // }
233 if (!divide(M.base, M.base, prim.gen1)) {
234 debug_handler_c("prime_ideal", "in ord(prime_ideal &, alg_ideal &)", 3,
235 std::cout << "returning(3) " << v << std::endl);
236 return v;
237 }
238 v++;
239 debug_handler_c("prime_ideal", "in ord(prime_ideal &, alg_ideal &)", 3,
240 std::cout << "Now(2) v is " << v << std::endl;
241 std::cout << "and M is " << M << std::endl);
242 } while (true);
243
244 }
245
246
247
248 long
ord(const prime_ideal & prim,const alg_number & a)249 ord (const prime_ideal & prim, const alg_number & a)
250 {
251 return ord(prim, a*order(prim.gen2.which_base()));
252 }
253
254
255
256 void
power(alg_ideal & c,const prime_ideal & a,const bigint & b)257 power (alg_ideal & c, const prime_ideal & a, const bigint & b)
258 {
259 debug_handler("prime_ideal", "in function power(alg_ideal &, "
260 "const prime_ideal &, const bigint &)");
261
262 const nf_base * O = a.gen2.which_base();
263 if (b.is_negative())
264 power(c, inverse(alg_ideal(alg_number(a.gen1, a.gen2.which_base()),
265 a.gen2)), -b);
266 else if (b.is_zero())
267 c = order(O);
268 else {
269 bigint gen1;
270 bigint expo;
271 div_rem(expo, gen1, b, a.e);
272 if (!gen1.is_zero()) ++expo;
273 power(gen1, a.gen1, expo);
274 if (!a.gen2.is_zero()) {
275 alg_number gen2;
276 power(gen2, a.gen2, b);
277 c = alg_ideal(gen1, gen2);
278 }
279 else
280 c = alg_ideal(gen1, alg_number(0, O));
281 }
282 }
283
284
285 // Input/Output:
286 std::ostream &
operator <<(std::ostream & s,const prime_ideal & a)287 operator << (std::ostream & s, const prime_ideal & a)
288 {
289 s << "<";
290 s << a.gen1;
291 if (!a.gen2.is_zero()) {
292 s << ", ";
293 s << a.gen2;
294 }
295 s << ">" << std::flush;
296 return s;
297 }
298
299
300
301 std::istream &
operator >>(std::istream & s,prime_ideal & p)302 operator >> (std::istream & s, prime_ideal & p)
303 {
304 bigint a;
305 alg_number b;
306 char c;
307
308 do {
309 s.get(c);
310 } while (isspace(c));
311 if (c != '<')
312 lidia_error_handler("prime_ideal", "operator >>:: A prime ideal must be"
313 "of form `<prime, alg_number>'");
314 s >> a;
315
316 do {
317 s.get(c);
318 } while (isspace(c));
319 if (c != ',')
320 lidia_error_handler("prime_ideal", "operator >>:: A prime ideal must be"
321 "of form `<prime, alg_number>'");
322
323 s >> b;
324 do {
325 s.get(c);
326 } while (isspace(c));
327 if (c != '>')
328 lidia_error_handler("prime_ideal", "operator >>:: A prime ideal must be"
329 "of form `<prime, alg_number>'");
330
331 p = prime_ideal(a, b);
332 return s;
333 }
334
335
336
337 #ifdef LIDIA_NAMESPACE
338 } // end of namespace LiDIA
339 #endif
340