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

spielwiese
Last change on this file since afda22 was afda22, checked in by Hans Schoenemann <hannes@…>, 12 years ago
fix: #329, part 1: alg. ext
  • Property mode set to 100644
File size: 25.6 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) or (deg(a)==deg(b) && lc(a)>lc(b))) */
278BOOLEAN naGreater(number a, number b, const coeffs cf)
279{
280  naTest(a); naTest(b);
281  if (naIsZero(a, cf))
282  {
283    if (naIsZero(b, cf)) return FALSE;
284    return !n_GreaterZero(pGetCoeff((poly)b),cf);
285  }
286  if (naIsZero(b, cf))
287  {
288    return n_GreaterZero(pGetCoeff((poly)a),cf);
289  }
290  int aDeg = p_Totaldegree((poly)a, naRing);
291  int bDeg = p_Totaldegree((poly)b, naRing);
292  if (aDeg>bDeg) return TRUE;
293  if (aDeg<bDeg) return FALSE;
294  return n_Greater(pGetCoeff((poly)a),pGetCoeff((poly)b),cf);
295}
296
297/* TRUE iff a != 0 and (LC(a) > 0 or deg(a) > 0) */
298BOOLEAN naGreaterZero(number a, const coeffs cf)
299{
300  naTest(a);
301  if (a == NULL)                                            return FALSE;
302  if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE;
303  if (p_Totaldegree((poly)a, naRing) > 0)                   return TRUE;
304  return FALSE;
305}
306
307void naCoeffWrite(const coeffs cf, BOOLEAN details)
308{
309  assume( cf != NULL );
310
311  const ring A = cf->extRing;
312
313  assume( A != NULL );
314  assume( A->cf != NULL );
315
316  n_CoeffWrite(A->cf, details);
317
318//  rWrite(A);
319
320  const int P = rVar(A);
321  assume( P > 0 );
322
323  Print("//   %d parameter    : ", P);
324
325  for (int nop=0; nop < P; nop ++)
326    Print("%s ", rRingVar(nop, A));
327
328  PrintLn();
329
330  const ideal I = A->qideal;
331
332  assume( I != NULL );
333  assume( IDELEMS(I) == 1 );
334
335
336  if ( details )
337  {
338    PrintS("//   minpoly        : (");
339    p_Write0( I->m[0], A);
340    PrintS(")");
341  }
342  else
343    PrintS("//   minpoly        : ...");
344
345  PrintLn();
346
347/*
348  char *x = rRingVar(0, A);
349
350  Print("//   Coefficients live in the extension field K[%s]/<f(%s)>\n", x, x);
351  Print("//   with the minimal polynomial f(%s) = %s\n", x,
352        p_String(A->qideal->m[0], A));
353  PrintS("//   and K: ");
354*/
355}
356
357number naAdd(number a, number b, const coeffs cf)
358{
359  naTest(a); naTest(b);
360  if (a == NULL) return naCopy(b, cf);
361  if (b == NULL) return naCopy(a, cf);
362  poly aPlusB = p_Add_q(p_Copy((poly)a, naRing),
363                        p_Copy((poly)b, naRing), naRing);
364  definiteReduce(aPlusB, naMinpoly, cf);
365  return (number)aPlusB;
366}
367
368number naSub(number a, number b, const coeffs cf)
369{
370  naTest(a); naTest(b);
371  if (b == NULL) return naCopy(a, cf);
372  poly minusB = p_Neg(p_Copy((poly)b, naRing), naRing);
373  if (a == NULL) return (number)minusB;
374  poly aMinusB = p_Add_q(p_Copy((poly)a, naRing), minusB, naRing);
375  definiteReduce(aMinusB, naMinpoly, cf);
376  return (number)aMinusB;
377}
378
379number naMult(number a, number b, const coeffs cf)
380{
381  naTest(a); naTest(b);
382  if (a == NULL) return NULL;
383  if (b == NULL) return NULL;
384  poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing),
385                          p_Copy((poly)b, naRing), naRing);
386  definiteReduce(aTimesB, naMinpoly, cf);
387  return (number)aTimesB;
388}
389
390number naDiv(number a, number b, const coeffs cf)
391{
392  naTest(a); naTest(b);
393  if (b == NULL) WerrorS(nDivBy0);
394  if (a == NULL) return NULL;
395  poly bInverse = (poly)naInvers(b, cf);
396  poly aDivB = p_Mult_q(p_Copy((poly)a, naRing), bInverse, naRing);
397  definiteReduce(aDivB, naMinpoly, cf);
398  return (number)aDivB;
399}
400
401/* 0^0 = 0;
402   for |exp| <= 7 compute power by a simple multiplication loop;
403   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
404      p^13 = p^1 * p^4 * p^8, where we utilise that
405      p^(2^(k+1)) = p^(2^k) * p^(2^k);
406   intermediate reduction modulo the minimal polynomial is controlled by
407   the in-place method heuristicReduce(poly, poly, coeffs); see there.
408*/
409void naPower(number a, int exp, number *b, const coeffs cf)
410{
411  naTest(a);
412
413  /* special cases first */
414  if (a == NULL)
415  {
416    if (exp >= 0) *b = NULL;
417    else          WerrorS(nDivBy0);
418  }
419  else if (exp ==  0) { *b = naInit(1, cf); return; }
420  else if (exp ==  1) { *b = naCopy(a, cf); return; }
421  else if (exp == -1) { *b = naInvers(a, cf); return; }
422
423  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
424
425  /* now compute a^expAbs */
426  poly pow; poly aAsPoly = (poly)a;
427  if (expAbs <= 7)
428  {
429    pow = p_Copy(aAsPoly, naRing);
430    for (int i = 2; i <= expAbs; i++)
431    {
432      pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing);
433      heuristicReduce(pow, naMinpoly, cf);
434    }
435    definiteReduce(pow, naMinpoly, cf);
436  }
437  else
438  {
439    pow = p_ISet(1, naRing);
440    poly factor = p_Copy(aAsPoly, naRing);
441    while (expAbs != 0)
442    {
443      if (expAbs & 1)
444      {
445        pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing);
446        heuristicReduce(pow, naMinpoly, cf);
447      }
448      expAbs = expAbs / 2;
449      if (expAbs != 0)
450      {
451        factor = p_Mult_q(factor, p_Copy(factor, naRing), naRing);
452        heuristicReduce(factor, naMinpoly, cf);
453      }
454    }
455    p_Delete(&factor, naRing);
456    definiteReduce(pow, naMinpoly, cf);
457  }
458
459  /* invert if original exponent was negative */
460  number n = (number)pow;
461  if (exp < 0)
462  {
463    number m = naInvers(n, cf);
464    naDelete(&n, cf);
465    n = m;
466  }
467  *b = n;
468}
469
470/* may reduce p modulo the reducer by calling definiteReduce;
471   the decision is made based on the following heuristic
472   (which should also only be changed here in this method):
473      if (deg(p) > 10*deg(reducer) then perform reduction;
474   modifies p */
475void heuristicReduce(poly &p, poly reducer, const coeffs cf)
476{
477  #ifdef LDEBUG
478  p_Test((poly)p, naRing);
479  p_Test((poly)reducer, naRing);
480  #endif
481  if (p_Totaldegree(p, naRing) > 10 * p_Totaldegree(reducer, naRing))
482    definiteReduce(p, reducer, cf);
483}
484
485void naWriteLong(number &a, const coeffs cf)
486{
487  naTest(a);
488  if (a == NULL)
489    StringAppendS("0");
490  else
491  {
492    poly aAsPoly = (poly)a;
493    /* basically, just write aAsPoly using p_Write,
494       but use brackets around the output, if a is not
495       a constant living in naCoeffs = cf->extRing->cf */
496    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
497    if (useBrackets) StringAppendS("(");
498    p_String0Long(aAsPoly, naRing, naRing);
499    if (useBrackets) StringAppendS(")");
500  }
501}
502
503void naWriteShort(number &a, const coeffs cf)
504{
505  naTest(a);
506  if (a == NULL)
507    StringAppendS("0");
508  else
509  {
510    poly aAsPoly = (poly)a;
511    /* basically, just write aAsPoly using p_Write,
512       but use brackets around the output, if a is not
513       a constant living in naCoeffs = cf->extRing->cf */
514    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
515    if (useBrackets) StringAppendS("(");
516    p_String0Short(aAsPoly, naRing, naRing);
517    if (useBrackets) StringAppendS(")");
518  }
519}
520
521const char * naRead(const char *s, number *a, const coeffs cf)
522{
523  poly aAsPoly;
524  const char * result = p_Read(s, aAsPoly, naRing);
525  definiteReduce(aAsPoly, naMinpoly, cf);
526  *a = (number)aAsPoly;
527  return result;
528}
529
530#if 0
531/* implemented by the rule lcm(a, b) = a * b / gcd(a, b) */
532number naLcm(number a, number b, const coeffs cf)
533{
534  naTest(a); naTest(b);
535  if (a == NULL) return NULL;
536  if (b == NULL) return NULL;
537  number theProduct = (number)p_Mult_q(p_Copy((poly)a, naRing),
538                                       p_Copy((poly)b, naRing), naRing);
539  /* note that theProduct needs not be reduced w.r.t. naMinpoly;
540     but the final division will take care of the necessary reduction */
541  number theGcd = naGcd(a, b, cf);
542  return naDiv(theProduct, theGcd, cf);
543}
544#endif
545
546/* expects *param to be castable to AlgExtInfo */
547static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
548{
549  if (ID != n) return FALSE;
550  AlgExtInfo *e = (AlgExtInfo *)param;
551  /* for extension coefficient fields we expect the underlying
552     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
553     this expectation is based on the assumption that we have properly
554     registered cf and perform reference counting rather than creating
555     multiple copies of the same coefficient field/domain/ring */
556  if (naRing == e->r)
557    return TRUE;
558  /* (Note that then also the minimal ideals will necessarily be
559     the same, as they are attached to the ring.) */
560
561  // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)...
562  if( rEqual(naRing, e->r, TRUE) ) // also checks the equality of qideals
563  {
564    const ideal mi = naRing->qideal; 
565    assume( IDELEMS(mi) == 1 );
566    const ideal ii = e->r->qideal;
567    assume( IDELEMS(ii) == 1 );
568
569    // TODO: the following should be extended for 2 *equal* rings...
570    assume( p_EqualPolys(mi->m[0], ii->m[0], naRing, e->r) );
571   
572    rDelete(e->r);
573   
574    return TRUE;
575  }
576
577  return FALSE;
578
579}
580
581int naSize(number a, const coeffs cf)
582{
583  if (a == NULL) return -1;
584  /* this has been taken from the old implementation of field extensions,
585     where we computed the sum of the degree and the number of terms in
586     (poly)a; so we leave it at that, for the time being;
587     maybe, the number of terms alone is a better measure? */
588  poly aAsPoly = (poly)a;
589  int theDegree = 0; int noOfTerms = 0;
590  while (aAsPoly != NULL)
591  {
592    noOfTerms++;
593    int d = p_GetExp(aAsPoly, 1, naRing);
594    if (d > theDegree) theDegree = d;
595    pIter(aAsPoly);
596  }
597  return theDegree + noOfTerms;
598}
599
600/* performs polynomial division and overrides p by the remainder
601   of division of p by the reducer;
602   modifies p */
603void definiteReduce(poly &p, poly reducer, const coeffs cf)
604{
605  #ifdef LDEBUG
606  p_Test((poly)p, naRing);
607  p_Test((poly)reducer, naRing);
608  #endif
609  if ((p!=NULL) && (p_GetExp(p,1,naRing)>=p_GetExp(reducer,1,naRing)))
610  {
611    p_PolyDiv(p, reducer, FALSE, naRing);
612  }
613}
614
615void  naNormalize(number &a, const coeffs cf)
616{
617  poly aa=(poly)a;
618  if (aa!=naMinpoly)
619    definiteReduce(aa,naMinpoly,cf);
620  a=(number)aa;
621}
622
623#ifdef HAVE_FACTORY
624number naConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
625{
626  if (n.isZero()) return NULL;
627  poly p=convFactoryPSingP(n,naRing);
628  return (number)p;
629}
630CanonicalForm naConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
631{
632  naTest(n);
633  if (n==NULL) return CanonicalForm(0);
634
635  return convSingPFactoryP((poly)n,naRing);
636}
637#endif
638
639
640/* IMPORTANT NOTE: Since an algebraic field extension is again a field,
641                   the gcd of two elements is not very interesting. (It
642                   is actually any unit in the field, i.e., any non-
643                   zero element.) Note that the below method does not operate
644                   in this strong sense but rather computes the gcd of
645                   two given elements in the underlying polynomial ring. */
646number naGcd(number a, number b, const coeffs cf)
647{
648  naTest(a); naTest(b);
649  if ((a == NULL) && (b == NULL)) WerrorS(nDivBy0);
650  return (number)p_Gcd((poly)a, (poly)b, naRing);
651}
652
653number naInvers(number a, const coeffs cf)
654{
655  naTest(a);
656  if (a == NULL) WerrorS(nDivBy0);
657  poly aFactor = NULL; poly mFactor = NULL;
658  poly theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
659  naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor);
660  /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */
661  assume(naIsOne((number)theGcd, cf));
662  p_Delete(&theGcd, naRing);
663  p_Delete(&mFactor, naRing);
664  return (number)(aFactor);
665}
666
667/* assumes that src = Q, dst = Q(a) */
668number naMap00(number a, const coeffs src, const coeffs dst)
669{
670  if (n_IsZero(a, src)) return NULL;
671  assume(src == dst->extRing->cf);
672  poly result = p_One(dst->extRing);
673  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
674  return (number)result;
675}
676
677/* assumes that src = Z/p, dst = Q(a) */
678number naMapP0(number a, const coeffs src, const coeffs dst)
679{
680  if (n_IsZero(a, src)) return NULL;
681  /* mapping via intermediate int: */
682  int n = n_Int(a, src);
683  number q = n_Init(n, dst->extRing->cf);
684  poly result = p_One(dst->extRing);
685  p_SetCoeff(result, q, dst->extRing);
686  return (number)result;
687}
688
689#if 0
690/* assumes that either src = Q(a), dst = Q(a), or
691                       src = Z/p(a), dst = Z/p(a)     */
692number naCopyMap(number a, const coeffs src, const coeffs dst)
693{
694  return naCopy(a, dst);
695}
696#endif
697
698number naCopyExt(number a, const coeffs src, const coeffs)
699{
700  fraction fa=(fraction)a;
701  return (number)p_Copy(NUM(fa),src->extRing);
702}
703
704/* assumes that src = Q, dst = Z/p(a) */
705number naMap0P(number a, const coeffs src, const coeffs dst)
706{
707  if (n_IsZero(a, src)) return NULL;
708  int p = rChar(dst->extRing);
709
710  number q = nlModP(a, src, dst->extRing->cf);
711
712  poly result = p_NSet(q, dst->extRing);
713
714  return (number)result;
715}
716
717/* assumes that src = Z/p, dst = Z/p(a) */
718number naMapPP(number a, const coeffs src, const coeffs dst)
719{
720  if (n_IsZero(a, src)) return NULL;
721  assume(src == dst->extRing->cf);
722  poly result = p_One(dst->extRing);
723  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
724  return (number)result;
725}
726
727/* assumes that src = Z/u, dst = Z/p(a), where u != p */
728number naMapUP(number a, const coeffs src, const coeffs dst)
729{
730  if (n_IsZero(a, src)) return NULL;
731  /* mapping via intermediate int: */
732  int n = n_Int(a, src);
733  number q = n_Init(n, dst->extRing->cf);
734  poly result = p_One(dst->extRing);
735  p_SetCoeff(result, q, dst->extRing);
736  return (number)result;
737}
738
739nMapFunc naSetMap(const coeffs src, const coeffs dst)
740{
741  /* dst is expected to be an algebraic field extension */
742  assume(getCoeffType(dst) == ID);
743
744  if( src == dst ) return ndCopyMap;
745
746  int h = 0; /* the height of the extension tower given by dst */
747  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
748  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
749
750  /* for the time being, we only provide maps if h = 1 and if b is Q or
751     some field Z/pZ: */
752  if (h==0)
753  {
754    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
755      return naMap00;                            /// Q     -->  Q(a)
756    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
757      return naMapP0;                            /// Z/p   -->  Q(a)
758    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
759      return naMap0P;                            /// Q      --> Z/p(a)
760    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
761    {
762      if (src->ch == dst->ch) return naMapPP;    /// Z/p    --> Z/p(a)
763      else return naMapUP;                       /// Z/u    --> Z/p(a)
764    }
765  }
766  if (h != 1) return NULL;
767  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
768  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
769
770  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
771  {
772    if (strcmp(rRingVar(0, src->extRing),
773               rRingVar(0, dst->extRing)) == 0)
774    {
775      if (src->type==n_algExt)
776         return ndCopyMap; // naCopyMap;         /// Q(a)   --> Q(a)
777      else
778         return naCopyExt;
779    }
780    else
781      return NULL;                               /// Q(b)   --> Q(a)
782  }
783
784  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
785  {
786    if (strcmp(rRingVar(0,src->extRing),rRingVar(0,dst->extRing))==0)
787    {
788      if (src->type==n_algExt)
789        return ndCopyMap; // naCopyMap;          /// Z/p(a) --> Z/p(a)
790      else
791         return naCopyExt;
792    }
793    else
794      return NULL;                               /// Z/p(b) --> Z/p(a)
795  }
796
797  return NULL;                                           /// default
798}
799
800static int naParDeg(number a, const coeffs cf)
801{
802  if (a == NULL) return -1;
803  poly aa=(poly)a;
804  return cf->extRing->pFDeg(aa,cf->extRing);
805}
806
807/// return the specified parameter as a number in the given alg. field
808static number naParameter(const int iParameter, const coeffs cf)
809{
810  assume(getCoeffType(cf) == ID);
811
812  const ring R = cf->extRing;
813  assume( R != NULL );
814  assume( 0 < iParameter && iParameter <= rVar(R) );
815
816  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
817
818  return (number) p;
819}
820
821
822/// if m == var(i)/1 => return i,
823int naIsParam(number m, const coeffs cf)
824{
825  assume(getCoeffType(cf) == ID);
826
827  const ring R = cf->extRing;
828  assume( R != NULL );
829
830  return p_Var( (poly)m, R );
831}
832
833BOOLEAN naInitChar(coeffs cf, void * infoStruct)
834{
835  assume( infoStruct != NULL );
836
837  AlgExtInfo *e = (AlgExtInfo *)infoStruct;
838  /// first check whether cf->extRing != NULL and delete old ring???
839
840  assume(e->r                     != NULL);      // extRing;
841  assume(e->r->cf                 != NULL);      // extRing->cf;
842
843  assume((e->r->qideal            != NULL) &&    // minideal has one
844         (IDELEMS(e->r->qideal)   == 1)    &&    // non-zero generator
845         (e->r->qideal->m[0]      != NULL)    ); // at m[0];
846
847  assume( cf != NULL );
848  assume(getCoeffType(cf) == ID);                     // coeff type;
849
850  e->r->ref ++; // increase the ref.counter for the ground poly. ring!
851  const ring R = e->r; // no copy!
852  assume( R->qideal == e->r->qideal );
853  cf->extRing  = R;
854
855  /* propagate characteristic up so that it becomes
856     directly accessible in cf: */
857  cf->ch = R->cf->ch;
858
859  #ifdef LDEBUG
860  p_Test((poly)naMinpoly, naRing);
861  #endif
862
863  cf->cfGreaterZero  = naGreaterZero;
864  cf->cfGreater      = naGreater;
865  cf->cfEqual        = naEqual;
866  cf->cfIsZero       = naIsZero;
867  cf->cfIsOne        = naIsOne;
868  cf->cfIsMOne       = naIsMOne;
869  cf->cfInit         = naInit;
870  cf->cfInit_bigint  = naInit_bigint;
871  cf->cfInt          = naInt;
872  cf->cfNeg          = naNeg;
873  cf->cfAdd          = naAdd;
874  cf->cfSub          = naSub;
875  cf->cfMult         = naMult;
876  cf->cfDiv          = naDiv;
877  cf->cfExactDiv     = naDiv;
878  cf->cfPower        = naPower;
879  cf->cfCopy         = naCopy;
880
881  cf->cfWriteLong        = naWriteLong;
882
883  if( rCanShortOut(naRing) )
884    cf->cfWriteShort = naWriteShort;
885  else
886    cf->cfWriteShort = naWriteLong;
887
888  cf->cfRead         = naRead;
889  cf->cfDelete       = naDelete;
890  cf->cfSetMap       = naSetMap;
891  cf->cfGetDenom     = naGetDenom;
892  cf->cfGetNumerator = naGetNumerator;
893  cf->cfRePart       = naCopy;
894  cf->cfImPart       = naImPart;
895  cf->cfCoeffWrite   = naCoeffWrite;
896  cf->cfNormalize    = naNormalize;
897#ifdef LDEBUG
898  cf->cfDBTest       = naDBTest;
899#endif
900  cf->cfGcd          = naGcd;
901  //cf->cfLcm          = naLcm;
902  cf->cfSize         = naSize;
903  cf->nCoeffIsEqual  = naCoeffIsEqual;
904  cf->cfInvers       = naInvers;
905  cf->cfIntDiv       = naDiv;
906#ifdef HAVE_FACTORY
907  cf->convFactoryNSingN=naConvFactoryNSingN;
908  cf->convSingNFactoryN=naConvSingNFactoryN;
909#endif
910  cf->cfParDeg = naParDeg;
911
912  cf->iNumberOfParameters = rVar(R);
913  cf->pParameterNames = R->names;
914  cf->cfParameter = naParameter;
915
916  return FALSE;
917}
Note: See TracBrowser for help on using the repository browser.