source: git/libpolys/polys/ext_fields/algext.cc @ eea50a

spielwiese
Last change on this file since eea50a was bca341, checked in by Hans Schoenemann <hannes@…>, 12 years ago
fix: naIsOne, naIsMOne improved: algext.cc::definiteReduce
  • Property mode set to 100644
File size: 25.4 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
56#define naTest(a) naDBTest(a,__FILE__,__LINE__,cf)
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
70 * coefficients of our numbers, e.g., use naCoeffs->nAdd to add
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
172  int aDeg = 0;
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);
181  naDelete(&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 ((a==NULL) || (!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 ((a==NULL) || (!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;
283  int aDeg = 0;
284  if (a != NULL) aDeg = p_Totaldegree((poly)a, naRing);
285  int bDeg = 0;
286  if (b != NULL) bDeg = p_Totaldegree((poly)b, naRing);
287  return (aDeg > bDeg);
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);
390  definiteReduce(aDivB, naMinpoly, cf);
391  return (number)aDivB;
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);
457    naDelete(&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);
535  return naDiv(theProduct, theGcd, 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  if ((p!=NULL) && (p_GetExp(p,1,naRing)>=p_GetExp(reducer,1,naRing)))
603  {
604    p_PolyDiv(p, reducer, FALSE, naRing);
605  }
606}
607
608void  naNormalize(number &a, const coeffs cf)
609{
610  poly aa=(poly)a;
611  if (aa!=naMinpoly)
612    definiteReduce(aa,naMinpoly,cf);
613  a=(number)aa;
614}
615
616#ifdef HAVE_FACTORY
617number naConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
618{
619  if (n.isZero()) return NULL;
620  poly p=convFactoryPSingP(n,naRing);
621  return (number)p;
622}
623CanonicalForm naConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
624{
625  naTest(n);
626  if (n==NULL) return CanonicalForm(0);
627
628  return convSingPFactoryP((poly)n,naRing);
629}
630#endif
631
632
633/* IMPORTANT NOTE: Since an algebraic field extension is again a field,
634                   the gcd of two elements is not very interesting. (It
635                   is actually any unit in the field, i.e., any non-
636                   zero element.) Note that the below method does not operate
637                   in this strong sense but rather computes the gcd of
638                   two given elements in the underlying polynomial ring. */
639number naGcd(number a, number b, const coeffs cf)
640{
641  naTest(a); naTest(b);
642  if ((a == NULL) && (b == NULL)) WerrorS(nDivBy0);
643  return (number)p_Gcd((poly)a, (poly)b, naRing);
644}
645
646number naInvers(number a, const coeffs cf)
647{
648  naTest(a);
649  if (a == NULL) WerrorS(nDivBy0);
650  poly aFactor = NULL; poly mFactor = NULL;
651  poly theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
652  naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor);
653  /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */
654  assume(naIsOne((number)theGcd, cf));
655  p_Delete(&theGcd, naRing);
656  p_Delete(&mFactor, naRing);
657  return (number)(aFactor);
658}
659
660/* assumes that src = Q, dst = Q(a) */
661number naMap00(number a, const coeffs src, const coeffs dst)
662{
663  if (n_IsZero(a, src)) return NULL;
664  assume(src == dst->extRing->cf);
665  poly result = p_One(dst->extRing);
666  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
667  return (number)result;
668}
669
670/* assumes that src = Z/p, dst = Q(a) */
671number naMapP0(number a, const coeffs src, const coeffs dst)
672{
673  if (n_IsZero(a, src)) return NULL;
674  /* mapping via intermediate int: */
675  int n = n_Int(a, src);
676  number q = n_Init(n, dst->extRing->cf);
677  poly result = p_One(dst->extRing);
678  p_SetCoeff(result, q, dst->extRing);
679  return (number)result;
680}
681
682#if 0
683/* assumes that either src = Q(a), dst = Q(a), or
684                       src = Z/p(a), dst = Z/p(a)     */
685number naCopyMap(number a, const coeffs src, const coeffs dst)
686{
687  return naCopy(a, dst);
688}
689#endif
690
691number naCopyExt(number a, const coeffs src, const coeffs)
692{
693  fraction fa=(fraction)a;
694  return (number)p_Copy(NUM(fa),src->extRing);
695}
696
697/* assumes that src = Q, dst = Z/p(a) */
698number naMap0P(number a, const coeffs src, const coeffs dst)
699{
700  if (n_IsZero(a, src)) return NULL;
701  int p = rChar(dst->extRing);
702
703  number q = nlModP(a, src, dst->extRing->cf);
704
705  poly result = p_NSet(q, dst->extRing);
706
707  return (number)result;
708}
709
710/* assumes that src = Z/p, dst = Z/p(a) */
711number naMapPP(number a, const coeffs src, const coeffs dst)
712{
713  if (n_IsZero(a, src)) return NULL;
714  assume(src == dst->extRing->cf);
715  poly result = p_One(dst->extRing);
716  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
717  return (number)result;
718}
719
720/* assumes that src = Z/u, dst = Z/p(a), where u != p */
721number naMapUP(number a, const coeffs src, const coeffs dst)
722{
723  if (n_IsZero(a, src)) return NULL;
724  /* mapping via intermediate int: */
725  int n = n_Int(a, src);
726  number q = n_Init(n, dst->extRing->cf);
727  poly result = p_One(dst->extRing);
728  p_SetCoeff(result, q, dst->extRing);
729  return (number)result;
730}
731
732nMapFunc naSetMap(const coeffs src, const coeffs dst)
733{
734  /* dst is expected to be an algebraic field extension */
735  assume(getCoeffType(dst) == ID);
736
737  if( src == dst ) return ndCopyMap;
738
739  int h = 0; /* the height of the extension tower given by dst */
740  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
741  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
742
743  /* for the time being, we only provide maps if h = 1 and if b is Q or
744     some field Z/pZ: */
745  if (h==0)
746  {
747    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
748      return naMap00;                            /// Q     -->  Q(a)
749    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
750      return naMapP0;                            /// Z/p   -->  Q(a)
751    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
752      return naMap0P;                            /// Q      --> Z/p(a)
753    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
754    {
755      if (src->ch == dst->ch) return naMapPP;    /// Z/p    --> Z/p(a)
756      else return naMapUP;                       /// Z/u    --> Z/p(a)
757    }
758  }
759  if (h != 1) return NULL;
760  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
761  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
762
763  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
764  {
765    if (strcmp(rRingVar(0, src->extRing),
766               rRingVar(0, dst->extRing)) == 0)
767    {
768      if (src->type==n_algExt)
769         return ndCopyMap; // naCopyMap;         /// Q(a)   --> Q(a)
770      else
771         return naCopyExt;
772    }
773    else
774      return NULL;                               /// Q(b)   --> Q(a)
775  }
776
777  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
778  {
779    if (strcmp(rRingVar(0,src->extRing),rRingVar(0,dst->extRing))==0)
780    {
781      if (src->type==n_algExt)
782        return ndCopyMap; // naCopyMap;          /// Z/p(a) --> Z/p(a)
783      else
784         return naCopyExt;
785    }
786    else
787      return NULL;                               /// Z/p(b) --> Z/p(a)
788  }
789
790  return NULL;                                           /// default
791}
792
793static int naParDeg(number a, const coeffs cf)
794{
795  if (a == NULL) return -1;
796  poly aa=(poly)a;
797  return cf->extRing->pFDeg(aa,cf->extRing);
798}
799
800/// return the specified parameter as a number in the given alg. field
801static number naParameter(const int iParameter, const coeffs cf)
802{
803  assume(getCoeffType(cf) == ID);
804
805  const ring R = cf->extRing;
806  assume( R != NULL );
807  assume( 0 < iParameter && iParameter <= rVar(R) );
808
809  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
810
811  return (number) p;
812}
813
814
815/// if m == var(i)/1 => return i,
816int naIsParam(number m, const coeffs cf)
817{
818  assume(getCoeffType(cf) == ID);
819
820  const ring R = cf->extRing;
821  assume( R != NULL );
822
823  return p_Var( (poly)m, R );
824}
825
826BOOLEAN naInitChar(coeffs cf, void * infoStruct)
827{
828  assume( infoStruct != NULL );
829
830  AlgExtInfo *e = (AlgExtInfo *)infoStruct;
831  /// first check whether cf->extRing != NULL and delete old ring???
832
833  assume(e->r                     != NULL);      // extRing;
834  assume(e->r->cf                 != NULL);      // extRing->cf;
835
836  assume((e->r->qideal            != NULL) &&    // minideal has one
837         (IDELEMS(e->r->qideal)   == 1)    &&    // non-zero generator
838         (e->r->qideal->m[0]      != NULL)    ); // at m[0];
839
840  assume( cf != NULL );
841  assume(getCoeffType(cf) == ID);                     // coeff type;
842
843  e->r->ref ++; // increase the ref.counter for the ground poly. ring!
844  const ring R = e->r; // no copy!
845  assume( R->qideal == e->r->qideal );
846  cf->extRing  = R;
847
848  /* propagate characteristic up so that it becomes
849     directly accessible in cf: */
850  cf->ch = R->cf->ch;
851
852  #ifdef LDEBUG
853  p_Test((poly)naMinpoly, naRing);
854  #endif
855
856  cf->cfGreaterZero  = naGreaterZero;
857  cf->cfGreater      = naGreater;
858  cf->cfEqual        = naEqual;
859  cf->cfIsZero       = naIsZero;
860  cf->cfIsOne        = naIsOne;
861  cf->cfIsMOne       = naIsMOne;
862  cf->cfInit         = naInit;
863  cf->cfInit_bigint  = naInit_bigint;
864  cf->cfInt          = naInt;
865  cf->cfNeg          = naNeg;
866  cf->cfAdd          = naAdd;
867  cf->cfSub          = naSub;
868  cf->cfMult         = naMult;
869  cf->cfDiv          = naDiv;
870  cf->cfExactDiv     = naDiv;
871  cf->cfPower        = naPower;
872  cf->cfCopy         = naCopy;
873
874  cf->cfWriteLong        = naWriteLong;
875
876  if( rCanShortOut(naRing) )
877    cf->cfWriteShort = naWriteShort;
878  else
879    cf->cfWriteShort = naWriteLong;
880
881  cf->cfRead         = naRead;
882  cf->cfDelete       = naDelete;
883  cf->cfSetMap       = naSetMap;
884  cf->cfGetDenom     = naGetDenom;
885  cf->cfGetNumerator = naGetNumerator;
886  cf->cfRePart       = naCopy;
887  cf->cfImPart       = naImPart;
888  cf->cfCoeffWrite   = naCoeffWrite;
889  cf->cfNormalize    = naNormalize;
890#ifdef LDEBUG
891  cf->cfDBTest       = naDBTest;
892#endif
893  cf->cfGcd          = naGcd;
894  //cf->cfLcm          = naLcm;
895  cf->cfSize         = naSize;
896  cf->nCoeffIsEqual  = naCoeffIsEqual;
897  cf->cfInvers       = naInvers;
898  cf->cfIntDiv       = naDiv;
899#ifdef HAVE_FACTORY
900  cf->convFactoryNSingN=naConvFactoryNSingN;
901  cf->convSingNFactoryN=naConvSingNFactoryN;
902#endif
903  cf->cfParDeg = naParDeg;
904
905  cf->iNumberOfParameters = rVar(R);
906  cf->pParameterNames = R->names;
907  cf->cfParameter = naParameter;
908
909  return FALSE;
910}
Note: See TracBrowser for help on using the repository browser.