# source:git/libpolys/polys/ext_fields/algext.cc@7e9f12

spielwiese
Last change on this file since 7e9f12 was 7e9f12, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
removed minimal ideal from AlgExtInfo infavour of qring chg: algext. expects a qring (qideal != NULL) and uses it directly
• Property mode set to `100644`
File size: 25.3 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[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#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#ifdef HAVE_FACTORY
47#include <polys/clapconv.h>
48#include <factory/factory.h>
49#endif
50
51#include "ext_fields/algext.h"
52#define TRANSEXT_PRIVATES 1
53#include "ext_fields/transext.h"
54
55#ifdef LDEBUG
57BOOLEAN  naDBTest(number a, const char *f, const int l, const coeffs r);
58#else
59#define naTest(a) (TRUE)
60#endif
61
62/// Our own type!
63static const n_coeffType ID = n_algExt;
64
65/* polynomial ring in which our numbers live */
66#define naRing cf->extRing
67
68/* coeffs object in which the coefficients of our numbers live;
69 * methods attached to naCoeffs may be used to compute with the
71 * coefficients of our numbers */
72#define naCoeffs cf->extRing->cf
73
74/* minimal polynomial */
75#define naMinpoly naRing->qideal->m[0]
76
77/// forward declarations
78BOOLEAN  naGreaterZero(number a, const coeffs cf);
79BOOLEAN  naGreater(number a, number b, const coeffs cf);
80BOOLEAN  naEqual(number a, number b, const coeffs cf);
81BOOLEAN  naIsOne(number a, const coeffs cf);
82BOOLEAN  naIsMOne(number a, const coeffs cf);
83BOOLEAN  naIsZero(number a, const coeffs cf);
84number   naInit(long i, const coeffs cf);
85int      naInt(number &a, const coeffs cf);
86number   naNeg(number a, const coeffs cf);
87number   naInvers(number a, const coeffs cf);
88number   naAdd(number a, number b, const coeffs cf);
89number   naSub(number a, number b, const coeffs cf);
90number   naMult(number a, number b, const coeffs cf);
91number   naDiv(number a, number b, const coeffs cf);
92void     naPower(number a, int exp, number *b, const coeffs cf);
93number   naCopy(number a, const coeffs cf);
94void     naWriteLong(number &a, const coeffs cf);
95void     naWriteShort(number &a, const coeffs cf);
96number   naRePart(number a, const coeffs cf);
97number   naImPart(number a, const coeffs cf);
98number   naGetDenom(number &a, const coeffs cf);
99number   naGetNumerator(number &a, const coeffs cf);
100number   naGcd(number a, number b, const coeffs cf);
101//number   naLcm(number a, number b, const coeffs cf);
102int      naSize(number a, const coeffs cf);
103void     naDelete(number *a, const coeffs cf);
104void     naCoeffWrite(const coeffs cf, BOOLEAN details);
105number   naIntDiv(number a, number b, const coeffs cf);
106const char * naRead(const char *s, number *a, const coeffs cf);
107
108static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param);
109
110#ifdef LDEBUG
111BOOLEAN naDBTest(number a, const char *f, const int l, const coeffs cf)
112{
113  assume(getCoeffType(cf) == ID);
114  if (a == NULL) return TRUE;
115  p_Test((poly)a, naRing);
116  if((((poly)a)!=naMinpoly)
117  && p_Totaldegree((poly)a, naRing) >= p_Totaldegree(naMinpoly, naRing))
118  {
119    Print("deg >= deg(minpoly) in %s:%d\n",f,l);
120    return FALSE;
121  }
122  return TRUE;
123}
124#endif
125
126void heuristicReduce(poly &p, poly reducer, const coeffs cf);
127void definiteReduce(poly &p, poly reducer, const coeffs cf);
128
129/* returns the bottom field in this field extension tower; if the tower
130   is flat, i.e., if there is no extension, then r itself is returned;
131   as a side-effect, the counter 'height' is filled with the height of
132   the extension tower (in case the tower is flat, 'height' is zero) */
133static coeffs nCoeff_bottom(const coeffs r, int &height)
134{
135  assume(r != NULL);
136  coeffs cf = r;
137  height = 0;
138  while (nCoeff_is_Extension(cf))
139  {
140    assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL);
141    cf = cf->extRing->cf;
142    height++;
143  }
144  return cf;
145}
146
147BOOLEAN naIsZero(number a, const coeffs cf)
148{
149  naTest(a);
150  return (a == NULL);
151}
152
153void naDelete(number * a, const coeffs cf)
154{
155  if (*a == NULL) return;
156  if (((poly)*a)==naMinpoly) { *a=NULL;return;}
157  poly aAsPoly = (poly)(*a);
158  p_Delete(&aAsPoly, naRing);
159  *a = NULL;
160}
161
162BOOLEAN naEqual(number a, number b, const coeffs cf)
163{
164  naTest(a); naTest(b);
165
166  /// simple tests
167  if (a == b) return TRUE;
168  if ((a == NULL) && (b != NULL)) return FALSE;
169  if ((b == NULL) && (a != NULL)) return FALSE;
170
171  /// deg test
173  if (a != NULL) aDeg = p_Totaldegree((poly)a, naRing);
174  int bDeg = 0;
175  if (b != NULL) bDeg = p_Totaldegree((poly)b, naRing);
176  if (aDeg != bDeg) return FALSE;
177
178  /// subtraction test
179  number c = naSub(a, b, cf);
180  BOOLEAN result = naIsZero(c, cf);
182  return result;
183}
184
185number naCopy(number a, const coeffs cf)
186{
187  naTest(a);
188  if (a == NULL) return NULL;
189  if (((poly)a)==naMinpoly) return a;
190  return (number)p_Copy((poly)a, naRing);
191}
192
193number naGetNumerator(number &a, const coeffs cf)
194{
195  return naCopy(a, cf);
196}
197
198number naGetDenom(number &a, const coeffs cf)
199{
200  naTest(a);
201  return naInit(1, cf);
202}
203
204BOOLEAN naIsOne(number a, const coeffs cf)
205{
206  naTest(a);
207  poly aAsPoly = (poly)a;
208  if (!p_IsConstant(aAsPoly, naRing)) return FALSE;
209  return n_IsOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
210}
211
212BOOLEAN naIsMOne(number a, const coeffs cf)
213{
214  naTest(a);
215  poly aAsPoly = (poly)a;
216  if (!p_IsConstant(aAsPoly, naRing)) return FALSE;
217  return n_IsMOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
218}
219
220/// this is in-place, modifies a
221number naNeg(number a, const coeffs cf)
222{
223  naTest(a);
224  if (a != NULL) a = (number)p_Neg((poly)a, naRing);
225  return a;
226}
227
228number naImPart(number a, const coeffs cf)
229{
230  naTest(a);
231  return NULL;
232}
233
234number naInit_bigint(number longratBigIntNumber, const coeffs src, const coeffs cf)
235{
236  assume( cf != NULL );
237
238  const ring A = cf->extRing;
239
240  assume( A != NULL );
241
242  const coeffs C = A->cf;
243
244  assume( C != NULL );
245
246  number n = n_Init_bigint(longratBigIntNumber, src, C);
247
248  if ( n_IsZero(n, C) )
249  {
250    n_Delete(&n, C);
251    return NULL;
252  }
253
254  return (number)p_NSet(n, A);
255}
256
257
258
259number naInit(long i, const coeffs cf)
260{
261  if (i == 0) return NULL;
262  else        return (number)p_ISet(i, naRing);
263}
264
265int naInt(number &a, const coeffs cf)
266{
267  naTest(a);
268  poly aAsPoly = (poly)a;
269  if(aAsPoly == NULL)
270    return 0;
271  if (!p_IsConstant(aAsPoly, naRing))
272    return 0;
273  assume( aAsPoly != NULL );
274  return n_Int(p_GetCoeff(aAsPoly, naRing), naCoeffs);
275}
276
277/* TRUE iff (a != 0 and (b == 0 or deg(a) > deg(b))) */
278BOOLEAN naGreater(number a, number b, const coeffs cf)
279{
280  naTest(a); naTest(b);
281  if (naIsZero(a, cf)) return FALSE;
282  if (naIsZero(b, cf)) return TRUE;
284  if (a != NULL) aDeg = p_Totaldegree((poly)a, naRing);
285  int bDeg = 0;
286  if (b != NULL) bDeg = p_Totaldegree((poly)b, naRing);
288}
289
290/* TRUE iff a != 0 and (LC(a) > 0 or deg(a) > 0) */
291BOOLEAN naGreaterZero(number a, const coeffs cf)
292{
293  naTest(a);
294  if (a == NULL)                                            return FALSE;
295  if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE;
296  if (p_Totaldegree((poly)a, naRing) > 0)                   return TRUE;
297  return FALSE;
298}
299
300void naCoeffWrite(const coeffs cf, BOOLEAN details)
301{
302  assume( cf != NULL );
303
304  const ring A = cf->extRing;
305
306  assume( A != NULL );
307  assume( A->cf != NULL );
308
309  n_CoeffWrite(A->cf, details);
310
311//  rWrite(A);
312
313  const int P = rVar(A);
314  assume( P > 0 );
315
316  Print("//   %d parameter    : ", P);
317
318  for (int nop=0; nop < P; nop ++)
319    Print("%s ", rRingVar(nop, A));
320
321  PrintLn();
322
323  const ideal I = A->qideal;
324
325  assume( I != NULL );
326  assume( IDELEMS(I) == 1 );
327
328
329  if ( details )
330  {
331    PrintS("//   minpoly        : (");
332    p_Write0( I->m[0], A);
333    PrintS(")");
334  }
335  else
336    PrintS("//   minpoly        : ...");
337
338  PrintLn();
339
340/*
341  char *x = rRingVar(0, A);
342
343  Print("//   Coefficients live in the extension field K[%s]/<f(%s)>\n", x, x);
344  Print("//   with the minimal polynomial f(%s) = %s\n", x,
345        p_String(A->qideal->m[0], A));
346  PrintS("//   and K: ");
347*/
348}
349
350number naAdd(number a, number b, const coeffs cf)
351{
352  naTest(a); naTest(b);
353  if (a == NULL) return naCopy(b, cf);
354  if (b == NULL) return naCopy(a, cf);
355  poly aPlusB = p_Add_q(p_Copy((poly)a, naRing),
356                        p_Copy((poly)b, naRing), naRing);
357  definiteReduce(aPlusB, naMinpoly, cf);
358  return (number)aPlusB;
359}
360
361number naSub(number a, number b, const coeffs cf)
362{
363  naTest(a); naTest(b);
364  if (b == NULL) return naCopy(a, cf);
365  poly minusB = p_Neg(p_Copy((poly)b, naRing), naRing);
366  if (a == NULL) return (number)minusB;
367  poly aMinusB = p_Add_q(p_Copy((poly)a, naRing), minusB, naRing);
368  definiteReduce(aMinusB, naMinpoly, cf);
369  return (number)aMinusB;
370}
371
372number naMult(number a, number b, const coeffs cf)
373{
374  naTest(a); naTest(b);
375  if (a == NULL) return NULL;
376  if (b == NULL) return NULL;
377  poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing),
378                          p_Copy((poly)b, naRing), naRing);
379  definiteReduce(aTimesB, naMinpoly, cf);
380  return (number)aTimesB;
381}
382
383number naDiv(number a, number b, const coeffs cf)
384{
385  naTest(a); naTest(b);
386  if (b == NULL) WerrorS(nDivBy0);
387  if (a == NULL) return NULL;
388  poly bInverse = (poly)naInvers(b, cf);
389  poly aDivB = p_Mult_q(p_Copy((poly)a, naRing), bInverse, naRing);
392}
393
394/* 0^0 = 0;
395   for |exp| <= 7 compute power by a simple multiplication loop;
396   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
397      p^13 = p^1 * p^4 * p^8, where we utilise that
398      p^(2^(k+1)) = p^(2^k) * p^(2^k);
399   intermediate reduction modulo the minimal polynomial is controlled by
400   the in-place method heuristicReduce(poly, poly, coeffs); see there.
401*/
402void naPower(number a, int exp, number *b, const coeffs cf)
403{
404  naTest(a);
405
406  /* special cases first */
407  if (a == NULL)
408  {
409    if (exp >= 0) *b = NULL;
410    else          WerrorS(nDivBy0);
411  }
412  else if (exp ==  0) { *b = naInit(1, cf); return; }
413  else if (exp ==  1) { *b = naCopy(a, cf); return; }
414  else if (exp == -1) { *b = naInvers(a, cf); return; }
415
416  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
417
418  /* now compute a^expAbs */
419  poly pow; poly aAsPoly = (poly)a;
420  if (expAbs <= 7)
421  {
422    pow = p_Copy(aAsPoly, naRing);
423    for (int i = 2; i <= expAbs; i++)
424    {
425      pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing);
426      heuristicReduce(pow, naMinpoly, cf);
427    }
428    definiteReduce(pow, naMinpoly, cf);
429  }
430  else
431  {
432    pow = p_ISet(1, naRing);
433    poly factor = p_Copy(aAsPoly, naRing);
434    while (expAbs != 0)
435    {
436      if (expAbs & 1)
437      {
438        pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing);
439        heuristicReduce(pow, naMinpoly, cf);
440      }
441      expAbs = expAbs / 2;
442      if (expAbs != 0)
443      {
444        factor = p_Mult_q(factor, p_Copy(factor, naRing), naRing);
445        heuristicReduce(factor, naMinpoly, cf);
446      }
447    }
448    p_Delete(&factor, naRing);
449    definiteReduce(pow, naMinpoly, cf);
450  }
451
452  /* invert if original exponent was negative */
453  number n = (number)pow;
454  if (exp < 0)
455  {
456    number m = naInvers(n, cf);
458    n = m;
459  }
460  *b = n;
461}
462
463/* may reduce p modulo the reducer by calling definiteReduce;
464   the decision is made based on the following heuristic
465   (which should also only be changed here in this method):
466      if (deg(p) > 10*deg(reducer) then perform reduction;
467   modifies p */
468void heuristicReduce(poly &p, poly reducer, const coeffs cf)
469{
470  #ifdef LDEBUG
471  p_Test((poly)p, naRing);
472  p_Test((poly)reducer, naRing);
473  #endif
474  if (p_Totaldegree(p, naRing) > 10 * p_Totaldegree(reducer, naRing))
475    definiteReduce(p, reducer, cf);
476}
477
478void naWriteLong(number &a, const coeffs cf)
479{
480  naTest(a);
481  if (a == NULL)
482    StringAppendS("0");
483  else
484  {
485    poly aAsPoly = (poly)a;
486    /* basically, just write aAsPoly using p_Write,
487       but use brackets around the output, if a is not
488       a constant living in naCoeffs = cf->extRing->cf */
489    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
490    if (useBrackets) StringAppendS("(");
491    p_String0Long(aAsPoly, naRing, naRing);
492    if (useBrackets) StringAppendS(")");
493  }
494}
495
496void naWriteShort(number &a, const coeffs cf)
497{
498  naTest(a);
499  if (a == NULL)
500    StringAppendS("0");
501  else
502  {
503    poly aAsPoly = (poly)a;
504    /* basically, just write aAsPoly using p_Write,
505       but use brackets around the output, if a is not
506       a constant living in naCoeffs = cf->extRing->cf */
507    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
508    if (useBrackets) StringAppendS("(");
509    p_String0Short(aAsPoly, naRing, naRing);
510    if (useBrackets) StringAppendS(")");
511  }
512}
513
514const char * naRead(const char *s, number *a, const coeffs cf)
515{
516  poly aAsPoly;
517  const char * result = p_Read(s, aAsPoly, naRing);
518  definiteReduce(aAsPoly, naMinpoly, cf);
519  *a = (number)aAsPoly;
520  return result;
521}
522
523#if 0
524/* implemented by the rule lcm(a, b) = a * b / gcd(a, b) */
525number naLcm(number a, number b, const coeffs cf)
526{
527  naTest(a); naTest(b);
528  if (a == NULL) return NULL;
529  if (b == NULL) return NULL;
530  number theProduct = (number)p_Mult_q(p_Copy((poly)a, naRing),
531                                       p_Copy((poly)b, naRing), naRing);
532  /* note that theProduct needs not be reduced w.r.t. naMinpoly;
533     but the final division will take care of the necessary reduction */
534  number theGcd = naGcd(a, b, cf);
536}
537#endif
538
539/* expects *param to be castable to AlgExtInfo */
540static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
541{
542  if (ID != n) return FALSE;
543  AlgExtInfo *e = (AlgExtInfo *)param;
544  /* for extension coefficient fields we expect the underlying
545     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
546     this expectation is based on the assumption that we have properly
547     registered cf and perform reference counting rather than creating
548     multiple copies of the same coefficient field/domain/ring */
549  if (naRing == e->r)
550    return TRUE;
551  /* (Note that then also the minimal ideals will necessarily be
552     the same, as they are attached to the ring.) */
553
554  // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)...
555  if( rEqual(naRing, e->r, TRUE) ) // also checks the equality of qideals
556  {
557    const ideal mi = naRing->qideal;
558    assume( IDELEMS(mi) == 1 );
559    const ideal ii = e->r->qideal;
560    assume( IDELEMS(ii) == 1 );
561
562    // TODO: the following should be extended for 2 *equal* rings...
563    assume( p_EqualPolys(mi->m[0], ii->m[0], naRing, e->r) );
564
565    rDelete(e->r);
566
567    return TRUE;
568  }
569
570  return FALSE;
571
572}
573
574int naSize(number a, const coeffs cf)
575{
576  if (a == NULL) return -1;
577  /* this has been taken from the old implementation of field extensions,
578     where we computed the sum of the degree and the number of terms in
579     (poly)a; so we leave it at that, for the time being;
580     maybe, the number of terms alone is a better measure? */
581  poly aAsPoly = (poly)a;
582  int theDegree = 0; int noOfTerms = 0;
583  while (aAsPoly != NULL)
584  {
585    noOfTerms++;
586    int d = p_GetExp(aAsPoly, 1, naRing);
587    if (d > theDegree) theDegree = d;
588    pIter(aAsPoly);
589  }
590  return theDegree + noOfTerms;
591}
592
593/* performs polynomial division and overrides p by the remainder
594   of division of p by the reducer;
595   modifies p */
596void definiteReduce(poly &p, poly reducer, const coeffs cf)
597{
598  #ifdef LDEBUG
599  p_Test((poly)p, naRing);
600  p_Test((poly)reducer, naRing);
601  #endif
602  p_PolyDiv(p, reducer, FALSE, naRing);
603}
604
605void  naNormalize(number &a, const coeffs cf)
606{
607  poly aa=(poly)a;
608  if (aa!=naMinpoly)
609    definiteReduce(aa,naMinpoly,cf);
610  a=(number)aa;
611}
612
613#ifdef HAVE_FACTORY
614number naConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
615{
616  if (n.isZero()) return NULL;
617  poly p=convFactoryPSingP(n,naRing);
618  return (number)p;
619}
620CanonicalForm naConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
621{
622  naTest(n);
623  if (n==NULL) return CanonicalForm(0);
624
625  return convSingPFactoryP((poly)n,naRing);
626}
627#endif
628
629
630/* IMPORTANT NOTE: Since an algebraic field extension is again a field,
631                   the gcd of two elements is not very interesting. (It
632                   is actually any unit in the field, i.e., any non-
633                   zero element.) Note that the below method does not operate
634                   in this strong sense but rather computes the gcd of
635                   two given elements in the underlying polynomial ring. */
636number naGcd(number a, number b, const coeffs cf)
637{
638  naTest(a); naTest(b);
639  if ((a == NULL) && (b == NULL)) WerrorS(nDivBy0);
640  return (number)p_Gcd((poly)a, (poly)b, naRing);
641}
642
643number naInvers(number a, const coeffs cf)
644{
645  naTest(a);
646  if (a == NULL) WerrorS(nDivBy0);
647  poly aFactor = NULL; poly mFactor = NULL;
648  poly theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
649  naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor);
650  /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */
651  assume(naIsOne((number)theGcd, cf));
652  p_Delete(&theGcd, naRing);
653  p_Delete(&mFactor, naRing);
654  return (number)(aFactor);
655}
656
657/* assumes that src = Q, dst = Q(a) */
658number naMap00(number a, const coeffs src, const coeffs dst)
659{
660  if (n_IsZero(a, src)) return NULL;
661  assume(src == dst->extRing->cf);
662  poly result = p_One(dst->extRing);
663  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
664  return (number)result;
665}
666
667/* assumes that src = Z/p, dst = Q(a) */
668number naMapP0(number a, const coeffs src, const coeffs dst)
669{
670  if (n_IsZero(a, src)) return NULL;
671  /* mapping via intermediate int: */
672  int n = n_Int(a, src);
673  number q = n_Init(n, dst->extRing->cf);
674  poly result = p_One(dst->extRing);
675  p_SetCoeff(result, q, dst->extRing);
676  return (number)result;
677}
678
679#if 0
680/* assumes that either src = Q(a), dst = Q(a), or
681                       src = Z/p(a), dst = Z/p(a)     */
682number naCopyMap(number a, const coeffs src, const coeffs dst)
683{
684  return naCopy(a, dst);
685}
686#endif
687
688number naCopyExt(number a, const coeffs src, const coeffs)
689{
690  fraction fa=(fraction)a;
691  return (number)p_Copy(NUM(fa),src->extRing);
692}
693
694/* assumes that src = Q, dst = Z/p(a) */
695number naMap0P(number a, const coeffs src, const coeffs dst)
696{
697  if (n_IsZero(a, src)) return NULL;
698  int p = rChar(dst->extRing);
699
700  number q = nlModP(a, src, dst->extRing->cf);
701
702  poly result = p_NSet(q, dst->extRing);
703
704  return (number)result;
705}
706
707/* assumes that src = Z/p, dst = Z/p(a) */
708number naMapPP(number a, const coeffs src, const coeffs dst)
709{
710  if (n_IsZero(a, src)) return NULL;
711  assume(src == dst->extRing->cf);
712  poly result = p_One(dst->extRing);
713  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
714  return (number)result;
715}
716
717/* assumes that src = Z/u, dst = Z/p(a), where u != p */
718number naMapUP(number a, const coeffs src, const coeffs dst)
719{
720  if (n_IsZero(a, src)) return NULL;
721  /* mapping via intermediate int: */
722  int n = n_Int(a, src);
723  number q = n_Init(n, dst->extRing->cf);
724  poly result = p_One(dst->extRing);
725  p_SetCoeff(result, q, dst->extRing);
726  return (number)result;
727}
728
729nMapFunc naSetMap(const coeffs src, const coeffs dst)
730{
731  /* dst is expected to be an algebraic field extension */
732  assume(getCoeffType(dst) == ID);
733
734  if( src == dst ) return ndCopyMap;
735
736  int h = 0; /* the height of the extension tower given by dst */
737  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
738  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
739
740  /* for the time being, we only provide maps if h = 1 and if b is Q or
741     some field Z/pZ: */
742  if (h==0)
743  {
744    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
745      return naMap00;                            /// Q     -->  Q(a)
746    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
747      return naMapP0;                            /// Z/p   -->  Q(a)
748    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
749      return naMap0P;                            /// Q      --> Z/p(a)
750    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
751    {
752      if (src->ch == dst->ch) return naMapPP;    /// Z/p    --> Z/p(a)
753      else return naMapUP;                       /// Z/u    --> Z/p(a)
754    }
755  }
756  if (h != 1) return NULL;
757  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
758  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
759
760  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
761  {
762    if (strcmp(rRingVar(0, src->extRing),
763               rRingVar(0, dst->extRing)) == 0)
764    {
765      if (src->type==n_algExt)
766         return ndCopyMap; // naCopyMap;         /// Q(a)   --> Q(a)
767      else
768         return naCopyExt;
769    }
770    else
771      return NULL;                               /// Q(b)   --> Q(a)
772  }
773
774  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
775  {
776    if (strcmp(rRingVar(0,src->extRing),rRingVar(0,dst->extRing))==0)
777    {
778      if (src->type==n_algExt)
779        return ndCopyMap; // naCopyMap;          /// Z/p(a) --> Z/p(a)
780      else
781         return naCopyExt;
782    }
783    else
784      return NULL;                               /// Z/p(b) --> Z/p(a)
785  }
786
787  return NULL;                                           /// default
788}
789
790static int naParDeg(number a, const coeffs cf)
791{
792  if (a == NULL) return -1;
793  poly aa=(poly)a;
794  return cf->extRing->pFDeg(aa,cf->extRing);
795}
796
797/// return the specified parameter as a number in the given alg. field
798static number naParameter(const int iParameter, const coeffs cf)
799{
800  assume(getCoeffType(cf) == ID);
801
802  const ring R = cf->extRing;
803  assume( R != NULL );
804  assume( 0 < iParameter && iParameter <= rVar(R) );
805
806  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
807
808  return (number) p;
809}
810
811
812/// if m == var(i)/1 => return i,
813int naIsParam(number m, const coeffs cf)
814{
815  assume(getCoeffType(cf) == ID);
816
817  const ring R = cf->extRing;
818  assume( R != NULL );
819
820  return p_Var( (poly)m, R );
821}
822
823BOOLEAN naInitChar(coeffs cf, void * infoStruct)
824{
825  assume( infoStruct != NULL );
826
827  AlgExtInfo *e = (AlgExtInfo *)infoStruct;
828  /// first check whether cf->extRing != NULL and delete old ring???
829
830  assume(e->r                     != NULL);      // extRing;
831  assume(e->r->cf                 != NULL);      // extRing->cf;
832
833  assume((e->r->qideal            != NULL) &&    // minideal has one
834         (IDELEMS(e->r->qideal)   == 1)    &&    // non-zero generator
835         (e->r->qideal->m[0]      != NULL)    ); // at m[0];
836
837  assume( cf != NULL );
838  assume(getCoeffType(cf) == ID);                     // coeff type;
839
840  e->r->ref ++; // increase the ref.counter for the ground poly. ring!
841  const ring R = e->r; // no copy!
842  assume( R->qideal == e->r->qideal );
843  cf->extRing  = R;
844
845  /* propagate characteristic up so that it becomes
846     directly accessible in cf: */
847  cf->ch = R->cf->ch;
848
849  #ifdef LDEBUG
850  p_Test((poly)naMinpoly, naRing);
851  #endif
852
853  cf->cfGreaterZero  = naGreaterZero;
854  cf->cfGreater      = naGreater;
855  cf->cfEqual        = naEqual;
856  cf->cfIsZero       = naIsZero;
857  cf->cfIsOne        = naIsOne;
858  cf->cfIsMOne       = naIsMOne;
859  cf->cfInit         = naInit;
860  cf->cfInit_bigint  = naInit_bigint;
861  cf->cfInt          = naInt;
862  cf->cfNeg          = naNeg;
864  cf->cfSub          = naSub;
865  cf->cfMult         = naMult;
868  cf->cfPower        = naPower;
869  cf->cfCopy         = naCopy;
870
871  cf->cfWriteLong        = naWriteLong;
872
873  if( rCanShortOut(naRing) )
874    cf->cfWriteShort = naWriteShort;
875  else
876    cf->cfWriteShort = naWriteLong;
877
880  cf->cfSetMap       = naSetMap;
881  cf->cfGetDenom     = naGetDenom;
882  cf->cfGetNumerator = naGetNumerator;
883  cf->cfRePart       = naCopy;
884  cf->cfImPart       = naImPart;
885  cf->cfCoeffWrite   = naCoeffWrite;
886  cf->cfNormalize    = naNormalize;
887#ifdef LDEBUG
889#endif
890  cf->cfGcd          = naGcd;
891  //cf->cfLcm          = naLcm;
892  cf->cfSize         = naSize;
893  cf->nCoeffIsEqual  = naCoeffIsEqual;
894  cf->cfInvers       = naInvers;