# source:git/libpolys/polys/ext_fields/algext.cc@2466d22

spielwiese
Last change on this file since 2466d22 was 2466d22, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: naPower(0)
• Property mode set to `100644`
File size: 36.2 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) ((void)(TRUE))
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
302  /// simple tests
303  if (a == b) return TRUE;
304  if ((a == NULL) && (b != NULL)) return FALSE;
305  if ((b == NULL) && (a != NULL)) return FALSE;
306
307  /// deg test
309  if (a != NULL) aDeg = p_Totaldegree((poly)a, naRing);
310  int bDeg = 0;
311  if (b != NULL) bDeg = p_Totaldegree((poly)b, naRing);
312  if (aDeg != bDeg) return FALSE;
313
314  /// subtraction test
315  number c = naSub(a, b, cf);
316  BOOLEAN result = naIsZero(c, cf);
318  return result;
319}
320
321number naCopy(number a, const coeffs cf)
322{
323  naTest(a);
324  if (a == NULL) return NULL;
325  if (((poly)a)==naMinpoly) return a;
326  return (number)p_Copy((poly)a, naRing);
327}
328
329number naGetNumerator(number &a, const coeffs cf)
330{
331  return naCopy(a, cf);
332}
333
334number naGetDenom(number &a, const coeffs cf)
335{
336  naTest(a);
337  return naInit(1, cf);
338}
339
340BOOLEAN naIsOne(number a, const coeffs cf)
341{
342  naTest(a);
343  poly aAsPoly = (poly)a;
344  if ((a==NULL) || (!p_IsConstant(aAsPoly, naRing))) return FALSE;
345  return n_IsOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
346}
347
348BOOLEAN naIsMOne(number a, const coeffs cf)
349{
350  naTest(a);
351  poly aAsPoly = (poly)a;
352  if ((a==NULL) || (!p_IsConstant(aAsPoly, naRing))) return FALSE;
353  return n_IsMOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
354}
355
356/// this is in-place, modifies a
357number naNeg(number a, const coeffs cf)
358{
359  naTest(a);
360  if (a != NULL) a = (number)p_Neg((poly)a, naRing);
361  return a;
362}
363
364number naImPart(number a, const coeffs cf)
365{
366  naTest(a);
367  return NULL;
368}
369
370number naInit_bigint(number longratBigIntNumber, const coeffs src, const coeffs cf)
371{
372  assume( cf != NULL );
373
374  const ring A = cf->extRing;
375
376  assume( A != NULL );
377
378  const coeffs C = A->cf;
379
380  assume( C != NULL );
381
382  number n = n_Init_bigint(longratBigIntNumber, src, C);
383
384  if ( n_IsZero(n, C) )
385  {
386    n_Delete(&n, C);
387    return NULL;
388  }
389
390  return (number)p_NSet(n, A);
391}
392
393
394
395number naInit(long i, const coeffs cf)
396{
397  if (i == 0) return NULL;
398  else        return (number)p_ISet(i, naRing);
399}
400
401int naInt(number &a, const coeffs cf)
402{
403  naTest(a);
404  poly aAsPoly = (poly)a;
405  if(aAsPoly == NULL)
406    return 0;
407  if (!p_IsConstant(aAsPoly, naRing))
408    return 0;
409  assume( aAsPoly != NULL );
410  return n_Int(p_GetCoeff(aAsPoly, naRing), naCoeffs);
411}
412
413/* TRUE iff (a != 0 and (b == 0 or deg(a) > deg(b) or (deg(a)==deg(b) && lc(a)>lc(b))) */
414BOOLEAN naGreater(number a, number b, const coeffs cf)
415{
416  naTest(a); naTest(b);
417  if (naIsZero(a, cf))
418  {
419    if (naIsZero(b, cf)) return FALSE;
420    return !n_GreaterZero(pGetCoeff((poly)b),cf);
421  }
422  if (naIsZero(b, cf))
423  {
424    return n_GreaterZero(pGetCoeff((poly)a),cf);
425  }
426  int aDeg = p_Totaldegree((poly)a, naRing);
427  int bDeg = p_Totaldegree((poly)b, naRing);
430  return n_Greater(pGetCoeff((poly)a),pGetCoeff((poly)b),cf);
431}
432
433/* TRUE iff a != 0 and (LC(a) > 0 or deg(a) > 0) */
434BOOLEAN naGreaterZero(number a, const coeffs cf)
435{
436  naTest(a);
437  if (a == NULL)                                            return FALSE;
438  if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE;
439  if (p_Totaldegree((poly)a, naRing) > 0)                   return TRUE;
440  return FALSE;
441}
442
443void naCoeffWrite(const coeffs cf, BOOLEAN details)
444{
445  assume( cf != NULL );
446
447  const ring A = cf->extRing;
448
449  assume( A != NULL );
450  assume( A->cf != NULL );
451
452  n_CoeffWrite(A->cf, details);
453
454//  rWrite(A);
455
456  const int P = rVar(A);
457  assume( P > 0 );
458
459  Print("//   %d parameter    : ", P);
460
461  for (int nop=0; nop < P; nop ++)
462    Print("%s ", rRingVar(nop, A));
463
464  PrintLn();
465
466  const ideal I = A->qideal;
467
468  assume( I != NULL );
469  assume( IDELEMS(I) == 1 );
470
471
472  if ( details )
473  {
474    PrintS("//   minpoly        : (");
475    p_Write0( I->m, A);
476    PrintS(")");
477  }
478  else
479    PrintS("//   minpoly        : ...");
480
481  PrintLn();
482
483/*
484  char *x = rRingVar(0, A);
485
486  Print("//   Coefficients live in the extension field K[%s]/<f(%s)>\n", x, x);
487  Print("//   with the minimal polynomial f(%s) = %s\n", x,
488        p_String(A->qideal->m, A));
489  PrintS("//   and K: ");
490*/
491}
492
493number naAdd(number a, number b, const coeffs cf)
494{
495  naTest(a); naTest(b);
496  if (a == NULL) return naCopy(b, cf);
497  if (b == NULL) return naCopy(a, cf);
498  poly aPlusB = p_Add_q(p_Copy((poly)a, naRing),
499                        p_Copy((poly)b, naRing), naRing);
500  definiteReduce(aPlusB, naMinpoly, cf);
501  return (number)aPlusB;
502}
503
504number naSub(number a, number b, const coeffs cf)
505{
506  naTest(a); naTest(b);
507  if (b == NULL) return naCopy(a, cf);
508  poly minusB = p_Neg(p_Copy((poly)b, naRing), naRing);
509  if (a == NULL) return (number)minusB;
510  poly aMinusB = p_Add_q(p_Copy((poly)a, naRing), minusB, naRing);
511  definiteReduce(aMinusB, naMinpoly, cf);
512  return (number)aMinusB;
513}
514
515number naMult(number a, number b, const coeffs cf)
516{
517  naTest(a); naTest(b);
518  if (a == NULL) return NULL;
519  if (b == NULL) return NULL;
520  poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing),
521                          p_Copy((poly)b, naRing), naRing);
522  definiteReduce(aTimesB, naMinpoly, cf);
523  return (number)aTimesB;
524}
525
526number naDiv(number a, number b, const coeffs cf)
527{
528  naTest(a); naTest(b);
529  if (b == NULL) WerrorS(nDivBy0);
530  if (a == NULL) return NULL;
531  poly bInverse = (poly)naInvers(b, cf);
532  if(bInverse != NULL) // b is non-zero divisor!
533  {
534    poly aDivB = p_Mult_q(p_Copy((poly)a, naRing), bInverse, naRing);
537  }
538  return NULL;
539}
540
541/* 0^0 = 0;
542   for |exp| <= 7 compute power by a simple multiplication loop;
543   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
544      p^13 = p^1 * p^4 * p^8, where we utilise that
545      p^(2^(k+1)) = p^(2^k) * p^(2^k);
546   intermediate reduction modulo the minimal polynomial is controlled by
547   the in-place method heuristicReduce(poly, poly, coeffs); see there.
548*/
549void naPower(number a, int exp, number *b, const coeffs cf)
550{
551  naTest(a);
552
553  /* special cases first */
554  if (a == NULL)
555  {
556    if (exp >= 0) *b = NULL;
557    else          WerrorS(nDivBy0);
558    return;
559  }
560  else if (exp ==  0) { *b = naInit(1, cf); return; }
561  else if (exp ==  1) { *b = naCopy(a, cf); return; }
562  else if (exp == -1) { *b = naInvers(a, cf); return; }
563
564  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
565
566  /* now compute a^expAbs */
567  poly pow; poly aAsPoly = (poly)a;
568  if (expAbs <= 7)
569  {
570    pow = p_Copy(aAsPoly, naRing);
571    for (int i = 2; i <= expAbs; i++)
572    {
573      pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing);
574      heuristicReduce(pow, naMinpoly, cf);
575    }
576    definiteReduce(pow, naMinpoly, cf);
577  }
578  else
579  {
580    pow = p_ISet(1, naRing);
581    poly factor = p_Copy(aAsPoly, naRing);
582    while (expAbs != 0)
583    {
584      if (expAbs & 1)
585      {
586        pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing);
587        heuristicReduce(pow, naMinpoly, cf);
588      }
589      expAbs = expAbs / 2;
590      if (expAbs != 0)
591      {
592        factor = p_Mult_q(factor, p_Copy(factor, naRing), naRing);
593        heuristicReduce(factor, naMinpoly, cf);
594      }
595    }
596    p_Delete(&factor, naRing);
597    definiteReduce(pow, naMinpoly, cf);
598  }
599
600  /* invert if original exponent was negative */
601  number n = (number)pow;
602  if (exp < 0)
603  {
604    number m = naInvers(n, cf);
606    n = m;
607  }
608  *b = n;
609}
610
611/* may reduce p modulo the reducer by calling definiteReduce;
612   the decision is made based on the following heuristic
613   (which should also only be changed here in this method):
614      if (deg(p) > 10*deg(reducer) then perform reduction;
615   modifies p */
616void heuristicReduce(poly &p, poly reducer, const coeffs cf)
617{
618  #ifdef LDEBUG
619  p_Test((poly)p, naRing);
620  p_Test((poly)reducer, naRing);
621  #endif
622  if (p_Totaldegree(p, naRing) > 10 * p_Totaldegree(reducer, naRing))
623    definiteReduce(p, reducer, cf);
624}
625
626void naWriteLong(number &a, const coeffs cf)
627{
628  naTest(a);
629  if (a == NULL)
630    StringAppendS("0");
631  else
632  {
633    poly aAsPoly = (poly)a;
634    /* basically, just write aAsPoly using p_Write,
635       but use brackets around the output, if a is not
636       a constant living in naCoeffs = cf->extRing->cf */
637    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
638    if (useBrackets) StringAppendS("(");
639    p_String0Long(aAsPoly, naRing, naRing);
640    if (useBrackets) StringAppendS(")");
641  }
642}
643
644void naWriteShort(number &a, const coeffs cf)
645{
646  naTest(a);
647  if (a == NULL)
648    StringAppendS("0");
649  else
650  {
651    poly aAsPoly = (poly)a;
652    /* basically, just write aAsPoly using p_Write,
653       but use brackets around the output, if a is not
654       a constant living in naCoeffs = cf->extRing->cf */
655    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
656    if (useBrackets) StringAppendS("(");
657    p_String0Short(aAsPoly, naRing, naRing);
658    if (useBrackets) StringAppendS(")");
659  }
660}
661
662const char * naRead(const char *s, number *a, const coeffs cf)
663{
664  poly aAsPoly;
665  const char * result = p_Read(s, aAsPoly, naRing);
666  definiteReduce(aAsPoly, naMinpoly, cf);
667  *a = (number)aAsPoly;
668  return result;
669}
670
671#if 0
672/* implemented by the rule lcm(a, b) = a * b / gcd(a, b) */
673number naLcm(number a, number b, const coeffs cf)
674{
675  naTest(a); naTest(b);
676  if (a == NULL) return NULL;
677  if (b == NULL) return NULL;
678  number theProduct = (number)p_Mult_q(p_Copy((poly)a, naRing),
679                                       p_Copy((poly)b, naRing), naRing);
680  /* note that theProduct needs not be reduced w.r.t. naMinpoly;
681     but the final division will take care of the necessary reduction */
682  number theGcd = naGcd(a, b, cf);
684}
685#endif
686
687/* expects *param to be castable to AlgExtInfo */
688static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
689{
690  if (ID != n) return FALSE;
691  AlgExtInfo *e = (AlgExtInfo *)param;
692  /* for extension coefficient fields we expect the underlying
693     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
694     this expectation is based on the assumption that we have properly
695     registered cf and perform reference counting rather than creating
696     multiple copies of the same coefficient field/domain/ring */
697  if (naRing == e->r)
698    return TRUE;
699  /* (Note that then also the minimal ideals will necessarily be
700     the same, as they are attached to the ring.) */
701
702  // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)...
703  if( rEqual(naRing, e->r, TRUE) ) // also checks the equality of qideals
704  {
705    const ideal mi = naRing->qideal;
706    assume( IDELEMS(mi) == 1 );
707    const ideal ii = e->r->qideal;
708    assume( IDELEMS(ii) == 1 );
709
710    // TODO: the following should be extended for 2 *equal* rings...
711    assume( p_EqualPolys(mi->m, ii->m, naRing, e->r) );
712
713    rDelete(e->r);
714
715    return TRUE;
716  }
717
718  return FALSE;
719
720}
721
722int naSize(number a, const coeffs cf)
723{
724  if (a == NULL) return -1;
725  /* this has been taken from the old implementation of field extensions,
726     where we computed the sum of the degree and the number of terms in
727     (poly)a; so we leave it at that, for the time being;
728     maybe, the number of terms alone is a better measure? */
729  poly aAsPoly = (poly)a;
730  int theDegree = 0; int noOfTerms = 0;
731  while (aAsPoly != NULL)
732  {
733    noOfTerms++;
734    int d = p_GetExp(aAsPoly, 1, naRing);
735    if (d > theDegree) theDegree = d;
736    pIter(aAsPoly);
737  }
738  return theDegree + noOfTerms;
739}
740
741/* performs polynomial division and overrides p by the remainder
742   of division of p by the reducer;
743   modifies p */
744void definiteReduce(poly &p, poly reducer, const coeffs cf)
745{
746  #ifdef LDEBUG
747  p_Test((poly)p, naRing);
748  p_Test((poly)reducer, naRing);
749  #endif
750  if ((p!=NULL) && (p_GetExp(p,1,naRing)>=p_GetExp(reducer,1,naRing)))
751  {
752    p_PolyDiv(p, reducer, FALSE, naRing);
753  }
754}
755
756void  naNormalize(number &a, const coeffs cf)
757{
758  poly aa=(poly)a;
759  if (aa!=naMinpoly)
760    definiteReduce(aa,naMinpoly,cf);
761  a=(number)aa;
762}
763
764#ifdef HAVE_FACTORY
765number naConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
766{
767  if (n.isZero()) return NULL;
768  poly p=convFactoryPSingP(n,naRing);
769  return (number)p;
770}
771CanonicalForm naConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
772{
773  naTest(n);
774  if (n==NULL) return CanonicalForm(0);
775
776  return convSingPFactoryP((poly)n,naRing);
777}
778#endif
779
780
781/* IMPORTANT NOTE: Since an algebraic field extension is again a field,
782                   the gcd of two elements is not very interesting. (It
783                   is actually any unit in the field, i.e., any non-
784                   zero element.) Note that the below method does not operate
785                   in this strong sense but rather computes the gcd of
786                   two given elements in the underlying polynomial ring. */
787number naGcd(number a, number b, const coeffs cf)
788{
789  naTest(a); naTest(b);
790  if ((a == NULL) && (b == NULL)) WerrorS(nDivBy0);
791  const ring R = naRing;
792  return (number) singclap_gcd(p_Copy((poly)a, R), p_Copy((poly)b, R), R);
793//  return (number)p_Gcd((poly)a, (poly)b, naRing);
794}
795
796number naInvers(number a, const coeffs cf)
797{
798  naTest(a);
799  if (a == NULL) WerrorS(nDivBy0);
800
801  poly aFactor = NULL; poly mFactor = NULL; poly theGcd = NULL;
802// singclap_extgcd!
803  const BOOLEAN ret = singclap_extgcd ((poly)a, naMinpoly, theGcd, aFactor, mFactor, naRing);
804
805  assume( !ret );
806
807//  if( ret ) theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
808
809  naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor);
810  p_Delete(&mFactor, naRing);
811
812  //  /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */
813  //  assume(naIsOne((number)theGcd, cf));
814
815  if( !naIsOne((number)theGcd, cf) )
816  {
817    WerrorS("zero divisor found - your minpoly is not irreducible");
818    p_Delete(&aFactor, naRing); aFactor = NULL;
819  }
820  p_Delete(&theGcd, naRing);
821
822  return (number)(aFactor);
823}
824
825/* assumes that src = Q, dst = Q(a) */
826number naMap00(number a, const coeffs src, const coeffs dst)
827{
828  if (n_IsZero(a, src)) return NULL;
829  assume(src == dst->extRing->cf);
830  poly result = p_One(dst->extRing);
831  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
832  return (number)result;
833}
834
835/* assumes that src = Z/p, dst = Q(a) */
836number naMapP0(number a, const coeffs src, const coeffs dst)
837{
838  if (n_IsZero(a, src)) return NULL;
839  /* mapping via intermediate int: */
840  int n = n_Int(a, src);
841  number q = n_Init(n, dst->extRing->cf);
842  poly result = p_One(dst->extRing);
843  p_SetCoeff(result, q, dst->extRing);
844  return (number)result;
845}
846
847#if 0
848/* assumes that either src = Q(a), dst = Q(a), or
849                       src = Z/p(a), dst = Z/p(a)     */
850number naCopyMap(number a, const coeffs src, const coeffs dst)
851{
852  return naCopy(a, dst);
853}
854#endif
855
856number naCopyExt(number a, const coeffs src, const coeffs)
857{
858  fraction fa=(fraction)a;
859  return (number)p_Copy(NUM(fa),src->extRing);
860}
861
862/* assumes that src = Q, dst = Z/p(a) */
863number naMap0P(number a, const coeffs src, const coeffs dst)
864{
865  if (n_IsZero(a, src)) return NULL;
866  int p = rChar(dst->extRing);
867
868  number q = nlModP(a, src, dst->extRing->cf);
869
870  poly result = p_NSet(q, dst->extRing);
871
872  return (number)result;
873}
874
875/* assumes that src = Z/p, dst = Z/p(a) */
876number naMapPP(number a, const coeffs src, const coeffs dst)
877{
878  if (n_IsZero(a, src)) return NULL;
879  assume(src == dst->extRing->cf);
880  poly result = p_One(dst->extRing);
881  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
882  return (number)result;
883}
884
885/* assumes that src = Z/u, dst = Z/p(a), where u != p */
886number naMapUP(number a, const coeffs src, const coeffs dst)
887{
888  if (n_IsZero(a, src)) return NULL;
889  /* mapping via intermediate int: */
890  int n = n_Int(a, src);
891  number q = n_Init(n, dst->extRing->cf);
892  poly result = p_One(dst->extRing);
893  p_SetCoeff(result, q, dst->extRing);
894  return (number)result;
895}
896
897nMapFunc naSetMap(const coeffs src, const coeffs dst)
898{
899  /* dst is expected to be an algebraic field extension */
900  assume(getCoeffType(dst) == ID);
901
902  if( src == dst ) return ndCopyMap;
903
904  int h = 0; /* the height of the extension tower given by dst */
905  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
906  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
907
908  /* for the time being, we only provide maps if h = 1 and if b is Q or
909     some field Z/pZ: */
910  if (h==0)
911  {
912    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
913      return naMap00;                            /// Q     -->  Q(a)
914    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
915      return naMapP0;                            /// Z/p   -->  Q(a)
916    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
917      return naMap0P;                            /// Q      --> Z/p(a)
918    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
919    {
920      if (src->ch == dst->ch) return naMapPP;    /// Z/p    --> Z/p(a)
921      else return naMapUP;                       /// Z/u    --> Z/p(a)
922    }
923  }
924  if (h != 1) return NULL;
925  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
926  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
927
928  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
929  {
930    if (strcmp(rRingVar(0, src->extRing),
931               rRingVar(0, dst->extRing)) == 0)
932    {
933      if (src->type==n_algExt)
934         return ndCopyMap; // naCopyMap;         /// Q(a)   --> Q(a)
935      else
936         return naCopyExt;
937    }
938    else
939      return NULL;                               /// Q(b)   --> Q(a)
940  }
941
942  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
943  {
944    if (strcmp(rRingVar(0,src->extRing),rRingVar(0,dst->extRing))==0)
945    {
946      if (src->type==n_algExt)
947        return ndCopyMap; // naCopyMap;          /// Z/p(a) --> Z/p(a)
948      else
949         return naCopyExt;
950    }
951    else
952      return NULL;                               /// Z/p(b) --> Z/p(a)
953  }
954
955  return NULL;                                           /// default
956}
957
958static int naParDeg(number a, const coeffs cf)
959{
960  if (a == NULL) return -1;
961  poly aa=(poly)a;
962  return cf->extRing->pFDeg(aa,cf->extRing);
963}
964
965/// return the specified parameter as a number in the given alg. field
966static number naParameter(const int iParameter, const coeffs cf)
967{
968  assume(getCoeffType(cf) == ID);
969
970  const ring R = cf->extRing;
971  assume( R != NULL );
972  assume( 0 < iParameter && iParameter <= rVar(R) );
973
974  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
975
976  return (number) p;
977}
978
979
980/// if m == var(i)/1 => return i,
981int naIsParam(number m, const coeffs cf)
982{
983  assume(getCoeffType(cf) == ID);
984
985  const ring R = cf->extRing;
986  assume( R != NULL );
987
988  return p_Var( (poly)m, R );
989}
990
991
992static void naClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
993{
994  assume(cf != NULL);
995  assume(getCoeffType(cf) == ID);
996  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
997
998  const ring   R = cf->extRing;
999  assume(R != NULL);
1000  const coeffs Q = R->cf;
1001  assume(Q != NULL);
1002  assume(nCoeff_is_Q(Q));
1003
1004  numberCollectionEnumerator.Reset();
1005
1006  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
1007  {
1008    c = n_Init(1, cf);
1009    return;
1010  }
1011
1012  naTest(numberCollectionEnumerator.Current());
1013
1014  // part 1, find a small candidate for gcd
1015  int s1; int s=2147483647; // max. int
1016
1017  const BOOLEAN lc_is_pos=naGreaterZero(numberCollectionEnumerator.Current(),cf);
1018
1019  int normalcount = 0;
1020
1021  poly cand1, cand;
1022
1023  do
1024  {
1025    number& n = numberCollectionEnumerator.Current();
1026    naNormalize(n, cf); ++normalcount;
1027
1028    naTest(n);
1029
1030    cand1 = (poly)n;
1031
1032    s1 = p_Deg(cand1, R); // naSize?
1033    if (s>s1)
1034    {
1035      cand = cand1;
1036      s = s1;
1037    }
1038  } while (numberCollectionEnumerator.MoveNext() );
1039
1040//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
1041
1042  cand = p_Copy(cand, R);
1043  // part 2: compute gcd(cand,all coeffs)
1044
1045  numberCollectionEnumerator.Reset();
1046
1047  int length = 0;
1048  while (numberCollectionEnumerator.MoveNext() )
1049  {
1050    number& n = numberCollectionEnumerator.Current();
1051    ++length;
1052
1053    if( (--normalcount) <= 0)
1054      naNormalize(n, cf);
1055
1056    naTest(n);
1057
1058//    p_InpGcd(cand, (poly)n, R);
1059
1060    cand = singclap_gcd(cand, p_Copy((poly)n, R), R);
1061
1062//    cand1 = p_Gcd(cand,(poly)n, R); p_Delete(&cand, R); cand = cand1;
1063
1064    assume( naGreaterZero((number)cand, cf) ); // ???
1065/*
1066    if(p_IsConstant(cand,R))
1067    {
1068      c = cand;
1069
1070      if(!lc_is_pos)
1071      {
1072        // make the leading coeff positive
1073        c = nlNeg(c, cf);
1074        numberCollectionEnumerator.Reset();
1075
1076        while (numberCollectionEnumerator.MoveNext() )
1077        {
1078          number& nn = numberCollectionEnumerator.Current();
1079          nn = nlNeg(nn, cf);
1080        }
1081      }
1082      return;
1083    }
1084*/
1085
1086  }
1087
1088  // part3: all coeffs = all coeffs / cand
1089  if (!lc_is_pos)
1090    cand = p_Neg(cand, R);
1091
1092  c = (number)cand; naTest(c);
1093
1094  poly cInverse = (poly)naInvers(c, cf);
1095  assume(cInverse != NULL); // c is non-zero divisor!?
1096
1097
1098  numberCollectionEnumerator.Reset();
1099
1100
1101  while (numberCollectionEnumerator.MoveNext() )
1102  {
1103    number& n = numberCollectionEnumerator.Current();
1104
1105    assume( length > 0 );
1106
1107    if( --length > 0 )
1108    {
1109      assume( cInverse != NULL );
1110      n = (number) p_Mult_q(p_Copy(cInverse, R), (poly)n, R);
1111    }
1112    else
1113    {
1114      n = (number) p_Mult_q(cInverse, (poly)n, R);
1115      cInverse = NULL;
1116      assume(length == 0);
1117    }
1118
1119    definiteReduce((poly &)n, naMinpoly, cf);
1120  }
1121
1122  assume(length == 0);
1123  assume(cInverse == NULL); //   p_Delete(&cInverse, R);
1124
1125  // Quick and dirty fix for constant content clearing... !?
1126  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
1127
1128  number cc;
1129
1130  extern void nlClearContentNoPositiveLead(ICoeffsEnumerator&, number&, const coeffs);
1131
1132  nlClearContentNoPositiveLead(itr, cc, Q); // TODO: get rid of (-LC) normalization!?
1133
1134  // over alg. ext. of Q // takes over the input number
1135  c = (number) p_Mult_nn( (poly)c, cc, R);
1136//      p_Mult_q(p_NSet(cc, R), , R);
1137
1138  n_Delete(&cc, Q);
1139
1140  // TODO: the above is not enough! need GCD's of polynomial coeffs...!
1141/*
1142  // old and wrong part of p_Content
1143    if (rField_is_Q_a(r) && !CLEARENUMERATORS) // should not be used anymore if CLEARENUMERATORS is 1
1144    {
1145      // we only need special handling for alg. ext.
1146      if (getCoeffType(r->cf)==n_algExt)
1147      {
1148        number hzz = n_Init(1, r->cf->extRing->cf);
1149        p=ph;
1150        while (p!=NULL)
1151        { // each monom: coeff in Q_a
1152          poly c_n_n=(poly)pGetCoeff(p);
1153          poly c_n=c_n_n;
1154          while (c_n!=NULL)
1155          { // each monom: coeff in Q
1156            d=n_Lcm(hzz,pGetCoeff(c_n),r->cf->extRing->cf);
1157            n_Delete(&hzz,r->cf->extRing->cf);
1158            hzz=d;
1159            pIter(c_n);
1160          }
1161          pIter(p);
1162        }
1163        // hzz contains the 1/lcm of all denominators in c_n_n
1164        h=n_Invers(hzz,r->cf->extRing->cf);
1165        n_Delete(&hzz,r->cf->extRing->cf);
1166        n_Normalize(h,r->cf->extRing->cf);
1167        if(!n_IsOne(h,r->cf->extRing->cf))
1168        {
1169          p=ph;
1170          while (p!=NULL)
1171          { // each monom: coeff in Q_a
1172            poly c_n=(poly)pGetCoeff(p);
1173            while (c_n!=NULL)
1174            { // each monom: coeff in Q
1175              d=n_Mult(h,pGetCoeff(c_n),r->cf->extRing->cf);
1176              n_Normalize(d,r->cf->extRing->cf);
1177              n_Delete(&pGetCoeff(c_n),r->cf->extRing->cf);
1178              pGetCoeff(c_n)=d;
1179              pIter(c_n);
1180            }
1181            pIter(p);
1182          }
1183        }
1184        n_Delete(&h,r->cf->extRing->cf);
1185      }
1186    }
1187*/
1188
1189
1190//  c = n_Init(1, cf); assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
1191}
1192
1193
1194static void naClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1195{
1196  assume(cf != NULL);
1197  assume(getCoeffType(cf) == ID);
1198  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
1199
1200  assume(cf->extRing != NULL);
1201  const coeffs Q = cf->extRing->cf;
1202  assume(Q != NULL);
1203  assume(nCoeff_is_Q(Q));
1204  number n;
1205  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
1206
1207  extern void nlClearDenominatorsNoPositiveLead(ICoeffsEnumerator&, number&, const coeffs);
1208
1209  nlClearDenominatorsNoPositiveLead(itr, n, Q); // this should probably be fine...
1210  c = (number)p_NSet(n, cf->extRing); // over alg. ext. of Q // takes over the input number
1211}
1212
1213
1214BOOLEAN naInitChar(coeffs cf, void * infoStruct)
1215{
1216  assume( infoStruct != NULL );
1217
1218  AlgExtInfo *e = (AlgExtInfo *)infoStruct;
1219  /// first check whether cf->extRing != NULL and delete old ring???
1220
1221  assume(e->r                     != NULL);      // extRing;
1222  assume(e->r->cf                 != NULL);      // extRing->cf;
1223
1224  assume((e->r->qideal            != NULL) &&    // minideal has one
1225         (IDELEMS(e->r->qideal)   == 1)    &&    // non-zero generator
1226         (e->r->qideal->m      != NULL)    ); // at m;
1227
1228  assume( cf != NULL );
1229  assume(getCoeffType(cf) == ID);                     // coeff type;
1230
1231  e->r->ref ++; // increase the ref.counter for the ground poly. ring!
1232  const ring R = e->r; // no copy!
1233  assume( R->qideal == e->r->qideal );
1234  cf->extRing  = R;
1235
1236  /* propagate characteristic up so that it becomes
1237     directly accessible in cf: */
1238  cf->ch = R->cf->ch;
1239
1240  #ifdef LDEBUG
1241  p_Test((poly)naMinpoly, naRing);
1242  #endif
1243
1244  cf->cfGreaterZero  = naGreaterZero;
1245  cf->cfGreater      = naGreater;
1246  cf->cfEqual        = naEqual;
1247  cf->cfIsZero       = naIsZero;
1248  cf->cfIsOne        = naIsOne;
1249  cf->cfIsMOne       = naIsMOne;
1250  cf->cfInit         = naInit;
1251  cf->cfInit_bigint  = naInit_bigint;
1252  cf->cfInt          = naInt;
1253  cf->cfNeg          = naNeg;
1255  cf->cfSub          = naSub;
1256  cf->cfMult         = naMult;
1259  cf->cfPower        = naPower;
1260  cf->cfCopy         = naCopy;
1261
1262  cf->cfWriteLong        = naWriteLong;
1263
1264  if( rCanShortOut(naRing) )
1265    cf->cfWriteShort = naWriteShort;
1266  else
1267    cf->cfWriteShort = naWriteLong;
1268
1271  cf->cfSetMap       = naSetMap;
1272  cf->cfGetDenom     = naGetDenom;
1273  cf->cfGetNumerator = naGetNumerator;
1274  cf->cfRePart       = naCopy;
1275  cf->cfImPart       = naImPart;
1276  cf->cfCoeffWrite   = naCoeffWrite;
1277  cf->cfNormalize    = naNormalize;
1278#ifdef LDEBUG
1280#endif
1281  cf->cfGcd          = naGcd;
1282  //cf->cfLcm          = naLcm;
1283  cf->cfSize         = naSize;
1284  cf->nCoeffIsEqual  = naCoeffIsEqual;
1285  cf->cfInvers       = naInvers;
1286  cf->cfIntDiv       = naDiv; // ???
1287#ifdef HAVE_FACTORY
1288  cf->convFactoryNSingN=naConvFactoryNSingN;
1289  cf->convSingNFactoryN=naConvSingNFactoryN;
1290#endif
1291  cf->cfParDeg = naParDeg;
1292
1293  cf->iNumberOfParameters = rVar(R);
1294  cf->pParameterNames = R->names;
1295  cf->cfParameter = naParameter;
1296
1297  if( nCoeff_is_Q(R->cf) )
1298  {
1299    cf->cfClearContent = naClearContent;
1300    cf->cfClearDenominators = naClearDenominators;
1301  }
1302
1303  return FALSE;
1304}
Note: See TracBrowser for help on using the repository browser.