1 | /**************************************** |
---|
2 | * Computer Algebra System SINGULAR * |
---|
3 | ****************************************/ |
---|
4 | /** |
---|
5 | * ABSTRACT: numbers in an algebraic extension field K[a] / < f(a) > |
---|
6 | * Assuming that we have a coeffs object cf, then these numbers |
---|
7 | * are polynomials in the polynomial ring K[a] represented by |
---|
8 | * cf->extRing. |
---|
9 | * IMPORTANT ASSUMPTIONS: |
---|
10 | * 1.) So far we assume that cf->extRing is a valid polynomial |
---|
11 | * ring in exactly one variable, i.e., K[a], where K is allowed |
---|
12 | * to be any field (representable in SINGULAR and which may |
---|
13 | * itself be some extension field, thus allowing for extension |
---|
14 | * towers). |
---|
15 | * 2.) Moreover, this implementation assumes that |
---|
16 | * cf->extRing->qideal is not NULL but an ideal with at |
---|
17 | * least one non-zero generator which may be accessed by |
---|
18 | * cf->extRing->qideal->m[0] and which represents the minimal |
---|
19 | * polynomial f(a) of the extension variable 'a' in K[a]. |
---|
20 | * 3.) As soon as an std method for polynomial rings becomes |
---|
21 | * availabe, all reduction steps modulo f(a) should be replaced |
---|
22 | * by a call to std. Moreover, in this situation one can finally |
---|
23 | * move from K[a] / < f(a) > to |
---|
24 | * K[a_1, ..., a_s] / I, with I some zero-dimensional ideal |
---|
25 | * in K[a_1, ..., a_s] given by a lex |
---|
26 | * Gröbner basis. |
---|
27 | * The code in algext.h and algext.cc is then capable of |
---|
28 | * computing in K[a_1, ..., a_s] / I. |
---|
29 | **/ |
---|
30 | |
---|
31 | |
---|
32 | |
---|
33 | |
---|
34 | #include <misc/auxiliary.h> |
---|
35 | |
---|
36 | #include <omalloc/omalloc.h> |
---|
37 | |
---|
38 | #include <reporter/reporter.h> |
---|
39 | |
---|
40 | #include <coeffs/coeffs.h> |
---|
41 | #include <coeffs/numbers.h> |
---|
42 | #include <coeffs/longrat.h> |
---|
43 | |
---|
44 | #include <polys/monomials/ring.h> |
---|
45 | #include <polys/monomials/p_polys.h> |
---|
46 | #include <polys/simpleideals.h> |
---|
47 | |
---|
48 | #include <polys/PolyEnumerator.h> |
---|
49 | |
---|
50 | #include <factory/factory.h> |
---|
51 | #include <polys/clapconv.h> |
---|
52 | #include <polys/clapsing.h> |
---|
53 | #include <polys/prCopy.h> |
---|
54 | |
---|
55 | #include <polys/ext_fields/algext.h> |
---|
56 | #define TRANSEXT_PRIVATES 1 |
---|
57 | #include <polys/ext_fields/transext.h> |
---|
58 | |
---|
59 | #ifdef LDEBUG |
---|
60 | #define naTest(a) naDBTest(a,__FILE__,__LINE__,cf) |
---|
61 | BOOLEAN naDBTest(number a, const char *f, const int l, const coeffs r); |
---|
62 | #else |
---|
63 | #define naTest(a) do {} while (0) |
---|
64 | #endif |
---|
65 | |
---|
66 | /// Our own type! |
---|
67 | static const n_coeffType ID = n_algExt; |
---|
68 | |
---|
69 | /* polynomial ring in which our numbers live */ |
---|
70 | #define naRing cf->extRing |
---|
71 | |
---|
72 | /* coeffs object in which the coefficients of our numbers live; |
---|
73 | * methods attached to naCoeffs may be used to compute with the |
---|
74 | * coefficients of our numbers, e.g., use naCoeffs->nAdd to add |
---|
75 | * coefficients of our numbers */ |
---|
76 | #define naCoeffs cf->extRing->cf |
---|
77 | |
---|
78 | /* minimal polynomial */ |
---|
79 | #define naMinpoly naRing->qideal->m[0] |
---|
80 | |
---|
81 | /// forward declarations |
---|
82 | BOOLEAN naGreaterZero(number a, const coeffs cf); |
---|
83 | BOOLEAN naGreater(number a, number b, const coeffs cf); |
---|
84 | BOOLEAN naEqual(number a, number b, const coeffs cf); |
---|
85 | BOOLEAN naIsOne(number a, const coeffs cf); |
---|
86 | BOOLEAN naIsMOne(number a, const coeffs cf); |
---|
87 | BOOLEAN naIsZero(number a, const coeffs cf); |
---|
88 | number naInit(long i, const coeffs cf); |
---|
89 | int naInt(number &a, const coeffs cf); |
---|
90 | number naNeg(number a, const coeffs cf); |
---|
91 | number naInvers(number a, const coeffs cf); |
---|
92 | number naAdd(number a, number b, const coeffs cf); |
---|
93 | number naSub(number a, number b, const coeffs cf); |
---|
94 | number naMult(number a, number b, const coeffs cf); |
---|
95 | number naDiv(number a, number b, const coeffs cf); |
---|
96 | void naPower(number a, int exp, number *b, const coeffs cf); |
---|
97 | number naCopy(number a, const coeffs cf); |
---|
98 | void naWriteLong(number &a, const coeffs cf); |
---|
99 | void naWriteShort(number &a, const coeffs cf); |
---|
100 | number naGetDenom(number &a, const coeffs cf); |
---|
101 | number naGetNumerator(number &a, const coeffs cf); |
---|
102 | number naGcd(number a, number b, const coeffs cf); |
---|
103 | //number naLcm(number a, number b, const coeffs cf); |
---|
104 | int naSize(number a, const coeffs cf); |
---|
105 | void naDelete(number *a, const coeffs cf); |
---|
106 | void naCoeffWrite(const coeffs cf, BOOLEAN details); |
---|
107 | //number naIntDiv(number a, number b, const coeffs cf); |
---|
108 | const char * naRead(const char *s, number *a, const coeffs cf); |
---|
109 | |
---|
110 | static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param); |
---|
111 | |
---|
112 | |
---|
113 | /// returns NULL if p == NULL, otherwise makes p monic by dividing |
---|
114 | /// by its leading coefficient (only done if this is not already 1); |
---|
115 | /// this assumes that we are over a ground field so that division |
---|
116 | /// is well-defined; |
---|
117 | /// modifies p |
---|
118 | // void p_Monic(poly p, const ring r); |
---|
119 | |
---|
120 | /// assumes that p and q are univariate polynomials in r, |
---|
121 | /// mentioning the same variable; |
---|
122 | /// assumes a global monomial ordering in r; |
---|
123 | /// assumes that not both p and q are NULL; |
---|
124 | /// returns the gcd of p and q; |
---|
125 | /// leaves p and q unmodified |
---|
126 | // poly p_Gcd(const poly p, const poly q, const ring r); |
---|
127 | |
---|
128 | /* returns NULL if p == NULL, otherwise makes p monic by dividing |
---|
129 | by its leading coefficient (only done if this is not already 1); |
---|
130 | this assumes that we are over a ground field so that division |
---|
131 | is well-defined; |
---|
132 | modifies p */ |
---|
133 | static inline void p_Monic(poly p, const ring r) |
---|
134 | { |
---|
135 | if (p == NULL) return; |
---|
136 | number n = n_Init(1, r->cf); |
---|
137 | if (p->next==NULL) { p_SetCoeff(p,n,r); return; } |
---|
138 | poly pp = p; |
---|
139 | number lc = p_GetCoeff(p, r); |
---|
140 | if (n_IsOne(lc, r->cf)) return; |
---|
141 | number lcInverse = n_Invers(lc, r->cf); |
---|
142 | p_SetCoeff(p, n, r); // destroys old leading coefficient! |
---|
143 | pIter(p); |
---|
144 | while (p != NULL) |
---|
145 | { |
---|
146 | number n = n_Mult(p_GetCoeff(p, r), lcInverse, r->cf); |
---|
147 | n_Normalize(n,r->cf); |
---|
148 | p_SetCoeff(p, n, r); // destroys old leading coefficient! |
---|
149 | pIter(p); |
---|
150 | } |
---|
151 | n_Delete(&lcInverse, r->cf); |
---|
152 | p = pp; |
---|
153 | } |
---|
154 | |
---|
155 | /// see p_Gcd; |
---|
156 | /// additional assumption: deg(p) >= deg(q); |
---|
157 | /// must destroy p and q (unless one of them is returned) |
---|
158 | static inline poly p_GcdHelper(poly &p, poly &q, const ring r) |
---|
159 | { |
---|
160 | while (q != NULL) |
---|
161 | { |
---|
162 | p_PolyDiv(p, q, FALSE, r); |
---|
163 | // swap p and q: |
---|
164 | poly& t = q; |
---|
165 | q = p; |
---|
166 | p = t; |
---|
167 | |
---|
168 | } |
---|
169 | return p; |
---|
170 | } |
---|
171 | |
---|
172 | /* assumes that p and q are univariate polynomials in r, |
---|
173 | mentioning the same variable; |
---|
174 | assumes a global monomial ordering in r; |
---|
175 | assumes that not both p and q are NULL; |
---|
176 | returns the gcd of p and q; |
---|
177 | leaves p and q unmodified */ |
---|
178 | static inline poly p_Gcd(const poly p, const poly q, const ring r) |
---|
179 | { |
---|
180 | assume((p != NULL) || (q != NULL)); |
---|
181 | |
---|
182 | poly a = p; poly b = q; |
---|
183 | if (p_Deg(a, r) < p_Deg(b, r)) { a = q; b = p; } |
---|
184 | a = p_Copy(a, r); b = p_Copy(b, r); |
---|
185 | |
---|
186 | /* We have to make p monic before we return it, so that if the |
---|
187 | gcd is a unit in the ground field, we will actually return 1. */ |
---|
188 | a = p_GcdHelper(a, b, r); |
---|
189 | p_Monic(a, r); |
---|
190 | return a; |
---|
191 | } |
---|
192 | |
---|
193 | /* see p_ExtGcd; |
---|
194 | additional assumption: deg(p) >= deg(q); |
---|
195 | must destroy p and q (unless one of them is returned) */ |
---|
196 | static inline poly p_ExtGcdHelper(poly &p, poly &pFactor, poly &q, poly &qFactor, |
---|
197 | ring r) |
---|
198 | { |
---|
199 | if (q == NULL) |
---|
200 | { |
---|
201 | qFactor = NULL; |
---|
202 | pFactor = p_ISet(1, r); |
---|
203 | p_SetCoeff(pFactor, n_Invers(p_GetCoeff(p, r), r->cf), r); |
---|
204 | p_Monic(p, r); |
---|
205 | return p; |
---|
206 | } |
---|
207 | else |
---|
208 | { |
---|
209 | poly pDivQ = p_PolyDiv(p, q, TRUE, r); |
---|
210 | poly ppFactor = NULL; poly qqFactor = NULL; |
---|
211 | poly theGcd = p_ExtGcdHelper(q, qqFactor, p, ppFactor, r); |
---|
212 | pFactor = ppFactor; |
---|
213 | qFactor = p_Add_q(qqFactor, |
---|
214 | p_Neg(p_Mult_q(pDivQ, p_Copy(ppFactor, r), r), r), |
---|
215 | r); |
---|
216 | return theGcd; |
---|
217 | } |
---|
218 | } |
---|
219 | |
---|
220 | |
---|
221 | /* assumes that p and q are univariate polynomials in r, |
---|
222 | mentioning the same variable; |
---|
223 | assumes a global monomial ordering in r; |
---|
224 | assumes that not both p and q are NULL; |
---|
225 | returns the gcd of p and q; |
---|
226 | moreover, afterwards pFactor and qFactor contain appropriate |
---|
227 | factors such that gcd(p, q) = p * pFactor + q * qFactor; |
---|
228 | leaves p and q unmodified */ |
---|
229 | poly p_ExtGcd(poly p, poly &pFactor, poly q, poly &qFactor, ring r) |
---|
230 | { |
---|
231 | assume((p != NULL) || (q != NULL)); |
---|
232 | poly a = p; poly b = q; BOOLEAN aCorrespondsToP = TRUE; |
---|
233 | if (p_Deg(a, r) < p_Deg(b, r)) |
---|
234 | { a = q; b = p; aCorrespondsToP = FALSE; } |
---|
235 | a = p_Copy(a, r); b = p_Copy(b, r); |
---|
236 | poly aFactor = NULL; poly bFactor = NULL; |
---|
237 | poly theGcd = p_ExtGcdHelper(a, aFactor, b, bFactor, r); |
---|
238 | if (aCorrespondsToP) { pFactor = aFactor; qFactor = bFactor; } |
---|
239 | else { pFactor = bFactor; qFactor = aFactor; } |
---|
240 | return theGcd; |
---|
241 | } |
---|
242 | |
---|
243 | |
---|
244 | |
---|
245 | #ifdef LDEBUG |
---|
246 | BOOLEAN naDBTest(number a, const char *f, const int l, const coeffs cf) |
---|
247 | { |
---|
248 | if (a == NULL) return TRUE; |
---|
249 | p_Test((poly)a, naRing); |
---|
250 | if (getCoeffType(cf)==n_algExt) |
---|
251 | { |
---|
252 | if((((poly)a)!=naMinpoly) |
---|
253 | && p_Totaldegree((poly)a, naRing) >= p_Totaldegree(naMinpoly, naRing) |
---|
254 | && (p_Totaldegree((poly)a, naRing)> 1)) // allow to output par(1) |
---|
255 | { |
---|
256 | dReportError("deg >= deg(minpoly) in %s:%d\n",f,l); |
---|
257 | return FALSE; |
---|
258 | } |
---|
259 | } |
---|
260 | return TRUE; |
---|
261 | } |
---|
262 | #endif |
---|
263 | |
---|
264 | void heuristicReduce(poly &p, poly reducer, const coeffs cf); |
---|
265 | void definiteReduce(poly &p, poly reducer, const coeffs cf); |
---|
266 | |
---|
267 | /* returns the bottom field in this field extension tower; if the tower |
---|
268 | is flat, i.e., if there is no extension, then r itself is returned; |
---|
269 | as a side-effect, the counter 'height' is filled with the height of |
---|
270 | the extension tower (in case the tower is flat, 'height' is zero) */ |
---|
271 | static coeffs nCoeff_bottom(const coeffs r, int &height) |
---|
272 | { |
---|
273 | assume(r != NULL); |
---|
274 | coeffs cf = r; |
---|
275 | height = 0; |
---|
276 | while (nCoeff_is_Extension(cf)) |
---|
277 | { |
---|
278 | assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL); |
---|
279 | cf = cf->extRing->cf; |
---|
280 | height++; |
---|
281 | } |
---|
282 | return cf; |
---|
283 | } |
---|
284 | |
---|
285 | BOOLEAN naIsZero(number a, const coeffs cf) |
---|
286 | { |
---|
287 | naTest(a); |
---|
288 | return (a == NULL); |
---|
289 | } |
---|
290 | |
---|
291 | void naDelete(number * a, const coeffs cf) |
---|
292 | { |
---|
293 | if (*a == NULL) return; |
---|
294 | if (((poly)*a)==naMinpoly) { *a=NULL;return;} |
---|
295 | poly aAsPoly = (poly)(*a); |
---|
296 | p_Delete(&aAsPoly, naRing); |
---|
297 | *a = NULL; |
---|
298 | } |
---|
299 | |
---|
300 | BOOLEAN naEqual(number a, number b, const coeffs cf) |
---|
301 | { |
---|
302 | naTest(a); naTest(b); |
---|
303 | /// simple tests |
---|
304 | if (a == NULL) return (b == NULL); |
---|
305 | if (b == NULL) return (a == NULL); |
---|
306 | return p_EqualPolys((poly)a,(poly)b,naRing); |
---|
307 | } |
---|
308 | |
---|
309 | number naCopy(number a, const coeffs cf) |
---|
310 | { |
---|
311 | naTest(a); |
---|
312 | if (a == NULL) return NULL; |
---|
313 | if (((poly)a)==naMinpoly) return a; |
---|
314 | return (number)p_Copy((poly)a, naRing); |
---|
315 | } |
---|
316 | |
---|
317 | number naGetNumerator(number &a, const coeffs cf) |
---|
318 | { |
---|
319 | return naCopy(a, cf); |
---|
320 | } |
---|
321 | |
---|
322 | number naGetDenom(number &a, const coeffs cf) |
---|
323 | { |
---|
324 | naTest(a); |
---|
325 | return naInit(1, cf); |
---|
326 | } |
---|
327 | |
---|
328 | BOOLEAN naIsOne(number a, const coeffs cf) |
---|
329 | { |
---|
330 | naTest(a); |
---|
331 | poly aAsPoly = (poly)a; |
---|
332 | if ((a==NULL) || (!p_IsConstant(aAsPoly, naRing))) return FALSE; |
---|
333 | return n_IsOne(p_GetCoeff(aAsPoly, naRing), naCoeffs); |
---|
334 | } |
---|
335 | |
---|
336 | BOOLEAN naIsMOne(number a, const coeffs cf) |
---|
337 | { |
---|
338 | naTest(a); |
---|
339 | poly aAsPoly = (poly)a; |
---|
340 | if ((a==NULL) || (!p_IsConstant(aAsPoly, naRing))) return FALSE; |
---|
341 | return n_IsMOne(p_GetCoeff(aAsPoly, naRing), naCoeffs); |
---|
342 | } |
---|
343 | |
---|
344 | /// this is in-place, modifies a |
---|
345 | number naNeg(number a, const coeffs cf) |
---|
346 | { |
---|
347 | naTest(a); |
---|
348 | if (a != NULL) a = (number)p_Neg((poly)a, naRing); |
---|
349 | return a; |
---|
350 | } |
---|
351 | |
---|
352 | number naInit_bigint(number longratBigIntNumber, const coeffs src, const coeffs cf) |
---|
353 | { |
---|
354 | assume( cf != NULL ); |
---|
355 | |
---|
356 | const ring A = cf->extRing; |
---|
357 | |
---|
358 | assume( A != NULL ); |
---|
359 | |
---|
360 | const coeffs C = A->cf; |
---|
361 | |
---|
362 | assume( C != NULL ); |
---|
363 | |
---|
364 | number n = n_Init_bigint(longratBigIntNumber, src, C); |
---|
365 | |
---|
366 | if ( n_IsZero(n, C) ) |
---|
367 | { |
---|
368 | n_Delete(&n, C); |
---|
369 | return NULL; |
---|
370 | } |
---|
371 | |
---|
372 | return (number)p_NSet(n, A); |
---|
373 | } |
---|
374 | |
---|
375 | |
---|
376 | |
---|
377 | number naInit(long i, const coeffs cf) |
---|
378 | { |
---|
379 | if (i == 0) return NULL; |
---|
380 | else return (number)p_ISet(i, naRing); |
---|
381 | } |
---|
382 | |
---|
383 | int naInt(number &a, const coeffs cf) |
---|
384 | { |
---|
385 | naTest(a); |
---|
386 | poly aAsPoly = (poly)a; |
---|
387 | if(aAsPoly == NULL) |
---|
388 | return 0; |
---|
389 | if (!p_IsConstant(aAsPoly, naRing)) |
---|
390 | return 0; |
---|
391 | assume( aAsPoly != NULL ); |
---|
392 | return n_Int(p_GetCoeff(aAsPoly, naRing), naCoeffs); |
---|
393 | } |
---|
394 | |
---|
395 | /* TRUE iff (a != 0 and (b == 0 or deg(a) > deg(b) or (deg(a)==deg(b) && lc(a)>lc(b))) */ |
---|
396 | BOOLEAN naGreater(number a, number b, const coeffs cf) |
---|
397 | { |
---|
398 | naTest(a); naTest(b); |
---|
399 | if (naIsZero(a, cf)) |
---|
400 | { |
---|
401 | if (naIsZero(b, cf)) return FALSE; |
---|
402 | return !n_GreaterZero(pGetCoeff((poly)b),cf); |
---|
403 | } |
---|
404 | if (naIsZero(b, cf)) |
---|
405 | { |
---|
406 | return n_GreaterZero(pGetCoeff((poly)a),cf); |
---|
407 | } |
---|
408 | int aDeg = p_Totaldegree((poly)a, naRing); |
---|
409 | int bDeg = p_Totaldegree((poly)b, naRing); |
---|
410 | if (aDeg>bDeg) return TRUE; |
---|
411 | if (aDeg<bDeg) return FALSE; |
---|
412 | return n_Greater(pGetCoeff((poly)a),pGetCoeff((poly)b),naCoeffs); |
---|
413 | } |
---|
414 | |
---|
415 | /* TRUE iff a != 0 and (LC(a) > 0 or deg(a) > 0) */ |
---|
416 | BOOLEAN naGreaterZero(number a, const coeffs cf) |
---|
417 | { |
---|
418 | naTest(a); |
---|
419 | if (a == NULL) return FALSE; |
---|
420 | if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE; |
---|
421 | if (p_Totaldegree((poly)a, naRing) > 0) return TRUE; |
---|
422 | return FALSE; |
---|
423 | } |
---|
424 | |
---|
425 | void naCoeffWrite(const coeffs cf, BOOLEAN details) |
---|
426 | { |
---|
427 | assume( cf != NULL ); |
---|
428 | |
---|
429 | const ring A = cf->extRing; |
---|
430 | |
---|
431 | assume( A != NULL ); |
---|
432 | assume( A->cf != NULL ); |
---|
433 | |
---|
434 | n_CoeffWrite(A->cf, details); |
---|
435 | |
---|
436 | // rWrite(A); |
---|
437 | |
---|
438 | const int P = rVar(A); |
---|
439 | assume( P > 0 ); |
---|
440 | |
---|
441 | Print("// %d parameter : ", P); |
---|
442 | |
---|
443 | for (int nop=0; nop < P; nop ++) |
---|
444 | Print("%s ", rRingVar(nop, A)); |
---|
445 | |
---|
446 | PrintLn(); |
---|
447 | |
---|
448 | const ideal I = A->qideal; |
---|
449 | |
---|
450 | assume( I != NULL ); |
---|
451 | assume( IDELEMS(I) == 1 ); |
---|
452 | |
---|
453 | |
---|
454 | if ( details ) |
---|
455 | { |
---|
456 | PrintS("// minpoly : ("); |
---|
457 | p_Write0( I->m[0], A); |
---|
458 | PrintS(")"); |
---|
459 | } |
---|
460 | else |
---|
461 | PrintS("// minpoly : ..."); |
---|
462 | |
---|
463 | PrintLn(); |
---|
464 | |
---|
465 | /* |
---|
466 | char *x = rRingVar(0, A); |
---|
467 | |
---|
468 | Print("// Coefficients live in the extension field K[%s]/<f(%s)>\n", x, x); |
---|
469 | Print("// with the minimal polynomial f(%s) = %s\n", x, |
---|
470 | p_String(A->qideal->m[0], A)); |
---|
471 | PrintS("// and K: "); |
---|
472 | */ |
---|
473 | } |
---|
474 | |
---|
475 | number naAdd(number a, number b, const coeffs cf) |
---|
476 | { |
---|
477 | naTest(a); naTest(b); |
---|
478 | if (a == NULL) return naCopy(b, cf); |
---|
479 | if (b == NULL) return naCopy(a, cf); |
---|
480 | poly aPlusB = p_Add_q(p_Copy((poly)a, naRing), |
---|
481 | p_Copy((poly)b, naRing), naRing); |
---|
482 | //definiteReduce(aPlusB, naMinpoly, cf); |
---|
483 | return (number)aPlusB; |
---|
484 | } |
---|
485 | |
---|
486 | number naSub(number a, number b, const coeffs cf) |
---|
487 | { |
---|
488 | naTest(a); naTest(b); |
---|
489 | if (b == NULL) return naCopy(a, cf); |
---|
490 | poly minusB = p_Neg(p_Copy((poly)b, naRing), naRing); |
---|
491 | if (a == NULL) return (number)minusB; |
---|
492 | poly aMinusB = p_Add_q(p_Copy((poly)a, naRing), minusB, naRing); |
---|
493 | //definiteReduce(aMinusB, naMinpoly, cf); |
---|
494 | return (number)aMinusB; |
---|
495 | } |
---|
496 | |
---|
497 | number naMult(number a, number b, const coeffs cf) |
---|
498 | { |
---|
499 | naTest(a); naTest(b); |
---|
500 | if ((a == NULL)||(b == NULL)) return NULL; |
---|
501 | poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing), |
---|
502 | p_Copy((poly)b, naRing), naRing); |
---|
503 | definiteReduce(aTimesB, naMinpoly, cf); |
---|
504 | p_Normalize(aTimesB,naRing); |
---|
505 | return (number)aTimesB; |
---|
506 | } |
---|
507 | |
---|
508 | number naDiv(number a, number b, const coeffs cf) |
---|
509 | { |
---|
510 | naTest(a); naTest(b); |
---|
511 | if (b == NULL) WerrorS(nDivBy0); |
---|
512 | if (a == NULL) return NULL; |
---|
513 | poly bInverse = (poly)naInvers(b, cf); |
---|
514 | if(bInverse != NULL) // b is non-zero divisor! |
---|
515 | { |
---|
516 | poly aDivB = p_Mult_q(p_Copy((poly)a, naRing), bInverse, naRing); |
---|
517 | definiteReduce(aDivB, naMinpoly, cf); |
---|
518 | p_Normalize(aDivB,naRing); |
---|
519 | return (number)aDivB; |
---|
520 | } |
---|
521 | return NULL; |
---|
522 | } |
---|
523 | |
---|
524 | /* 0^0 = 0; |
---|
525 | for |exp| <= 7 compute power by a simple multiplication loop; |
---|
526 | for |exp| >= 8 compute power along binary presentation of |exp|, e.g. |
---|
527 | p^13 = p^1 * p^4 * p^8, where we utilise that |
---|
528 | p^(2^(k+1)) = p^(2^k) * p^(2^k); |
---|
529 | intermediate reduction modulo the minimal polynomial is controlled by |
---|
530 | the in-place method heuristicReduce(poly, poly, coeffs); see there. |
---|
531 | */ |
---|
532 | void naPower(number a, int exp, number *b, const coeffs cf) |
---|
533 | { |
---|
534 | naTest(a); |
---|
535 | |
---|
536 | /* special cases first */ |
---|
537 | if (a == NULL) |
---|
538 | { |
---|
539 | if (exp >= 0) *b = NULL; |
---|
540 | else WerrorS(nDivBy0); |
---|
541 | return; |
---|
542 | } |
---|
543 | else if (exp == 0) { *b = naInit(1, cf); return; } |
---|
544 | else if (exp == 1) { *b = naCopy(a, cf); return; } |
---|
545 | else if (exp == -1) { *b = naInvers(a, cf); return; } |
---|
546 | |
---|
547 | int expAbs = exp; if (expAbs < 0) expAbs = -expAbs; |
---|
548 | |
---|
549 | /* now compute a^expAbs */ |
---|
550 | poly pow; poly aAsPoly = (poly)a; |
---|
551 | if (expAbs <= 7) |
---|
552 | { |
---|
553 | pow = p_Copy(aAsPoly, naRing); |
---|
554 | for (int i = 2; i <= expAbs; i++) |
---|
555 | { |
---|
556 | pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing); |
---|
557 | heuristicReduce(pow, naMinpoly, cf); |
---|
558 | } |
---|
559 | definiteReduce(pow, naMinpoly, cf); |
---|
560 | } |
---|
561 | else |
---|
562 | { |
---|
563 | pow = p_ISet(1, naRing); |
---|
564 | poly factor = p_Copy(aAsPoly, naRing); |
---|
565 | while (expAbs != 0) |
---|
566 | { |
---|
567 | if (expAbs & 1) |
---|
568 | { |
---|
569 | pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing); |
---|
570 | heuristicReduce(pow, naMinpoly, cf); |
---|
571 | } |
---|
572 | expAbs = expAbs / 2; |
---|
573 | if (expAbs != 0) |
---|
574 | { |
---|
575 | factor = p_Mult_q(factor, p_Copy(factor, naRing), naRing); |
---|
576 | heuristicReduce(factor, naMinpoly, cf); |
---|
577 | } |
---|
578 | } |
---|
579 | p_Delete(&factor, naRing); |
---|
580 | definiteReduce(pow, naMinpoly, cf); |
---|
581 | } |
---|
582 | |
---|
583 | /* invert if original exponent was negative */ |
---|
584 | number n = (number)pow; |
---|
585 | if (exp < 0) |
---|
586 | { |
---|
587 | number m = naInvers(n, cf); |
---|
588 | naDelete(&n, cf); |
---|
589 | n = m; |
---|
590 | } |
---|
591 | *b = n; |
---|
592 | } |
---|
593 | |
---|
594 | /* may reduce p modulo the reducer by calling definiteReduce; |
---|
595 | the decision is made based on the following heuristic |
---|
596 | (which should also only be changed here in this method): |
---|
597 | if (deg(p) > 10*deg(reducer) then perform reduction; |
---|
598 | modifies p */ |
---|
599 | void heuristicReduce(poly &p, poly reducer, const coeffs cf) |
---|
600 | { |
---|
601 | #ifdef LDEBUG |
---|
602 | p_Test((poly)p, naRing); |
---|
603 | p_Test((poly)reducer, naRing); |
---|
604 | #endif |
---|
605 | if (p_Totaldegree(p, naRing) > 10 * p_Totaldegree(reducer, naRing)) |
---|
606 | definiteReduce(p, reducer, cf); |
---|
607 | } |
---|
608 | |
---|
609 | void naWriteLong(number &a, const coeffs cf) |
---|
610 | { |
---|
611 | naTest(a); |
---|
612 | if (a == NULL) |
---|
613 | StringAppendS("0"); |
---|
614 | else |
---|
615 | { |
---|
616 | poly aAsPoly = (poly)a; |
---|
617 | /* basically, just write aAsPoly using p_Write, |
---|
618 | but use brackets around the output, if a is not |
---|
619 | a constant living in naCoeffs = cf->extRing->cf */ |
---|
620 | BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing)); |
---|
621 | if (useBrackets) StringAppendS("("); |
---|
622 | p_String0Long(aAsPoly, naRing, naRing); |
---|
623 | if (useBrackets) StringAppendS(")"); |
---|
624 | } |
---|
625 | } |
---|
626 | |
---|
627 | void naWriteShort(number &a, const coeffs cf) |
---|
628 | { |
---|
629 | naTest(a); |
---|
630 | if (a == NULL) |
---|
631 | StringAppendS("0"); |
---|
632 | else |
---|
633 | { |
---|
634 | poly aAsPoly = (poly)a; |
---|
635 | /* basically, just write aAsPoly using p_Write, |
---|
636 | but use brackets around the output, if a is not |
---|
637 | a constant living in naCoeffs = cf->extRing->cf */ |
---|
638 | BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing)); |
---|
639 | if (useBrackets) StringAppendS("("); |
---|
640 | p_String0Short(aAsPoly, naRing, naRing); |
---|
641 | if (useBrackets) StringAppendS(")"); |
---|
642 | } |
---|
643 | } |
---|
644 | |
---|
645 | const char * naRead(const char *s, number *a, const coeffs cf) |
---|
646 | { |
---|
647 | poly aAsPoly; |
---|
648 | const char * result = p_Read(s, aAsPoly, naRing); |
---|
649 | if (aAsPoly!=NULL) definiteReduce(aAsPoly, naMinpoly, cf); |
---|
650 | *a = (number)aAsPoly; |
---|
651 | return result; |
---|
652 | } |
---|
653 | |
---|
654 | #if 0 |
---|
655 | /* implemented by the rule lcm(a, b) = a * b / gcd(a, b) */ |
---|
656 | number naLcm(number a, number b, const coeffs cf) |
---|
657 | { |
---|
658 | naTest(a); naTest(b); |
---|
659 | if (a == NULL) return NULL; |
---|
660 | if (b == NULL) return NULL; |
---|
661 | number theProduct = (number)p_Mult_q(p_Copy((poly)a, naRing), |
---|
662 | p_Copy((poly)b, naRing), naRing); |
---|
663 | /* note that theProduct needs not be reduced w.r.t. naMinpoly; |
---|
664 | but the final division will take care of the necessary reduction */ |
---|
665 | number theGcd = naGcd(a, b, cf); |
---|
666 | return naDiv(theProduct, theGcd, cf); |
---|
667 | } |
---|
668 | #endif |
---|
669 | number napLcm(number b, const coeffs cf) |
---|
670 | { |
---|
671 | number h=n_Init(1,naRing->cf); |
---|
672 | poly bb=(poly)b; |
---|
673 | number d; |
---|
674 | while(bb!=NULL) |
---|
675 | { |
---|
676 | d=n_Lcm(h,pGetCoeff(bb), naRing->cf); |
---|
677 | n_Delete(&h,naRing->cf); |
---|
678 | h=d; |
---|
679 | pIter(bb); |
---|
680 | } |
---|
681 | return h; |
---|
682 | } |
---|
683 | number naLcmContent(number a, number b, const coeffs cf) |
---|
684 | { |
---|
685 | if (nCoeff_is_Zp(naRing->cf)) return naCopy(a,cf); |
---|
686 | #if 0 |
---|
687 | else { |
---|
688 | number g = ndGcd(a, b, cf); |
---|
689 | return g; |
---|
690 | } |
---|
691 | #else |
---|
692 | { |
---|
693 | a=(number)p_Copy((poly)a,naRing); |
---|
694 | number t=napLcm(b,cf); |
---|
695 | if(!n_IsOne(t,naRing->cf)) |
---|
696 | { |
---|
697 | number bt, rr; |
---|
698 | poly xx=(poly)a; |
---|
699 | while (xx!=NULL) |
---|
700 | { |
---|
701 | bt = n_SubringGcd(t, pGetCoeff(xx), naRing->cf); |
---|
702 | rr = n_Mult(t, pGetCoeff(xx), naRing->cf); |
---|
703 | n_Delete(&pGetCoeff(xx),naRing->cf); |
---|
704 | pGetCoeff(xx) = n_Div(rr, bt, naRing->cf); |
---|
705 | n_Normalize(pGetCoeff(xx),naRing->cf); |
---|
706 | n_Delete(&bt,naRing->cf); |
---|
707 | n_Delete(&rr,naRing->cf); |
---|
708 | pIter(xx); |
---|
709 | } |
---|
710 | } |
---|
711 | n_Delete(&t,naRing->cf); |
---|
712 | return (number) a; |
---|
713 | } |
---|
714 | #endif |
---|
715 | } |
---|
716 | |
---|
717 | /* expects *param to be castable to AlgExtInfo */ |
---|
718 | static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param) |
---|
719 | { |
---|
720 | if (n_algExt != n) return FALSE; |
---|
721 | AlgExtInfo *e = (AlgExtInfo *)param; |
---|
722 | /* for extension coefficient fields we expect the underlying |
---|
723 | polynomial rings to be IDENTICAL, i.e. the SAME OBJECT; |
---|
724 | this expectation is based on the assumption that we have properly |
---|
725 | registered cf and perform reference counting rather than creating |
---|
726 | multiple copies of the same coefficient field/domain/ring */ |
---|
727 | if (naRing == e->r) |
---|
728 | return TRUE; |
---|
729 | /* (Note that then also the minimal ideals will necessarily be |
---|
730 | the same, as they are attached to the ring.) */ |
---|
731 | |
---|
732 | // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)... |
---|
733 | if( rEqual(naRing, e->r, TRUE) ) // also checks the equality of qideals |
---|
734 | { |
---|
735 | const ideal mi = naRing->qideal; |
---|
736 | assume( IDELEMS(mi) == 1 ); |
---|
737 | const ideal ii = e->r->qideal; |
---|
738 | assume( IDELEMS(ii) == 1 ); |
---|
739 | |
---|
740 | // TODO: the following should be extended for 2 *equal* rings... |
---|
741 | assume( p_EqualPolys(mi->m[0], ii->m[0], naRing, e->r) ); |
---|
742 | |
---|
743 | rDelete(e->r); |
---|
744 | |
---|
745 | return TRUE; |
---|
746 | } |
---|
747 | |
---|
748 | return FALSE; |
---|
749 | |
---|
750 | } |
---|
751 | |
---|
752 | int naSize(number a, const coeffs cf) |
---|
753 | { |
---|
754 | if (a == NULL) return -1; |
---|
755 | /* this has been taken from the old implementation of field extensions, |
---|
756 | where we computed the sum of the degree and the number of terms in |
---|
757 | (poly)a; so we leave it at that, for the time being; |
---|
758 | maybe, the number of terms alone is a better measure? */ |
---|
759 | poly aAsPoly = (poly)a; |
---|
760 | int theDegree = 0; int noOfTerms = 0; |
---|
761 | while (aAsPoly != NULL) |
---|
762 | { |
---|
763 | noOfTerms++; |
---|
764 | int d = p_GetExp(aAsPoly, 1, naRing); |
---|
765 | if (d > theDegree) theDegree = d; |
---|
766 | pIter(aAsPoly); |
---|
767 | } |
---|
768 | return theDegree + noOfTerms; |
---|
769 | } |
---|
770 | |
---|
771 | /* performs polynomial division and overrides p by the remainder |
---|
772 | of division of p by the reducer; |
---|
773 | modifies p */ |
---|
774 | void definiteReduce(poly &p, poly reducer, const coeffs cf) |
---|
775 | { |
---|
776 | #ifdef LDEBUG |
---|
777 | p_Test((poly)p, naRing); |
---|
778 | p_Test((poly)reducer, naRing); |
---|
779 | #endif |
---|
780 | if ((p!=NULL) && (p_GetExp(p,1,naRing)>=p_GetExp(reducer,1,naRing))) |
---|
781 | { |
---|
782 | p_PolyDiv(p, reducer, FALSE, naRing); |
---|
783 | } |
---|
784 | } |
---|
785 | |
---|
786 | void naNormalize(number &a, const coeffs cf) |
---|
787 | { |
---|
788 | poly aa=(poly)a; |
---|
789 | if (aa!=naMinpoly) |
---|
790 | definiteReduce(aa,naMinpoly,cf); |
---|
791 | a=(number)aa; |
---|
792 | } |
---|
793 | |
---|
794 | number naConvFactoryNSingN( const CanonicalForm n, const coeffs cf) |
---|
795 | { |
---|
796 | if (n.isZero()) return NULL; |
---|
797 | poly p=convFactoryPSingP(n,naRing); |
---|
798 | return (number)p; |
---|
799 | } |
---|
800 | CanonicalForm naConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf ) |
---|
801 | { |
---|
802 | naTest(n); |
---|
803 | if (n==NULL) return CanonicalForm(0); |
---|
804 | |
---|
805 | return convSingPFactoryP((poly)n,naRing); |
---|
806 | } |
---|
807 | |
---|
808 | /* IMPORTANT NOTE: Since an algebraic field extension is again a field, |
---|
809 | the gcd of two elements is not very interesting. (It |
---|
810 | is actually any unit in the field, i.e., any non- |
---|
811 | zero element.) Note that the below method does not operate |
---|
812 | in this strong sense but rather computes the gcd of |
---|
813 | two given elements in the underlying polynomial ring. */ |
---|
814 | number naGcd(number a, number b, const coeffs cf) |
---|
815 | { |
---|
816 | if (a==NULL) return naCopy(b,cf); |
---|
817 | if (b==NULL) return naCopy(a,cf); |
---|
818 | |
---|
819 | poly ax=(poly)a; |
---|
820 | poly bx=(poly)b; |
---|
821 | if (pNext(ax)!=NULL) |
---|
822 | return (number)p_Copy(ax, naRing); |
---|
823 | else |
---|
824 | { |
---|
825 | if(nCoeff_is_Zp(naRing->cf)) |
---|
826 | return naInit(1,cf); |
---|
827 | else |
---|
828 | { |
---|
829 | number x = n_Copy(pGetCoeff((poly)a),naRing->cf); |
---|
830 | if (n_IsOne(x,naRing->cf)) |
---|
831 | return (number)p_NSet(x,naRing); |
---|
832 | while (pNext(ax)!=NULL) |
---|
833 | { |
---|
834 | pIter(ax); |
---|
835 | number y = n_SubringGcd(x, pGetCoeff(ax), naRing->cf); |
---|
836 | n_Delete(&x,naRing->cf); |
---|
837 | x = y; |
---|
838 | if (n_IsOne(x,naRing->cf)) |
---|
839 | return (number)p_NSet(x,naRing); |
---|
840 | } |
---|
841 | do |
---|
842 | { |
---|
843 | number y = n_SubringGcd(x, pGetCoeff(bx), naRing->cf); |
---|
844 | n_Delete(&x,naRing->cf); |
---|
845 | x = y; |
---|
846 | if (n_IsOne(x,naRing->cf)) |
---|
847 | return (number)p_NSet(x,naRing); |
---|
848 | pIter(bx); |
---|
849 | } |
---|
850 | while (bx!=NULL); |
---|
851 | return (number)p_NSet(x,naRing); |
---|
852 | } |
---|
853 | } |
---|
854 | #if 0 |
---|
855 | naTest(a); naTest(b); |
---|
856 | const ring R = naRing; |
---|
857 | return (number) singclap_gcd(p_Copy((poly)a, R), p_Copy((poly)b, R), R); |
---|
858 | #endif |
---|
859 | // return (number)p_Gcd((poly)a, (poly)b, naRing); |
---|
860 | } |
---|
861 | |
---|
862 | number naInvers(number a, const coeffs cf) |
---|
863 | { |
---|
864 | naTest(a); |
---|
865 | if (a == NULL) WerrorS(nDivBy0); |
---|
866 | |
---|
867 | poly aFactor = NULL; poly mFactor = NULL; poly theGcd = NULL; |
---|
868 | // singclap_extgcd! |
---|
869 | const BOOLEAN ret = singclap_extgcd ((poly)a, naMinpoly, theGcd, aFactor, mFactor, naRing); |
---|
870 | |
---|
871 | assume( !ret ); |
---|
872 | |
---|
873 | // if( ret ) theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing); |
---|
874 | |
---|
875 | naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor); |
---|
876 | p_Delete(&mFactor, naRing); |
---|
877 | |
---|
878 | // /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */ |
---|
879 | // assume(naIsOne((number)theGcd, cf)); |
---|
880 | |
---|
881 | if( !naIsOne((number)theGcd, cf) ) |
---|
882 | { |
---|
883 | WerrorS("zero divisor found - your minpoly is not irreducible"); |
---|
884 | p_Delete(&aFactor, naRing); aFactor = NULL; |
---|
885 | } |
---|
886 | p_Delete(&theGcd, naRing); |
---|
887 | |
---|
888 | return (number)(aFactor); |
---|
889 | } |
---|
890 | |
---|
891 | /* assumes that src = Q, dst = Q(a) */ |
---|
892 | number naMap00(number a, const coeffs src, const coeffs dst) |
---|
893 | { |
---|
894 | if (n_IsZero(a, src)) return NULL; |
---|
895 | assume(src == dst->extRing->cf); |
---|
896 | poly result = p_One(dst->extRing); |
---|
897 | p_SetCoeff(result, n_Copy(a, src), dst->extRing); |
---|
898 | return (number)result; |
---|
899 | } |
---|
900 | |
---|
901 | /* assumes that src = Z/p, dst = Q(a) */ |
---|
902 | number naMapP0(number a, const coeffs src, const coeffs dst) |
---|
903 | { |
---|
904 | if (n_IsZero(a, src)) return NULL; |
---|
905 | /* mapping via intermediate int: */ |
---|
906 | int n = n_Int(a, src); |
---|
907 | number q = n_Init(n, dst->extRing->cf); |
---|
908 | poly result = p_One(dst->extRing); |
---|
909 | p_SetCoeff(result, q, dst->extRing); |
---|
910 | return (number)result; |
---|
911 | } |
---|
912 | |
---|
913 | #if 0 |
---|
914 | /* assumes that either src = Q(a), dst = Q(a), or |
---|
915 | src = Z/p(a), dst = Z/p(a) */ |
---|
916 | number naCopyMap(number a, const coeffs src, const coeffs dst) |
---|
917 | { |
---|
918 | return naCopy(a, dst); |
---|
919 | } |
---|
920 | #endif |
---|
921 | |
---|
922 | number naCopyTrans2AlgExt(number a, const coeffs src, const coeffs dst) |
---|
923 | { |
---|
924 | assume (nCoeff_is_transExt (src)); |
---|
925 | assume (nCoeff_is_algExt (dst)); |
---|
926 | fraction fa=(fraction)a; |
---|
927 | poly p, q; |
---|
928 | if (rSamePolyRep(src->extRing, dst->extRing)) |
---|
929 | { |
---|
930 | p = p_Copy(NUM(fa),src->extRing); |
---|
931 | if (!DENIS1(fa)) |
---|
932 | { |
---|
933 | q = p_Copy(DEN(fa),src->extRing); |
---|
934 | assume (q != NULL); |
---|
935 | } |
---|
936 | } |
---|
937 | else |
---|
938 | { |
---|
939 | assume ((strcmp(rRingVar(0,src->extRing),rRingVar(0,dst->extRing))==0) && (rVar (src->extRing) == rVar (dst->extRing))); |
---|
940 | |
---|
941 | nMapFunc nMap= n_SetMap (src->extRing->cf, dst->extRing->cf); |
---|
942 | |
---|
943 | assume (nMap != NULL); |
---|
944 | p= p_PermPoly (NUM (fa), NULL, src->extRing, dst->extRing,nMap, NULL,rVar (src->extRing)); |
---|
945 | if (!DENIS1(fa)) |
---|
946 | { |
---|
947 | q= p_PermPoly (DEN (fa), NULL, src->extRing, dst->extRing,nMap, NULL,rVar (src->extRing)); |
---|
948 | assume (q != NULL); |
---|
949 | } |
---|
950 | } |
---|
951 | definiteReduce(p, dst->extRing->qideal->m[0], dst); |
---|
952 | assume (p_Test (p, dst->extRing)); |
---|
953 | if (!DENIS1(fa)) |
---|
954 | { |
---|
955 | definiteReduce(q, dst->extRing->qideal->m[0], dst); |
---|
956 | assume (p_Test (q, dst->extRing)); |
---|
957 | if (q != NULL) |
---|
958 | { |
---|
959 | number t= naDiv ((number)p,(number)q, dst); |
---|
960 | p_Delete (&p, dst->extRing); |
---|
961 | p_Delete (&q, dst->extRing); |
---|
962 | return t; |
---|
963 | } |
---|
964 | WerrorS ("mapping denominator to zero"); |
---|
965 | } |
---|
966 | return (number) p; |
---|
967 | } |
---|
968 | |
---|
969 | /* assumes that src = Q, dst = Z/p(a) */ |
---|
970 | number naMap0P(number a, const coeffs src, const coeffs dst) |
---|
971 | { |
---|
972 | if (n_IsZero(a, src)) return NULL; |
---|
973 | // int p = rChar(dst->extRing); |
---|
974 | |
---|
975 | number q = nlModP(a, src, dst->extRing->cf); |
---|
976 | |
---|
977 | poly result = p_NSet(q, dst->extRing); |
---|
978 | |
---|
979 | return (number)result; |
---|
980 | } |
---|
981 | |
---|
982 | /* assumes that src = Z/p, dst = Z/p(a) */ |
---|
983 | number naMapPP(number a, const coeffs src, const coeffs dst) |
---|
984 | { |
---|
985 | if (n_IsZero(a, src)) return NULL; |
---|
986 | assume(src == dst->extRing->cf); |
---|
987 | poly result = p_One(dst->extRing); |
---|
988 | p_SetCoeff(result, n_Copy(a, src), dst->extRing); |
---|
989 | return (number)result; |
---|
990 | } |
---|
991 | |
---|
992 | /* assumes that src = Z/u, dst = Z/p(a), where u != p */ |
---|
993 | number naMapUP(number a, const coeffs src, const coeffs dst) |
---|
994 | { |
---|
995 | if (n_IsZero(a, src)) return NULL; |
---|
996 | /* mapping via intermediate int: */ |
---|
997 | int n = n_Int(a, src); |
---|
998 | number q = n_Init(n, dst->extRing->cf); |
---|
999 | poly result = p_One(dst->extRing); |
---|
1000 | p_SetCoeff(result, q, dst->extRing); |
---|
1001 | return (number)result; |
---|
1002 | } |
---|
1003 | |
---|
1004 | number naGenMap(number a, const coeffs cf, const coeffs dst) |
---|
1005 | { |
---|
1006 | if (a==NULL) return NULL; |
---|
1007 | |
---|
1008 | const ring rSrc = cf->extRing; |
---|
1009 | const ring rDst = dst->extRing; |
---|
1010 | |
---|
1011 | const nMapFunc nMap=n_SetMap(rSrc->cf,rDst->cf); |
---|
1012 | poly f = (poly)a; |
---|
1013 | poly g = prMapR(f, nMap, rSrc, rDst); |
---|
1014 | |
---|
1015 | assume(n_Test((number)g, dst)); |
---|
1016 | return (number)g; |
---|
1017 | } |
---|
1018 | |
---|
1019 | number naGenTrans2AlgExt(number a, const coeffs cf, const coeffs dst) |
---|
1020 | { |
---|
1021 | if (a==NULL) return NULL; |
---|
1022 | |
---|
1023 | const ring rSrc = cf->extRing; |
---|
1024 | const ring rDst = dst->extRing; |
---|
1025 | |
---|
1026 | const nMapFunc nMap=n_SetMap(rSrc->cf,rDst->cf); |
---|
1027 | fraction f = (fraction)a; |
---|
1028 | poly g = prMapR(NUM(f), nMap, rSrc, rDst); |
---|
1029 | |
---|
1030 | number result=NULL; |
---|
1031 | poly h = NULL; |
---|
1032 | |
---|
1033 | if (!DENIS1(f)) |
---|
1034 | h = prMapR(DEN(f), nMap, rSrc, rDst); |
---|
1035 | |
---|
1036 | if (h!=NULL) |
---|
1037 | { |
---|
1038 | result=naDiv((number)g,(number)h,dst); |
---|
1039 | p_Delete(&g,dst->extRing); |
---|
1040 | p_Delete(&h,dst->extRing); |
---|
1041 | } |
---|
1042 | else |
---|
1043 | result=(number)g; |
---|
1044 | |
---|
1045 | assume(n_Test((number)result, dst)); |
---|
1046 | return (number)result; |
---|
1047 | } |
---|
1048 | |
---|
1049 | nMapFunc naSetMap(const coeffs src, const coeffs dst) |
---|
1050 | { |
---|
1051 | /* dst is expected to be an algebraic field extension */ |
---|
1052 | assume(getCoeffType(dst) == ID); |
---|
1053 | |
---|
1054 | if( src == dst ) return ndCopyMap; |
---|
1055 | |
---|
1056 | int h = 0; /* the height of the extension tower given by dst */ |
---|
1057 | coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */ |
---|
1058 | coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */ |
---|
1059 | |
---|
1060 | /* for the time being, we only provide maps if h = 1 or 0 */ |
---|
1061 | if (h==0) |
---|
1062 | { |
---|
1063 | if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst)) |
---|
1064 | return naMap00; /// Q --> Q(a) |
---|
1065 | if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst)) |
---|
1066 | return naMapP0; /// Z/p --> Q(a) |
---|
1067 | if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst)) |
---|
1068 | return naMap0P; /// Q --> Z/p(a) |
---|
1069 | if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst)) |
---|
1070 | { |
---|
1071 | if (src->ch == dst->ch) return naMapPP; /// Z/p --> Z/p(a) |
---|
1072 | else return naMapUP; /// Z/u --> Z/p(a) |
---|
1073 | } |
---|
1074 | } |
---|
1075 | if (h != 1) return NULL; |
---|
1076 | if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL; |
---|
1077 | if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL; |
---|
1078 | |
---|
1079 | nMapFunc nMap=n_SetMap(src->extRing->cf,dst->extRing->cf); |
---|
1080 | if (rSamePolyRep(src->extRing, dst->extRing) && (strcmp(rRingVar(0, src->extRing), rRingVar(0, dst->extRing)) == 0)) |
---|
1081 | { |
---|
1082 | if (src->type==n_algExt) |
---|
1083 | return ndCopyMap; // naCopyMap; /// K(a) --> K(a) |
---|
1084 | else |
---|
1085 | return naCopyTrans2AlgExt; |
---|
1086 | } |
---|
1087 | else if ((nMap!=NULL) && (strcmp(rRingVar(0,src->extRing),rRingVar(0,dst->extRing))==0) && (rVar (src->extRing) == rVar (dst->extRing))) |
---|
1088 | { |
---|
1089 | if (src->type==n_algExt) |
---|
1090 | return naGenMap; // naCopyMap; /// K(a) --> K'(a) |
---|
1091 | else |
---|
1092 | return naGenTrans2AlgExt; |
---|
1093 | } |
---|
1094 | |
---|
1095 | return NULL; /// default |
---|
1096 | } |
---|
1097 | |
---|
1098 | static int naParDeg(number a, const coeffs cf) |
---|
1099 | { |
---|
1100 | if (a == NULL) return -1; |
---|
1101 | poly aa=(poly)a; |
---|
1102 | return cf->extRing->pFDeg(aa,cf->extRing); |
---|
1103 | } |
---|
1104 | |
---|
1105 | /// return the specified parameter as a number in the given alg. field |
---|
1106 | static number naParameter(const int iParameter, const coeffs cf) |
---|
1107 | { |
---|
1108 | assume(getCoeffType(cf) == ID); |
---|
1109 | |
---|
1110 | const ring R = cf->extRing; |
---|
1111 | assume( R != NULL ); |
---|
1112 | assume( 0 < iParameter && iParameter <= rVar(R) ); |
---|
1113 | |
---|
1114 | poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R); |
---|
1115 | |
---|
1116 | return (number) p; |
---|
1117 | } |
---|
1118 | |
---|
1119 | |
---|
1120 | /// if m == var(i)/1 => return i, |
---|
1121 | int naIsParam(number m, const coeffs cf) |
---|
1122 | { |
---|
1123 | assume(getCoeffType(cf) == ID); |
---|
1124 | |
---|
1125 | const ring R = cf->extRing; |
---|
1126 | assume( R != NULL ); |
---|
1127 | |
---|
1128 | return p_Var( (poly)m, R ); |
---|
1129 | } |
---|
1130 | |
---|
1131 | |
---|
1132 | static void naClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf) |
---|
1133 | { |
---|
1134 | assume(cf != NULL); |
---|
1135 | assume(getCoeffType(cf) == ID); |
---|
1136 | assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) ! |
---|
1137 | |
---|
1138 | const ring R = cf->extRing; |
---|
1139 | assume(R != NULL); |
---|
1140 | const coeffs Q = R->cf; |
---|
1141 | assume(Q != NULL); |
---|
1142 | assume(nCoeff_is_Q(Q)); |
---|
1143 | |
---|
1144 | numberCollectionEnumerator.Reset(); |
---|
1145 | |
---|
1146 | if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial? |
---|
1147 | { |
---|
1148 | c = n_Init(1, cf); |
---|
1149 | return; |
---|
1150 | } |
---|
1151 | |
---|
1152 | naTest(numberCollectionEnumerator.Current()); |
---|
1153 | |
---|
1154 | // part 1, find a small candidate for gcd |
---|
1155 | int s1; int s=2147483647; // max. int |
---|
1156 | |
---|
1157 | const BOOLEAN lc_is_pos=naGreaterZero(numberCollectionEnumerator.Current(),cf); |
---|
1158 | |
---|
1159 | int normalcount = 0; |
---|
1160 | |
---|
1161 | poly cand1, cand; |
---|
1162 | |
---|
1163 | do |
---|
1164 | { |
---|
1165 | number& n = numberCollectionEnumerator.Current(); |
---|
1166 | naNormalize(n, cf); ++normalcount; |
---|
1167 | |
---|
1168 | naTest(n); |
---|
1169 | |
---|
1170 | cand1 = (poly)n; |
---|
1171 | |
---|
1172 | s1 = p_Deg(cand1, R); // naSize? |
---|
1173 | if (s>s1) |
---|
1174 | { |
---|
1175 | cand = cand1; |
---|
1176 | s = s1; |
---|
1177 | } |
---|
1178 | } while (numberCollectionEnumerator.MoveNext() ); |
---|
1179 | |
---|
1180 | // assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer! |
---|
1181 | |
---|
1182 | cand = p_Copy(cand, R); |
---|
1183 | // part 2: compute gcd(cand,all coeffs) |
---|
1184 | |
---|
1185 | numberCollectionEnumerator.Reset(); |
---|
1186 | |
---|
1187 | int length = 0; |
---|
1188 | while (numberCollectionEnumerator.MoveNext() ) |
---|
1189 | { |
---|
1190 | number& n = numberCollectionEnumerator.Current(); |
---|
1191 | ++length; |
---|
1192 | |
---|
1193 | if( (--normalcount) <= 0) |
---|
1194 | naNormalize(n, cf); |
---|
1195 | |
---|
1196 | naTest(n); |
---|
1197 | |
---|
1198 | // p_InpGcd(cand, (poly)n, R); |
---|
1199 | |
---|
1200 | cand = singclap_gcd(cand, p_Copy((poly)n, R), R); |
---|
1201 | |
---|
1202 | // cand1 = p_Gcd(cand,(poly)n, R); p_Delete(&cand, R); cand = cand1; |
---|
1203 | |
---|
1204 | assume( naGreaterZero((number)cand, cf) ); // ??? |
---|
1205 | /* |
---|
1206 | if(p_IsConstant(cand,R)) |
---|
1207 | { |
---|
1208 | c = cand; |
---|
1209 | |
---|
1210 | if(!lc_is_pos) |
---|
1211 | { |
---|
1212 | // make the leading coeff positive |
---|
1213 | c = nlNeg(c, cf); |
---|
1214 | numberCollectionEnumerator.Reset(); |
---|
1215 | |
---|
1216 | while (numberCollectionEnumerator.MoveNext() ) |
---|
1217 | { |
---|
1218 | number& nn = numberCollectionEnumerator.Current(); |
---|
1219 | nn = nlNeg(nn, cf); |
---|
1220 | } |
---|
1221 | } |
---|
1222 | return; |
---|
1223 | } |
---|
1224 | */ |
---|
1225 | |
---|
1226 | } |
---|
1227 | |
---|
1228 | |
---|
1229 | // part3: all coeffs = all coeffs / cand |
---|
1230 | if (!lc_is_pos) |
---|
1231 | cand = p_Neg(cand, R); |
---|
1232 | |
---|
1233 | c = (number)cand; naTest(c); |
---|
1234 | |
---|
1235 | poly cInverse = (poly)naInvers(c, cf); |
---|
1236 | assume(cInverse != NULL); // c is non-zero divisor!? |
---|
1237 | |
---|
1238 | |
---|
1239 | numberCollectionEnumerator.Reset(); |
---|
1240 | |
---|
1241 | |
---|
1242 | while (numberCollectionEnumerator.MoveNext() ) |
---|
1243 | { |
---|
1244 | number& n = numberCollectionEnumerator.Current(); |
---|
1245 | |
---|
1246 | assume( length > 0 ); |
---|
1247 | |
---|
1248 | if( --length > 0 ) |
---|
1249 | { |
---|
1250 | assume( cInverse != NULL ); |
---|
1251 | n = (number) p_Mult_q(p_Copy(cInverse, R), (poly)n, R); |
---|
1252 | } |
---|
1253 | else |
---|
1254 | { |
---|
1255 | n = (number) p_Mult_q(cInverse, (poly)n, R); |
---|
1256 | cInverse = NULL; |
---|
1257 | assume(length == 0); |
---|
1258 | } |
---|
1259 | |
---|
1260 | definiteReduce((poly &)n, naMinpoly, cf); |
---|
1261 | } |
---|
1262 | |
---|
1263 | assume(length == 0); |
---|
1264 | assume(cInverse == NULL); // p_Delete(&cInverse, R); |
---|
1265 | |
---|
1266 | // Quick and dirty fix for constant content clearing... !? |
---|
1267 | CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys! |
---|
1268 | |
---|
1269 | number cc; |
---|
1270 | |
---|
1271 | n_ClearContent(itr, cc, Q); // TODO: get rid of (-LC) normalization!? |
---|
1272 | |
---|
1273 | // over alg. ext. of Q // takes over the input number |
---|
1274 | c = (number) p_Mult_nn( (poly)c, cc, R); |
---|
1275 | // p_Mult_q(p_NSet(cc, R), , R); |
---|
1276 | |
---|
1277 | n_Delete(&cc, Q); |
---|
1278 | |
---|
1279 | // TODO: the above is not enough! need GCD's of polynomial coeffs...! |
---|
1280 | /* |
---|
1281 | // old and wrong part of p_Content |
---|
1282 | if (rField_is_Q_a(r) && !CLEARENUMERATORS) // should not be used anymore if CLEARENUMERATORS is 1 |
---|
1283 | { |
---|
1284 | // we only need special handling for alg. ext. |
---|
1285 | if (getCoeffType(r->cf)==n_algExt) |
---|
1286 | { |
---|
1287 | number hzz = n_Init(1, r->cf->extRing->cf); |
---|
1288 | p=ph; |
---|
1289 | while (p!=NULL) |
---|
1290 | { // each monom: coeff in Q_a |
---|
1291 | poly c_n_n=(poly)pGetCoeff(p); |
---|
1292 | poly c_n=c_n_n; |
---|
1293 | while (c_n!=NULL) |
---|
1294 | { // each monom: coeff in Q |
---|
1295 | d=n_Lcm(hzz,pGetCoeff(c_n),r->cf->extRing->cf); |
---|
1296 | n_Delete(&hzz,r->cf->extRing->cf); |
---|
1297 | hzz=d; |
---|
1298 | pIter(c_n); |
---|
1299 | } |
---|
1300 | pIter(p); |
---|
1301 | } |
---|
1302 | // hzz contains the 1/lcm of all denominators in c_n_n |
---|
1303 | h=n_Invers(hzz,r->cf->extRing->cf); |
---|
1304 | n_Delete(&hzz,r->cf->extRing->cf); |
---|
1305 | n_Normalize(h,r->cf->extRing->cf); |
---|
1306 | if(!n_IsOne(h,r->cf->extRing->cf)) |
---|
1307 | { |
---|
1308 | p=ph; |
---|
1309 | while (p!=NULL) |
---|
1310 | { // each monom: coeff in Q_a |
---|
1311 | poly c_n=(poly)pGetCoeff(p); |
---|
1312 | while (c_n!=NULL) |
---|
1313 | { // each monom: coeff in Q |
---|
1314 | d=n_Mult(h,pGetCoeff(c_n),r->cf->extRing->cf); |
---|
1315 | n_Normalize(d,r->cf->extRing->cf); |
---|
1316 | n_Delete(&pGetCoeff(c_n),r->cf->extRing->cf); |
---|
1317 | pGetCoeff(c_n)=d; |
---|
1318 | pIter(c_n); |
---|
1319 | } |
---|
1320 | pIter(p); |
---|
1321 | } |
---|
1322 | } |
---|
1323 | n_Delete(&h,r->cf->extRing->cf); |
---|
1324 | } |
---|
1325 | } |
---|
1326 | */ |
---|
1327 | |
---|
1328 | |
---|
1329 | // c = n_Init(1, cf); assume(FALSE); // TODO: NOT YET IMPLEMENTED!!! |
---|
1330 | } |
---|
1331 | |
---|
1332 | |
---|
1333 | static void naClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf) |
---|
1334 | { |
---|
1335 | assume(cf != NULL); |
---|
1336 | assume(getCoeffType(cf) == ID); |
---|
1337 | assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) ! |
---|
1338 | |
---|
1339 | assume(cf->extRing != NULL); |
---|
1340 | const coeffs Q = cf->extRing->cf; |
---|
1341 | assume(Q != NULL); |
---|
1342 | assume(nCoeff_is_Q(Q)); |
---|
1343 | number n; |
---|
1344 | CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys! |
---|
1345 | n_ClearDenominators(itr, n, Q); // this should probably be fine... |
---|
1346 | c = (number)p_NSet(n, cf->extRing); // over alg. ext. of Q // takes over the input number |
---|
1347 | } |
---|
1348 | |
---|
1349 | void naKillChar(coeffs cf) |
---|
1350 | { |
---|
1351 | if ((--cf->extRing->ref) == 0) |
---|
1352 | rDelete(cf->extRing); |
---|
1353 | } |
---|
1354 | |
---|
1355 | char* naCoeffString(const coeffs r) // currently also for tranext. |
---|
1356 | { |
---|
1357 | const char* const* p=n_ParameterNames(r); |
---|
1358 | int l=0; |
---|
1359 | int i; |
---|
1360 | for(i=0; i<n_NumberOfParameters(r);i++) |
---|
1361 | { |
---|
1362 | l+=(strlen(p[i])+1); |
---|
1363 | } |
---|
1364 | char *s=(char *)omAlloc(l+10+1); |
---|
1365 | s[0]='\0'; |
---|
1366 | snprintf(s,10+1,"%d",r->ch); /* Fp(a) or Q(a) */ |
---|
1367 | char tt[2]; |
---|
1368 | tt[0]=','; |
---|
1369 | tt[1]='\0'; |
---|
1370 | for(i=0; i<n_NumberOfParameters(r);i++) |
---|
1371 | { |
---|
1372 | strcat(s,tt); |
---|
1373 | strcat(s,p[i]); |
---|
1374 | } |
---|
1375 | return s; |
---|
1376 | } |
---|
1377 | |
---|
1378 | number naChineseRemainder(number *x, number *q,int rl, BOOLEAN sym,const coeffs cf) |
---|
1379 | { |
---|
1380 | poly *P=(poly*)omAlloc(rl*sizeof(poly*)); |
---|
1381 | number *X=(number *)omAlloc(rl*sizeof(number)); |
---|
1382 | int i; |
---|
1383 | for(i=0;i<rl;i++) P[i]=p_Copy((poly)(x[i]),cf->extRing); |
---|
1384 | poly result=p_ChineseRemainder(P,X,q,rl,cf->extRing); |
---|
1385 | omFreeSize(X,rl*sizeof(number)); |
---|
1386 | omFreeSize(P,rl*sizeof(poly*)); |
---|
1387 | return ((number)result); |
---|
1388 | } |
---|
1389 | |
---|
1390 | number naFarey(number p, number n, const coeffs cf) |
---|
1391 | { |
---|
1392 | // n is really a bigint |
---|
1393 | poly result=p_Farey(p_Copy((poly)p,cf->extRing),n,cf->extRing); |
---|
1394 | return ((number)result); |
---|
1395 | } |
---|
1396 | |
---|
1397 | |
---|
1398 | BOOLEAN naInitChar(coeffs cf, void * infoStruct) |
---|
1399 | { |
---|
1400 | assume( infoStruct != NULL ); |
---|
1401 | |
---|
1402 | AlgExtInfo *e = (AlgExtInfo *)infoStruct; |
---|
1403 | /// first check whether cf->extRing != NULL and delete old ring??? |
---|
1404 | |
---|
1405 | assume(e->r != NULL); // extRing; |
---|
1406 | assume(e->r->cf != NULL); // extRing->cf; |
---|
1407 | |
---|
1408 | assume((e->r->qideal != NULL) && // minideal has one |
---|
1409 | (IDELEMS(e->r->qideal) == 1) && // non-zero generator |
---|
1410 | (e->r->qideal->m[0] != NULL) ); // at m[0]; |
---|
1411 | |
---|
1412 | assume( cf != NULL ); |
---|
1413 | assume(getCoeffType(cf) == ID); // coeff type; |
---|
1414 | |
---|
1415 | e->r->ref ++; // increase the ref.counter for the ground poly. ring! |
---|
1416 | const ring R = e->r; // no copy! |
---|
1417 | cf->extRing = R; |
---|
1418 | |
---|
1419 | /* propagate characteristic up so that it becomes |
---|
1420 | directly accessible in cf: */ |
---|
1421 | cf->ch = R->cf->ch; |
---|
1422 | |
---|
1423 | cf->is_field=TRUE; |
---|
1424 | cf->is_domain=TRUE; |
---|
1425 | |
---|
1426 | #ifdef LDEBUG |
---|
1427 | p_Test((poly)naMinpoly, naRing); |
---|
1428 | #endif |
---|
1429 | |
---|
1430 | cf->cfCoeffString = naCoeffString; |
---|
1431 | |
---|
1432 | cf->cfGreaterZero = naGreaterZero; |
---|
1433 | cf->cfGreater = naGreater; |
---|
1434 | cf->cfEqual = naEqual; |
---|
1435 | cf->cfIsZero = naIsZero; |
---|
1436 | cf->cfIsOne = naIsOne; |
---|
1437 | cf->cfIsMOne = naIsMOne; |
---|
1438 | cf->cfInit = naInit; |
---|
1439 | cf->cfInit_bigint = naInit_bigint; |
---|
1440 | cf->cfFarey = naFarey; |
---|
1441 | cf->cfChineseRemainder= naChineseRemainder; |
---|
1442 | cf->cfInt = naInt; |
---|
1443 | cf->cfInpNeg = naNeg; |
---|
1444 | cf->cfAdd = naAdd; |
---|
1445 | cf->cfSub = naSub; |
---|
1446 | cf->cfMult = naMult; |
---|
1447 | cf->cfDiv = naDiv; |
---|
1448 | cf->cfExactDiv = naDiv; |
---|
1449 | cf->cfPower = naPower; |
---|
1450 | cf->cfCopy = naCopy; |
---|
1451 | |
---|
1452 | cf->cfWriteLong = naWriteLong; |
---|
1453 | |
---|
1454 | if( rCanShortOut(naRing) ) |
---|
1455 | cf->cfWriteShort = naWriteShort; |
---|
1456 | else |
---|
1457 | cf->cfWriteShort = naWriteLong; |
---|
1458 | |
---|
1459 | cf->cfRead = naRead; |
---|
1460 | cf->cfDelete = naDelete; |
---|
1461 | cf->cfSetMap = naSetMap; |
---|
1462 | cf->cfGetDenom = naGetDenom; |
---|
1463 | cf->cfGetNumerator = naGetNumerator; |
---|
1464 | cf->cfRePart = naCopy; |
---|
1465 | cf->cfCoeffWrite = naCoeffWrite; |
---|
1466 | cf->cfNormalize = naNormalize; |
---|
1467 | cf->cfKillChar = naKillChar; |
---|
1468 | #ifdef LDEBUG |
---|
1469 | cf->cfDBTest = naDBTest; |
---|
1470 | #endif |
---|
1471 | cf->cfGcd = naGcd; |
---|
1472 | cf->cfLcm = naLcmContent; |
---|
1473 | cf->cfSize = naSize; |
---|
1474 | cf->nCoeffIsEqual = naCoeffIsEqual; |
---|
1475 | cf->cfInvers = naInvers; |
---|
1476 | cf->convFactoryNSingN=naConvFactoryNSingN; |
---|
1477 | cf->convSingNFactoryN=naConvSingNFactoryN; |
---|
1478 | cf->cfParDeg = naParDeg; |
---|
1479 | |
---|
1480 | cf->iNumberOfParameters = rVar(R); |
---|
1481 | cf->pParameterNames = (const char**)R->names; |
---|
1482 | cf->cfParameter = naParameter; |
---|
1483 | cf->has_simple_Inverse= R->cf->has_simple_Inverse; |
---|
1484 | /* cf->has_simple_Alloc= FALSE; */ |
---|
1485 | |
---|
1486 | if( nCoeff_is_Q(R->cf) ) |
---|
1487 | { |
---|
1488 | cf->cfClearContent = naClearContent; |
---|
1489 | cf->cfClearDenominators = naClearDenominators; |
---|
1490 | } |
---|
1491 | |
---|
1492 | return FALSE; |
---|
1493 | } |
---|
1494 | |
---|
1495 | template class CRecursivePolyCoeffsEnumerator<NAConverter>; |
---|
1496 | |
---|
1497 | template class IAccessor<snumber*>; |
---|
1498 | |
---|
1499 | /* --------------------------------------------------------------------*/ |
---|
1500 | #if 0 |
---|
1501 | void npolyCoeffWrite(const coeffs cf, BOOLEAN details) |
---|
1502 | { |
---|
1503 | assume( cf != NULL ); |
---|
1504 | |
---|
1505 | const ring A = cf->extRing; |
---|
1506 | |
---|
1507 | assume( A != NULL ); |
---|
1508 | Print("// polynomial ring as coefficient ring :\n"); |
---|
1509 | rWrite(A); |
---|
1510 | PrintLn(); |
---|
1511 | } |
---|
1512 | number npolyMult(number a, number b, const coeffs cf) |
---|
1513 | { |
---|
1514 | naTest(a); naTest(b); |
---|
1515 | if ((a == NULL)||(b == NULL)) return NULL; |
---|
1516 | poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing), |
---|
1517 | p_Copy((poly)b, naRing), naRing); |
---|
1518 | return (number)aTimesB; |
---|
1519 | } |
---|
1520 | |
---|
1521 | void npolyPower(number a, int exp, number *b, const coeffs cf) |
---|
1522 | { |
---|
1523 | naTest(a); |
---|
1524 | |
---|
1525 | /* special cases first */ |
---|
1526 | if (a == NULL) |
---|
1527 | { |
---|
1528 | if (exp >= 0) *b = NULL; |
---|
1529 | else WerrorS(nDivBy0); |
---|
1530 | return; |
---|
1531 | } |
---|
1532 | else if (exp == 0) { *b = naInit(1, cf); return; } |
---|
1533 | else if (exp == 1) { *b = naCopy(a, cf); return; } |
---|
1534 | else if (exp == -1) { *b = naInvers(a, cf); return; } |
---|
1535 | |
---|
1536 | int expAbs = exp; if (expAbs < 0) expAbs = -expAbs; |
---|
1537 | |
---|
1538 | /* now compute a^expAbs */ |
---|
1539 | poly pow; poly aAsPoly = (poly)a; |
---|
1540 | if (expAbs <= 7) |
---|
1541 | { |
---|
1542 | pow = p_Copy(aAsPoly, naRing); |
---|
1543 | for (int i = 2; i <= expAbs; i++) |
---|
1544 | { |
---|
1545 | pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing); |
---|
1546 | } |
---|
1547 | } |
---|
1548 | else |
---|
1549 | { |
---|
1550 | pow = p_ISet(1, naRing); |
---|
1551 | poly factor = p_Copy(aAsPoly, naRing); |
---|
1552 | while (expAbs != 0) |
---|
1553 | { |
---|
1554 | if (expAbs & 1) |
---|
1555 | { |
---|
1556 | pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing); |
---|
1557 | } |
---|
1558 | expAbs = expAbs / 2; |
---|
1559 | if (expAbs != 0) |
---|
1560 | { |
---|
1561 | factor = p_Mult_q(factor, p_Copy(factor, naRing), naRing); |
---|
1562 | } |
---|
1563 | } |
---|
1564 | p_Delete(&factor, naRing); |
---|
1565 | } |
---|
1566 | |
---|
1567 | /* invert if original exponent was negative */ |
---|
1568 | number n = (number)pow; |
---|
1569 | if (exp < 0) |
---|
1570 | { |
---|
1571 | number m = npolyInvers(n, cf); |
---|
1572 | naDelete(&n, cf); |
---|
1573 | n = m; |
---|
1574 | } |
---|
1575 | *b = n; |
---|
1576 | } |
---|
1577 | |
---|
1578 | number npolyDiv(number a, number b, const coeffs cf) |
---|
1579 | { |
---|
1580 | naTest(a); naTest(b); |
---|
1581 | if (b == NULL) WerrorS(nDivBy0); |
---|
1582 | if (a == NULL) return NULL; |
---|
1583 | poly p=singclap_pdivide((poly)a,(poly)b,naRing); |
---|
1584 | return (number)p; |
---|
1585 | } |
---|
1586 | |
---|
1587 | |
---|
1588 | BOOLEAN npolyInitChar(coeffs cf, void * infoStruct) |
---|
1589 | { |
---|
1590 | assume( infoStruct != NULL ); |
---|
1591 | |
---|
1592 | AlgExtInfo *e = (AlgExtInfo *)infoStruct; |
---|
1593 | /// first check whether cf->extRing != NULL and delete old ring??? |
---|
1594 | |
---|
1595 | assume(e->r != NULL); // extRing; |
---|
1596 | assume(e->r->cf != NULL); // extRing->cf; |
---|
1597 | |
---|
1598 | assume( cf != NULL ); |
---|
1599 | |
---|
1600 | e->r->ref ++; // increase the ref.counter for the ground poly. ring! |
---|
1601 | const ring R = e->r; // no copy! |
---|
1602 | cf->extRing = R; |
---|
1603 | |
---|
1604 | /* propagate characteristic up so that it becomes |
---|
1605 | directly accessible in cf: */ |
---|
1606 | cf->ch = R->cf->ch; |
---|
1607 | cf->is_field=FALSE; |
---|
1608 | cf->is_domain=TRUE; |
---|
1609 | |
---|
1610 | cf->cfCoeffString = naCoeffString; |
---|
1611 | |
---|
1612 | cf->cfGreaterZero = naGreaterZero; |
---|
1613 | cf->cfGreater = naGreater; |
---|
1614 | cf->cfEqual = naEqual; |
---|
1615 | cf->cfIsZero = naIsZero; |
---|
1616 | cf->cfIsOne = naIsOne; |
---|
1617 | cf->cfIsMOne = naIsMOne; |
---|
1618 | cf->cfInit = naInit; |
---|
1619 | cf->cfInit_bigint = naInit_bigint; |
---|
1620 | cf->cfFarey = naFarey; |
---|
1621 | cf->cfChineseRemainder= naChineseRemainder; |
---|
1622 | cf->cfInt = naInt; |
---|
1623 | cf->cfInpNeg = naNeg; |
---|
1624 | cf->cfAdd = naAdd; |
---|
1625 | cf->cfSub = naSub; |
---|
1626 | cf->cfMult = npolyMult; |
---|
1627 | cf->cfDiv = npolyDiv; |
---|
1628 | cf->cfPower = naPower; |
---|
1629 | cf->cfCopy = naCopy; |
---|
1630 | |
---|
1631 | cf->cfWriteLong = naWriteLong; |
---|
1632 | |
---|
1633 | if( rCanShortOut(naRing) ) |
---|
1634 | cf->cfWriteShort = naWriteShort; |
---|
1635 | else |
---|
1636 | cf->cfWriteShort = naWriteLong; |
---|
1637 | |
---|
1638 | cf->cfRead = naRead; |
---|
1639 | cf->cfDelete = naDelete; |
---|
1640 | cf->cfSetMap = naSetMap; |
---|
1641 | cf->cfGetDenom = naGetDenom; |
---|
1642 | cf->cfGetNumerator = naGetNumerator; |
---|
1643 | cf->cfRePart = naCopy; |
---|
1644 | cf->cfCoeffWrite = npolyCoeffWrite; |
---|
1645 | cf->cfNormalize = npolyNormalize; |
---|
1646 | cf->cfKillChar = naKillChar; |
---|
1647 | #ifdef LDEBUG |
---|
1648 | cf->cfDBTest = naDBTest; |
---|
1649 | #endif |
---|
1650 | cf->cfGcd = naGcd; |
---|
1651 | cf->cfLcm = naLcmContent; |
---|
1652 | cf->cfSize = naSize; |
---|
1653 | cf->nCoeffIsEqual = naCoeffIsEqual; |
---|
1654 | cf->cfInvers = npolyInvers; |
---|
1655 | cf->convFactoryNSingN=naConvFactoryNSingN; |
---|
1656 | cf->convSingNFactoryN=naConvSingNFactoryN; |
---|
1657 | cf->cfParDeg = naParDeg; |
---|
1658 | |
---|
1659 | cf->iNumberOfParameters = rVar(R); |
---|
1660 | cf->pParameterNames = (const char**)R->names; |
---|
1661 | cf->cfParameter = naParameter; |
---|
1662 | cf->has_simple_Inverse=FALSE; |
---|
1663 | /* cf->has_simple_Alloc= FALSE; */ |
---|
1664 | |
---|
1665 | if( nCoeff_is_Q(R->cf) ) |
---|
1666 | { |
---|
1667 | cf->cfClearContent = naClearContent; |
---|
1668 | cf->cfClearDenominators = naClearDenominators; |
---|
1669 | } |
---|
1670 | |
---|
1671 | return FALSE; |
---|
1672 | } |
---|
1673 | #endif |
---|