# source:git/libpolys/polys/ext_fields/algext.cc@558f3cc

spielwiese
Last change on this file since 558f3cc was 558f3cc, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Fix: make positive leading coeff at normalizing!
• Property mode set to `100644`
File size: 35.7 KB
Line
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 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#include "config.h"
32#include <misc/auxiliary.h>
33
34#include <omalloc/omalloc.h>
35
36#include <reporter/reporter.h>
37
38#include <coeffs/coeffs.h>
39#include <coeffs/numbers.h>
40#include <coeffs/longrat.h>
41
42#include <polys/monomials/ring.h>
43#include <polys/monomials/p_polys.h>
44#include <polys/simpleideals.h>
45
46#include <polys/PolyEnumerator.h>
47
48#ifdef HAVE_FACTORY
49#include <factory/factory.h>
50#include <polys/clapconv.h>
51#include <polys/clapsing.h>
52#endif
53
54
55#include <polys/ext_fields/algext.h>
56#define TRANSEXT_PRIVATES 1
57#include <polys/ext_fields/transext.h>
58
59#ifdef LDEBUG
61BOOLEAN  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!
67static 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
75 * coefficients of our numbers */
76#define naCoeffs cf->extRing->cf
77
78/* minimal polynomial */
79#define naMinpoly naRing->qideal->m
80
81/// forward declarations
82BOOLEAN  naGreaterZero(number a, const coeffs cf);
83BOOLEAN  naGreater(number a, number b, const coeffs cf);
84BOOLEAN  naEqual(number a, number b, const coeffs cf);
85BOOLEAN  naIsOne(number a, const coeffs cf);
86BOOLEAN  naIsMOne(number a, const coeffs cf);
87BOOLEAN  naIsZero(number a, const coeffs cf);
88number   naInit(long i, const coeffs cf);
89int      naInt(number &a, const coeffs cf);
90number   naNeg(number a, const coeffs cf);
91number   naInvers(number a, const coeffs cf);
92number   naAdd(number a, number b, const coeffs cf);
93number   naSub(number a, number b, const coeffs cf);
94number   naMult(number a, number b, const coeffs cf);
95number   naDiv(number a, number b, const coeffs cf);
96void     naPower(number a, int exp, number *b, const coeffs cf);
97number   naCopy(number a, const coeffs cf);
98void     naWriteLong(number &a, const coeffs cf);
99void     naWriteShort(number &a, const coeffs cf);
100number   naRePart(number a, const coeffs cf);
101number   naImPart(number a, const coeffs cf);
102number   naGetDenom(number &a, const coeffs cf);
103number   naGetNumerator(number &a, const coeffs cf);
104number   naGcd(number a, number b, const coeffs cf);
105//number   naLcm(number a, number b, const coeffs cf);
106int      naSize(number a, const coeffs cf);
107void     naDelete(number *a, const coeffs cf);
108void     naCoeffWrite(const coeffs cf, BOOLEAN details);
109//number   naIntDiv(number a, number b, const coeffs cf);
110const char * naRead(const char *s, number *a, const coeffs cf);
111
112static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param);
113
114
115/// returns NULL if p == NULL, otherwise makes p monic by dividing
116///   by its leading coefficient (only done if this is not already 1);
117///   this assumes that we are over a ground field so that division
118///   is well-defined;
119///   modifies p
120// void      p_Monic(poly p, const ring r);
121
122///   assumes that p and q are univariate polynomials in r,
123///   mentioning the same variable;
124///   assumes a global monomial ordering in r;
125///   assumes that not both p and q are NULL;
126///   returns the gcd of p and q;
127///   leaves p and q unmodified
128// poly      p_Gcd(const poly p, const poly q, const ring r);
129
130/* returns NULL if p == NULL, otherwise makes p monic by dividing
131   by its leading coefficient (only done if this is not already 1);
132   this assumes that we are over a ground field so that division
133   is well-defined;
134   modifies p */
135static inline void p_Monic(poly p, const ring r)
136{
137  if (p == NULL) return;
138  number n = n_Init(1, r->cf);
139  if (p->next==NULL) { p_SetCoeff(p,n,r); return; }
140  poly pp = p;
141  number lc = p_GetCoeff(p, r);
142  if (n_IsOne(lc, r->cf)) return;
143  number lcInverse = n_Invers(lc, r->cf);
144  p_SetCoeff(p, n, r);   // destroys old leading coefficient!
145  pIter(p);
146  while (p != NULL)
147  {
148    number n = n_Mult(p_GetCoeff(p, r), lcInverse, r->cf);
149    n_Normalize(n,r->cf);
150    p_SetCoeff(p, n, r);   // destroys old leading coefficient!
151    pIter(p);
152  }
153  n_Delete(&lcInverse, r->cf);
154  p = pp;
155}
156
157/// see p_Gcd;
158///   additional assumption: deg(p) >= deg(q);
159///   must destroy p and q (unless one of them is returned)
160static inline poly p_GcdHelper(poly &p, poly &q, const ring r)
161{
162  while (q != NULL)
163  {
164    p_PolyDiv(p, q, FALSE, r);
165    // swap p and q:
166    poly& t = q;
167    q = p;
168    p = t;
169
170  }
171  return p;
172}
173
174/* assumes that p and q are univariate polynomials in r,
175   mentioning the same variable;
176   assumes a global monomial ordering in r;
177   assumes that not both p and q are NULL;
178   returns the gcd of p and q;
179   leaves p and q unmodified */
180static inline poly      p_Gcd(const poly p, const poly q, const ring r)
181{
182  assume((p != NULL) || (q != NULL));
183
184  poly a = p; poly b = q;
185  if (p_Deg(a, r) < p_Deg(b, r)) { a = q; b = p; }
186  a = p_Copy(a, r); b = p_Copy(b, r);
187
188  /* We have to make p monic before we return it, so that if the
189     gcd is a unit in the ground field, we will actually return 1. */
190  a = p_GcdHelper(a, b, r);
191  p_Monic(a, r);
192  return a;
193}
194
195/* see p_ExtGcd;
196   additional assumption: deg(p) >= deg(q);
197   must destroy p and q (unless one of them is returned) */
198static inline poly p_ExtGcdHelper(poly &p, poly &pFactor, poly &q, poly &qFactor,
199                                  ring r)
200{
201  if (q == NULL)
202  {
203    qFactor = NULL;
204    pFactor = p_ISet(1, r);
205    p_SetCoeff(pFactor, n_Invers(p_GetCoeff(p, r), r->cf), r);
206    p_Monic(p, r);
207    return p;
208  }
209  else
210  {
211    poly pDivQ = p_PolyDiv(p, q, TRUE, r);
212    poly ppFactor = NULL; poly qqFactor = NULL;
213    poly theGcd = p_ExtGcdHelper(q, qqFactor, p, ppFactor, r);
214    pFactor = ppFactor;
216                      p_Neg(p_Mult_q(pDivQ, p_Copy(ppFactor, r), r), r),
217                      r);
218    return theGcd;
219  }
220}
221
222/* assumes that p and q are univariate polynomials in r,
223   mentioning the same variable;
224   assumes a global monomial ordering in r;
225   assumes that not both p and q are NULL;
226   returns the gcd of p and q;
227   moreover, afterwards pFactor and qFactor contain appropriate
228   factors such that gcd(p, q) = p * pFactor + q * qFactor;
229   leaves p and q unmodified */
230poly p_ExtGcd(poly p, poly &pFactor, poly q, poly &qFactor, ring r)
231{
232  assume((p != NULL) || (q != NULL));
233  poly a = p; poly b = q; BOOLEAN aCorrespondsToP = TRUE;
234  if (p_Deg(a, r) < p_Deg(b, r))
235    { a = q; b = p; aCorrespondsToP = FALSE; }
236  a = p_Copy(a, r); b = p_Copy(b, r);
237  poly aFactor = NULL; poly bFactor = NULL;
238  poly theGcd = p_ExtGcdHelper(a, aFactor, b, bFactor, r);
239  if (aCorrespondsToP) { pFactor = aFactor; qFactor = bFactor; }
240  else                 { pFactor = bFactor; qFactor = aFactor; }
241  return theGcd;
242}
243
244
245
246#ifdef LDEBUG
247BOOLEAN naDBTest(number a, const char *f, const int l, const coeffs cf)
248{
249  assume(getCoeffType(cf) == ID);
250  if (a == NULL) return TRUE;
251  p_Test((poly)a, naRing);
252  if((((poly)a)!=naMinpoly)
253  && p_Totaldegree((poly)a, naRing) >= p_Totaldegree(naMinpoly, naRing))
254  {
255    Print("deg >= deg(minpoly) in %s:%d\n",f,l);
256    return FALSE;
257  }
258  return TRUE;
259}
260#endif
261
262void heuristicReduce(poly &p, poly reducer, const coeffs cf);
263void definiteReduce(poly &p, poly reducer, const coeffs cf);
264
265/* returns the bottom field in this field extension tower; if the tower
266   is flat, i.e., if there is no extension, then r itself is returned;
267   as a side-effect, the counter 'height' is filled with the height of
268   the extension tower (in case the tower is flat, 'height' is zero) */
269static coeffs nCoeff_bottom(const coeffs r, int &height)
270{
271  assume(r != NULL);
272  coeffs cf = r;
273  height = 0;
274  while (nCoeff_is_Extension(cf))
275  {
276    assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL);
277    cf = cf->extRing->cf;
278    height++;
279  }
280  return cf;
281}
282
283BOOLEAN naIsZero(number a, const coeffs cf)
284{
285  naTest(a);
286  return (a == NULL);
287}
288
289void naDelete(number * a, const coeffs cf)
290{
291  if (*a == NULL) return;
292  if (((poly)*a)==naMinpoly) { *a=NULL;return;}
293  poly aAsPoly = (poly)(*a);
294  p_Delete(&aAsPoly, naRing);
295  *a = NULL;
296}
297
298BOOLEAN naEqual(number a, number b, const coeffs cf)
299{
300  naTest(a); naTest(b);
301  /// simple tests
302  if (a == NULL) return (b == NULL);
303  if (b == NULL) return (a == NULL);
304  return p_EqualPolys((poly)a,(poly)b,naRing);
305}
306
307number naCopy(number a, const coeffs cf)
308{
309  naTest(a);
310  if (a == NULL) return NULL;
311  if (((poly)a)==naMinpoly) return a;
312  return (number)p_Copy((poly)a, naRing);
313}
314
315number naGetNumerator(number &a, const coeffs cf)
316{
317  return naCopy(a, cf);
318}
319
320number naGetDenom(number &a, const coeffs cf)
321{
322  naTest(a);
323  return naInit(1, cf);
324}
325
326BOOLEAN naIsOne(number a, const coeffs cf)
327{
328  naTest(a);
329  poly aAsPoly = (poly)a;
330  if ((a==NULL) || (!p_IsConstant(aAsPoly, naRing))) return FALSE;
331  return n_IsOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
332}
333
334BOOLEAN naIsMOne(number a, const coeffs cf)
335{
336  naTest(a);
337  poly aAsPoly = (poly)a;
338  if ((a==NULL) || (!p_IsConstant(aAsPoly, naRing))) return FALSE;
339  return n_IsMOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
340}
341
342/// this is in-place, modifies a
343number naNeg(number a, const coeffs cf)
344{
345  naTest(a);
346  if (a != NULL) a = (number)p_Neg((poly)a, naRing);
347  return a;
348}
349
350number naImPart(number a, const coeffs cf)
351{
352  naTest(a);
353  return NULL;
354}
355
356number naInit_bigint(number longratBigIntNumber, const coeffs src, const coeffs cf)
357{
358  assume( cf != NULL );
359
360  const ring A = cf->extRing;
361
362  assume( A != NULL );
363
364  const coeffs C = A->cf;
365
366  assume( C != NULL );
367
368  number n = n_Init_bigint(longratBigIntNumber, src, C);
369
370  if ( n_IsZero(n, C) )
371  {
372    n_Delete(&n, C);
373    return NULL;
374  }
375
376  return (number)p_NSet(n, A);
377}
378
379
380
381number naInit(long i, const coeffs cf)
382{
383  if (i == 0) return NULL;
384  else        return (number)p_ISet(i, naRing);
385}
386
387int naInt(number &a, const coeffs cf)
388{
389  naTest(a);
390  poly aAsPoly = (poly)a;
391  if(aAsPoly == NULL)
392    return 0;
393  if (!p_IsConstant(aAsPoly, naRing))
394    return 0;
395  assume( aAsPoly != NULL );
396  return n_Int(p_GetCoeff(aAsPoly, naRing), naCoeffs);
397}
398
399/* TRUE iff (a != 0 and (b == 0 or deg(a) > deg(b) or (deg(a)==deg(b) && lc(a)>lc(b))) */
400BOOLEAN naGreater(number a, number b, const coeffs cf)
401{
402  naTest(a); naTest(b);
403  if (naIsZero(a, cf))
404  {
405    if (naIsZero(b, cf)) return FALSE;
406    return !n_GreaterZero(pGetCoeff((poly)b),cf);
407  }
408  if (naIsZero(b, cf))
409  {
410    return n_GreaterZero(pGetCoeff((poly)a),cf);
411  }
412  int aDeg = p_Totaldegree((poly)a, naRing);
413  int bDeg = p_Totaldegree((poly)b, naRing);
416  return n_Greater(pGetCoeff((poly)a),pGetCoeff((poly)b),cf);
417}
418
419/* TRUE iff a != 0 and (LC(a) > 0 or deg(a) > 0) */
420BOOLEAN naGreaterZero(number a, const coeffs cf)
421{
422  naTest(a);
423  if (a == NULL)                                            return FALSE;
424  if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE;
425  if (p_Totaldegree((poly)a, naRing) > 0)                   return TRUE;
426  return FALSE;
427}
428
429void naCoeffWrite(const coeffs cf, BOOLEAN details)
430{
431  assume( cf != NULL );
432
433  const ring A = cf->extRing;
434
435  assume( A != NULL );
436  assume( A->cf != NULL );
437
438  n_CoeffWrite(A->cf, details);
439
440//  rWrite(A);
441
442  const int P = rVar(A);
443  assume( P > 0 );
444
445  Print("//   %d parameter    : ", P);
446
447  for (int nop=0; nop < P; nop ++)
448    Print("%s ", rRingVar(nop, A));
449
450  PrintLn();
451
452  const ideal I = A->qideal;
453
454  assume( I != NULL );
455  assume( IDELEMS(I) == 1 );
456
457
458  if ( details )
459  {
460    PrintS("//   minpoly        : (");
461    p_Write0( I->m, A);
462    PrintS(")");
463  }
464  else
465    PrintS("//   minpoly        : ...");
466
467  PrintLn();
468
469/*
470  char *x = rRingVar(0, A);
471
472  Print("//   Coefficients live in the extension field K[%s]/<f(%s)>\n", x, x);
473  Print("//   with the minimal polynomial f(%s) = %s\n", x,
474        p_String(A->qideal->m, A));
475  PrintS("//   and K: ");
476*/
477}
478
479number naAdd(number a, number b, const coeffs cf)
480{
481  naTest(a); naTest(b);
482  if (a == NULL) return naCopy(b, cf);
483  if (b == NULL) return naCopy(a, cf);
484  poly aPlusB = p_Add_q(p_Copy((poly)a, naRing),
485                        p_Copy((poly)b, naRing), naRing);
486  definiteReduce(aPlusB, naMinpoly, cf);
487  return (number)aPlusB;
488}
489
490number naSub(number a, number b, const coeffs cf)
491{
492  naTest(a); naTest(b);
493  if (b == NULL) return naCopy(a, cf);
494  poly minusB = p_Neg(p_Copy((poly)b, naRing), naRing);
495  if (a == NULL) return (number)minusB;
496  poly aMinusB = p_Add_q(p_Copy((poly)a, naRing), minusB, naRing);
497  definiteReduce(aMinusB, naMinpoly, cf);
498  return (number)aMinusB;
499}
500
501number naMult(number a, number b, const coeffs cf)
502{
503  naTest(a); naTest(b);
504  if ((a == NULL)||(b == NULL)) return NULL;
505  poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing),
506                          p_Copy((poly)b, naRing), naRing);
507  p_Normalize(aTimesB,naRing);
508  definiteReduce(aTimesB, naMinpoly, cf);
509  return (number)aTimesB;
510}
511
512number naDiv(number a, number b, const coeffs cf)
513{
514  naTest(a); naTest(b);
515  if (b == NULL) WerrorS(nDivBy0);
516  if (a == NULL) return NULL;
517  poly bInverse = (poly)naInvers(b, cf);
518  if(bInverse != NULL) // b is non-zero divisor!
519  {
520    poly aDivB = p_Mult_q(p_Copy((poly)a, naRing), bInverse, naRing);
523  }
524  return NULL;
525}
526
527/* 0^0 = 0;
528   for |exp| <= 7 compute power by a simple multiplication loop;
529   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
530      p^13 = p^1 * p^4 * p^8, where we utilise that
531      p^(2^(k+1)) = p^(2^k) * p^(2^k);
532   intermediate reduction modulo the minimal polynomial is controlled by
533   the in-place method heuristicReduce(poly, poly, coeffs); see there.
534*/
535void naPower(number a, int exp, number *b, const coeffs cf)
536{
537  naTest(a);
538
539  /* special cases first */
540  if (a == NULL)
541  {
542    if (exp >= 0) *b = NULL;
543    else          WerrorS(nDivBy0);
544    return;
545  }
546  else if (exp ==  0) { *b = naInit(1, cf); return; }
547  else if (exp ==  1) { *b = naCopy(a, cf); return; }
548  else if (exp == -1) { *b = naInvers(a, cf); return; }
549
550  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
551
552  /* now compute a^expAbs */
553  poly pow; poly aAsPoly = (poly)a;
554  if (expAbs <= 7)
555  {
556    pow = p_Copy(aAsPoly, naRing);
557    for (int i = 2; i <= expAbs; i++)
558    {
559      pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing);
560      heuristicReduce(pow, naMinpoly, cf);
561    }
562    definiteReduce(pow, naMinpoly, cf);
563  }
564  else
565  {
566    pow = p_ISet(1, naRing);
567    poly factor = p_Copy(aAsPoly, naRing);
568    while (expAbs != 0)
569    {
570      if (expAbs & 1)
571      {
572        pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing);
573        heuristicReduce(pow, naMinpoly, cf);
574      }
575      expAbs = expAbs / 2;
576      if (expAbs != 0)
577      {
578        factor = p_Mult_q(factor, p_Copy(factor, naRing), naRing);
579        heuristicReduce(factor, naMinpoly, cf);
580      }
581    }
582    p_Delete(&factor, naRing);
583    definiteReduce(pow, naMinpoly, cf);
584  }
585
586  /* invert if original exponent was negative */
587  number n = (number)pow;
588  if (exp < 0)
589  {
590    number m = naInvers(n, cf);
592    n = m;
593  }
594  *b = n;
595}
596
597/* may reduce p modulo the reducer by calling definiteReduce;
598   the decision is made based on the following heuristic
599   (which should also only be changed here in this method):
600      if (deg(p) > 10*deg(reducer) then perform reduction;
601   modifies p */
602void heuristicReduce(poly &p, poly reducer, const coeffs cf)
603{
604  #ifdef LDEBUG
605  p_Test((poly)p, naRing);
606  p_Test((poly)reducer, naRing);
607  #endif
608  if (p_Totaldegree(p, naRing) > 10 * p_Totaldegree(reducer, naRing))
609    definiteReduce(p, reducer, cf);
610}
611
612void naWriteLong(number &a, const coeffs cf)
613{
614  naTest(a);
615  if (a == NULL)
616    StringAppendS("0");
617  else
618  {
619    poly aAsPoly = (poly)a;
620    /* basically, just write aAsPoly using p_Write,
621       but use brackets around the output, if a is not
622       a constant living in naCoeffs = cf->extRing->cf */
623    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
624    if (useBrackets) StringAppendS("(");
625    p_String0Long(aAsPoly, naRing, naRing);
626    if (useBrackets) StringAppendS(")");
627  }
628}
629
630void naWriteShort(number &a, const coeffs cf)
631{
632  naTest(a);
633  if (a == NULL)
634    StringAppendS("0");
635  else
636  {
637    poly aAsPoly = (poly)a;
638    /* basically, just write aAsPoly using p_Write,
639       but use brackets around the output, if a is not
640       a constant living in naCoeffs = cf->extRing->cf */
641    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
642    if (useBrackets) StringAppendS("(");
643    p_String0Short(aAsPoly, naRing, naRing);
644    if (useBrackets) StringAppendS(")");
645  }
646}
647
648const char * naRead(const char *s, number *a, const coeffs cf)
649{
650  poly aAsPoly;
651  const char * result = p_Read(s, aAsPoly, naRing);
652  definiteReduce(aAsPoly, naMinpoly, cf);
653  *a = (number)aAsPoly;
654  return result;
655}
656
657#if 0
658/* implemented by the rule lcm(a, b) = a * b / gcd(a, b) */
659number naLcm(number a, number b, const coeffs cf)
660{
661  naTest(a); naTest(b);
662  if (a == NULL) return NULL;
663  if (b == NULL) return NULL;
664  number theProduct = (number)p_Mult_q(p_Copy((poly)a, naRing),
665                                       p_Copy((poly)b, naRing), naRing);
666  /* note that theProduct needs not be reduced w.r.t. naMinpoly;
667     but the final division will take care of the necessary reduction */
668  number theGcd = naGcd(a, b, cf);
670}
671#endif
672
673/* expects *param to be castable to AlgExtInfo */
674static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
675{
676  if (ID != n) return FALSE;
677  AlgExtInfo *e = (AlgExtInfo *)param;
678  /* for extension coefficient fields we expect the underlying
679     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
680     this expectation is based on the assumption that we have properly
681     registered cf and perform reference counting rather than creating
682     multiple copies of the same coefficient field/domain/ring */
683  if (naRing == e->r)
684    return TRUE;
685  /* (Note that then also the minimal ideals will necessarily be
686     the same, as they are attached to the ring.) */
687
688  // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)...
689  if( rEqual(naRing, e->r, TRUE) ) // also checks the equality of qideals
690  {
691    const ideal mi = naRing->qideal;
692    assume( IDELEMS(mi) == 1 );
693    const ideal ii = e->r->qideal;
694    assume( IDELEMS(ii) == 1 );
695
696    // TODO: the following should be extended for 2 *equal* rings...
697    assume( p_EqualPolys(mi->m, ii->m, naRing, e->r) );
698
699    rDelete(e->r);
700
701    return TRUE;
702  }
703
704  return FALSE;
705
706}
707
708int naSize(number a, const coeffs cf)
709{
710  if (a == NULL) return -1;
711  /* this has been taken from the old implementation of field extensions,
712     where we computed the sum of the degree and the number of terms in
713     (poly)a; so we leave it at that, for the time being;
714     maybe, the number of terms alone is a better measure? */
715  poly aAsPoly = (poly)a;
716  int theDegree = 0; int noOfTerms = 0;
717  while (aAsPoly != NULL)
718  {
719    noOfTerms++;
720    int d = p_GetExp(aAsPoly, 1, naRing);
721    if (d > theDegree) theDegree = d;
722    pIter(aAsPoly);
723  }
724  return theDegree + noOfTerms;
725}
726
727/* performs polynomial division and overrides p by the remainder
728   of division of p by the reducer;
729   modifies p */
730void definiteReduce(poly &p, poly reducer, const coeffs cf)
731{
732  #ifdef LDEBUG
733  p_Test((poly)p, naRing);
734  p_Test((poly)reducer, naRing);
735  #endif
736  if ((p!=NULL) && (p_GetExp(p,1,naRing)>=p_GetExp(reducer,1,naRing)))
737  {
738    p_PolyDiv(p, reducer, FALSE, naRing);
739  }
740}
741
742void  naNormalize(number &a, const coeffs cf)
743{
744  poly aa=(poly)a;
745  if (aa!=naMinpoly)
746    definiteReduce(aa,naMinpoly,cf);
747  a=(number)aa;
748}
749
750#ifdef HAVE_FACTORY
751number naConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
752{
753  if (n.isZero()) return NULL;
754  poly p=convFactoryPSingP(n,naRing);
755  return (number)p;
756}
757CanonicalForm naConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
758{
759  naTest(n);
760  if (n==NULL) return CanonicalForm(0);
761
762  return convSingPFactoryP((poly)n,naRing);
763}
764#endif
765
766
767/* IMPORTANT NOTE: Since an algebraic field extension is again a field,
768                   the gcd of two elements is not very interesting. (It
769                   is actually any unit in the field, i.e., any non-
770                   zero element.) Note that the below method does not operate
771                   in this strong sense but rather computes the gcd of
772                   two given elements in the underlying polynomial ring. */
773number naGcd(number a, number b, const coeffs cf)
774{
775  naTest(a); naTest(b);
776  if ((a == NULL) && (b == NULL)) WerrorS(nDivBy0);
777  const ring R = naRing;
778  return (number) singclap_gcd(p_Copy((poly)a, R), p_Copy((poly)b, R), R);
779//  return (number)p_Gcd((poly)a, (poly)b, naRing);
780}
781
782number naInvers(number a, const coeffs cf)
783{
784  naTest(a);
785  if (a == NULL) WerrorS(nDivBy0);
786
787  poly aFactor = NULL; poly mFactor = NULL; poly theGcd = NULL;
788// singclap_extgcd!
789  const BOOLEAN ret = singclap_extgcd ((poly)a, naMinpoly, theGcd, aFactor, mFactor, naRing);
790
791  assume( !ret );
792
793//  if( ret ) theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
794
795  naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor);
796  p_Delete(&mFactor, naRing);
797
798  //  /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */
799  //  assume(naIsOne((number)theGcd, cf));
800
801  if( !naIsOne((number)theGcd, cf) )
802  {
803    WerrorS("zero divisor found - your minpoly is not irreducible");
804    p_Delete(&aFactor, naRing); aFactor = NULL;
805  }
806  p_Delete(&theGcd, naRing);
807
808  return (number)(aFactor);
809}
810
811/* assumes that src = Q, dst = Q(a) */
812number naMap00(number a, const coeffs src, const coeffs dst)
813{
814  if (n_IsZero(a, src)) return NULL;
815  assume(src == dst->extRing->cf);
816  poly result = p_One(dst->extRing);
817  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
818  return (number)result;
819}
820
821/* assumes that src = Z/p, dst = Q(a) */
822number naMapP0(number a, const coeffs src, const coeffs dst)
823{
824  if (n_IsZero(a, src)) return NULL;
825  /* mapping via intermediate int: */
826  int n = n_Int(a, src);
827  number q = n_Init(n, dst->extRing->cf);
828  poly result = p_One(dst->extRing);
829  p_SetCoeff(result, q, dst->extRing);
830  return (number)result;
831}
832
833#if 0
834/* assumes that either src = Q(a), dst = Q(a), or
835                       src = Z/p(a), dst = Z/p(a)     */
836number naCopyMap(number a, const coeffs src, const coeffs dst)
837{
838  return naCopy(a, dst);
839}
840#endif
841
842number naCopyExt(number a, const coeffs src, const coeffs)
843{
844  fraction fa=(fraction)a;
845  return (number)p_Copy(NUM(fa),src->extRing);
846}
847
848/* assumes that src = Q, dst = Z/p(a) */
849number naMap0P(number a, const coeffs src, const coeffs dst)
850{
851  if (n_IsZero(a, src)) return NULL;
852  int p = rChar(dst->extRing);
853
854  number q = nlModP(a, src, dst->extRing->cf);
855
856  poly result = p_NSet(q, dst->extRing);
857
858  return (number)result;
859}
860
861/* assumes that src = Z/p, dst = Z/p(a) */
862number naMapPP(number a, const coeffs src, const coeffs dst)
863{
864  if (n_IsZero(a, src)) return NULL;
865  assume(src == dst->extRing->cf);
866  poly result = p_One(dst->extRing);
867  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
868  return (number)result;
869}
870
871/* assumes that src = Z/u, dst = Z/p(a), where u != p */
872number naMapUP(number a, const coeffs src, const coeffs dst)
873{
874  if (n_IsZero(a, src)) return NULL;
875  /* mapping via intermediate int: */
876  int n = n_Int(a, src);
877  number q = n_Init(n, dst->extRing->cf);
878  poly result = p_One(dst->extRing);
879  p_SetCoeff(result, q, dst->extRing);
880  return (number)result;
881}
882
883nMapFunc naSetMap(const coeffs src, const coeffs dst)
884{
885  /* dst is expected to be an algebraic field extension */
886  assume(getCoeffType(dst) == ID);
887
888  if( src == dst ) return ndCopyMap;
889
890  int h = 0; /* the height of the extension tower given by dst */
891  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
892  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
893
894  /* for the time being, we only provide maps if h = 1 and if b is Q or
895     some field Z/pZ: */
896  if (h==0)
897  {
898    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
899      return naMap00;                            /// Q     -->  Q(a)
900    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
901      return naMapP0;                            /// Z/p   -->  Q(a)
902    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
903      return naMap0P;                            /// Q      --> Z/p(a)
904    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
905    {
906      if (src->ch == dst->ch) return naMapPP;    /// Z/p    --> Z/p(a)
907      else return naMapUP;                       /// Z/u    --> Z/p(a)
908    }
909  }
910  if (h != 1) return NULL;
911  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
912  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
913
914  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
915  {
916    if (strcmp(rRingVar(0, src->extRing),
917               rRingVar(0, dst->extRing)) == 0)
918    {
919      if (src->type==n_algExt)
920         return ndCopyMap; // naCopyMap;         /// Q(a)   --> Q(a)
921      else
922         return naCopyExt;
923    }
924    else
925      return NULL;                               /// Q(b)   --> Q(a)
926  }
927
928  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
929  {
930    if (strcmp(rRingVar(0,src->extRing),rRingVar(0,dst->extRing))==0)
931    {
932      if (src->type==n_algExt)
933        return ndCopyMap; // naCopyMap;          /// Z/p(a) --> Z/p(a)
934      else
935         return naCopyExt;
936    }
937    else
938      return NULL;                               /// Z/p(b) --> Z/p(a)
939  }
940
941  return NULL;                                           /// default
942}
943
944static int naParDeg(number a, const coeffs cf)
945{
946  if (a == NULL) return -1;
947  poly aa=(poly)a;
948  return cf->extRing->pFDeg(aa,cf->extRing);
949}
950
951/// return the specified parameter as a number in the given alg. field
952static number naParameter(const int iParameter, const coeffs cf)
953{
954  assume(getCoeffType(cf) == ID);
955
956  const ring R = cf->extRing;
957  assume( R != NULL );
958  assume( 0 < iParameter && iParameter <= rVar(R) );
959
960  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
961
962  return (number) p;
963}
964
965
966/// if m == var(i)/1 => return i,
967int naIsParam(number m, const coeffs cf)
968{
969  assume(getCoeffType(cf) == ID);
970
971  const ring R = cf->extRing;
972  assume( R != NULL );
973
974  return p_Var( (poly)m, R );
975}
976
977
978static void naClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
979{
980  assume(cf != NULL);
981  assume(getCoeffType(cf) == ID);
982  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
983
984  const ring   R = cf->extRing;
985  assume(R != NULL);
986  const coeffs Q = R->cf;
987  assume(Q != NULL);
988  assume(nCoeff_is_Q(Q));
989
990  numberCollectionEnumerator.Reset();
991
992  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
993  {
994    c = n_Init(1, cf);
995    return;
996  }
997
998  naTest(numberCollectionEnumerator.Current());
999
1000  // part 1, find a small candidate for gcd
1001  int s1; int s=2147483647; // max. int
1002
1003  const BOOLEAN lc_is_pos=naGreaterZero(numberCollectionEnumerator.Current(),cf);
1004
1005  int normalcount = 0;
1006
1007  poly cand1, cand;
1008
1009  do
1010  {
1011    number& n = numberCollectionEnumerator.Current();
1012    naNormalize(n, cf); ++normalcount;
1013
1014    naTest(n);
1015
1016    cand1 = (poly)n;
1017
1018    s1 = p_Deg(cand1, R); // naSize?
1019    if (s>s1)
1020    {
1021      cand = cand1;
1022      s = s1;
1023    }
1024  } while (numberCollectionEnumerator.MoveNext() );
1025
1026//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
1027
1028  cand = p_Copy(cand, R);
1029  // part 2: compute gcd(cand,all coeffs)
1030
1031  numberCollectionEnumerator.Reset();
1032
1033  int length = 0;
1034  while (numberCollectionEnumerator.MoveNext() )
1035  {
1036    number& n = numberCollectionEnumerator.Current();
1037    ++length;
1038
1039    if( (--normalcount) <= 0)
1040      naNormalize(n, cf);
1041
1042    naTest(n);
1043
1044//    p_InpGcd(cand, (poly)n, R);
1045
1046    cand = singclap_gcd(cand, p_Copy((poly)n, R), R);
1047
1048//    cand1 = p_Gcd(cand,(poly)n, R); p_Delete(&cand, R); cand = cand1;
1049
1050    assume( naGreaterZero((number)cand, cf) ); // ???
1051/*
1052    if(p_IsConstant(cand,R))
1053    {
1054      c = cand;
1055
1056      if(!lc_is_pos)
1057      {
1058        // make the leading coeff positive
1059        c = nlNeg(c, cf);
1060        numberCollectionEnumerator.Reset();
1061
1062        while (numberCollectionEnumerator.MoveNext() )
1063        {
1064          number& nn = numberCollectionEnumerator.Current();
1065          nn = nlNeg(nn, cf);
1066        }
1067      }
1068      return;
1069    }
1070*/
1071
1072  }
1073
1074  // part3: all coeffs = all coeffs / cand
1075  if (!lc_is_pos)
1076    cand = p_Neg(cand, R);
1077
1078  c = (number)cand; naTest(c);
1079
1080  poly cInverse = (poly)naInvers(c, cf);
1081  assume(cInverse != NULL); // c is non-zero divisor!?
1082
1083
1084  numberCollectionEnumerator.Reset();
1085
1086
1087  while (numberCollectionEnumerator.MoveNext() )
1088  {
1089    number& n = numberCollectionEnumerator.Current();
1090
1091    assume( length > 0 );
1092
1093    if( --length > 0 )
1094    {
1095      assume( cInverse != NULL );
1096      n = (number) p_Mult_q(p_Copy(cInverse, R), (poly)n, R);
1097    }
1098    else
1099    {
1100      n = (number) p_Mult_q(cInverse, (poly)n, R);
1101      cInverse = NULL;
1102      assume(length == 0);
1103    }
1104
1105    definiteReduce((poly &)n, naMinpoly, cf);
1106  }
1107
1108  assume(length == 0);
1109  assume(cInverse == NULL); //   p_Delete(&cInverse, R);
1110
1111  // Quick and dirty fix for constant content clearing... !?
1112  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
1113
1114  number cc;
1115
1116  n_ClearContent(itr, cc, Q); // TODO: get rid of (-LC) normalization!?
1117
1118  // over alg. ext. of Q // takes over the input number
1119  c = (number) p_Mult_nn( (poly)c, cc, R);
1120//      p_Mult_q(p_NSet(cc, R), , R);
1121
1122  n_Delete(&cc, Q);
1123
1124  // TODO: the above is not enough! need GCD's of polynomial coeffs...!
1125/*
1126  // old and wrong part of p_Content
1127    if (rField_is_Q_a(r) && !CLEARENUMERATORS) // should not be used anymore if CLEARENUMERATORS is 1
1128    {
1129      // we only need special handling for alg. ext.
1130      if (getCoeffType(r->cf)==n_algExt)
1131      {
1132        number hzz = n_Init(1, r->cf->extRing->cf);
1133        p=ph;
1134        while (p!=NULL)
1135        { // each monom: coeff in Q_a
1136          poly c_n_n=(poly)pGetCoeff(p);
1137          poly c_n=c_n_n;
1138          while (c_n!=NULL)
1139          { // each monom: coeff in Q
1140            d=n_Lcm(hzz,pGetCoeff(c_n),r->cf->extRing->cf);
1141            n_Delete(&hzz,r->cf->extRing->cf);
1142            hzz=d;
1143            pIter(c_n);
1144          }
1145          pIter(p);
1146        }
1147        // hzz contains the 1/lcm of all denominators in c_n_n
1148        h=n_Invers(hzz,r->cf->extRing->cf);
1149        n_Delete(&hzz,r->cf->extRing->cf);
1150        n_Normalize(h,r->cf->extRing->cf);
1151        if(!n_IsOne(h,r->cf->extRing->cf))
1152        {
1153          p=ph;
1154          while (p!=NULL)
1155          { // each monom: coeff in Q_a
1156            poly c_n=(poly)pGetCoeff(p);
1157            while (c_n!=NULL)
1158            { // each monom: coeff in Q
1159              d=n_Mult(h,pGetCoeff(c_n),r->cf->extRing->cf);
1160              n_Normalize(d,r->cf->extRing->cf);
1161              n_Delete(&pGetCoeff(c_n),r->cf->extRing->cf);
1162              pGetCoeff(c_n)=d;
1163              pIter(c_n);
1164            }
1165            pIter(p);
1166          }
1167        }
1168        n_Delete(&h,r->cf->extRing->cf);
1169      }
1170    }
1171*/
1172
1173
1174//  c = n_Init(1, cf); assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
1175}
1176
1177
1178static void naClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1179{
1180  assume(cf != NULL);
1181  assume(getCoeffType(cf) == ID);
1182  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
1183
1184  assume(cf->extRing != NULL);
1185  const coeffs Q = cf->extRing->cf;
1186  assume(Q != NULL);
1187  assume(nCoeff_is_Q(Q));
1188  number n;
1189  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
1190  n_ClearDenominators(itr, n, Q); // this should probably be fine...
1191  c = (number)p_NSet(n, cf->extRing); // over alg. ext. of Q // takes over the input number
1192}
1193
1194
1195BOOLEAN naInitChar(coeffs cf, void * infoStruct)
1196{
1197  assume( infoStruct != NULL );
1198
1199  AlgExtInfo *e = (AlgExtInfo *)infoStruct;
1200  /// first check whether cf->extRing != NULL and delete old ring???
1201
1202  assume(e->r                     != NULL);      // extRing;
1203  assume(e->r->cf                 != NULL);      // extRing->cf;
1204
1205  assume((e->r->qideal            != NULL) &&    // minideal has one
1206         (IDELEMS(e->r->qideal)   == 1)    &&    // non-zero generator
1207         (e->r->qideal->m      != NULL)    ); // at m;
1208
1209  assume( cf != NULL );
1210  assume(getCoeffType(cf) == ID);                     // coeff type;
1211
1212  e->r->ref ++; // increase the ref.counter for the ground poly. ring!
1213  const ring R = e->r; // no copy!
1214  assume( R->qideal == e->r->qideal );
1215  cf->extRing  = R;
1216
1217  /* propagate characteristic up so that it becomes
1218     directly accessible in cf: */
1219  cf->ch = R->cf->ch;
1220
1221  #ifdef LDEBUG
1222  p_Test((poly)naMinpoly, naRing);
1223  #endif
1224
1225  cf->cfGreaterZero  = naGreaterZero;
1226  cf->cfGreater      = naGreater;
1227  cf->cfEqual        = naEqual;
1228  cf->cfIsZero       = naIsZero;
1229  cf->cfIsOne        = naIsOne;
1230  cf->cfIsMOne       = naIsMOne;
1231  cf->cfInit         = naInit;
1232  cf->cfInit_bigint  = naInit_bigint;
1233  cf->cfInt          = naInt;
1234  cf->cfNeg          = naNeg;
1236  cf->cfSub          = naSub;
1237  cf->cfMult         = naMult;
1240  cf->cfPower        = naPower;
1241  cf->cfCopy         = naCopy;
1242
1243  cf->cfWriteLong        = naWriteLong;
1244
1245  if( rCanShortOut(naRing) )
1246    cf->cfWriteShort = naWriteShort;
1247  else
1248    cf->cfWriteShort = naWriteLong;
1249
1252  cf->cfSetMap       = naSetMap;
1253  cf->cfGetDenom     = naGetDenom;
1254  cf->cfGetNumerator = naGetNumerator;
1255  cf->cfRePart       = naCopy;
1256  cf->cfImPart       = naImPart;
1257  cf->cfCoeffWrite   = naCoeffWrite;
1258  cf->cfNormalize    = naNormalize;
1259#ifdef LDEBUG
1261#endif
1262  cf->cfGcd          = naGcd;
1263  //cf->cfLcm          = naLcm;
1264  cf->cfSize         = naSize;
1265  cf->nCoeffIsEqual  = naCoeffIsEqual;
1266  cf->cfInvers       = naInvers;
1267  cf->cfIntDiv       = naDiv; // ???
1268#ifdef HAVE_FACTORY
1269  cf->convFactoryNSingN=naConvFactoryNSingN;
1270  cf->convSingNFactoryN=naConvSingNFactoryN;
1271#endif
1272  cf->cfParDeg = naParDeg;
1273
1274  cf->iNumberOfParameters = rVar(R);
1275  cf->pParameterNames = R->names;
1276  cf->cfParameter = naParameter;
1277
1278  if( nCoeff_is_Q(R->cf) )
1279  {
1280    cf->cfClearContent = naClearContent;
1281    cf->cfClearDenominators = naClearDenominators;
1282  }
1283
1284  return FALSE;
1285}
Note: See TracBrowser for help on using the repository browser.