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

spielwiese
Last change on this file since fea494 was fea494, checked in by Hans Schoenemann <hannes@…>, 9 years ago
format
  • Property mode set to 100644
File size: 45.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
32
33
34#include <misc/auxiliary.h>
35
36#include <omalloc/omalloc.h>
37
38#include <reporter/reporter.h>
39
40#include <coeffs/coeffs.h>
41#include <coeffs/numbers.h>
42#include <coeffs/longrat.h>
43
44#include <polys/monomials/ring.h>
45#include <polys/monomials/p_polys.h>
46#include <polys/simpleideals.h>
47
48#include <polys/PolyEnumerator.h>
49
50#include <factory/factory.h>
51#include <polys/clapconv.h>
52#include <polys/clapsing.h>
53#include <polys/prCopy.h>
54
55#include <polys/ext_fields/algext.h>
56#define TRANSEXT_PRIVATES 1
57#include <polys/ext_fields/transext.h>
58
59#ifdef LDEBUG
60#define naTest(a) naDBTest(a,__FILE__,__LINE__,cf)
61BOOLEAN  naDBTest(number a, const char *f, const int l, const coeffs r);
62#else
63#define naTest(a) do {} while (0)
64#endif
65
66/// Our own type!
67static const n_coeffType ID = n_algExt;
68
69/* polynomial ring in which our numbers live */
70#define naRing cf->extRing
71
72/* coeffs object in which the coefficients of our numbers live;
73 * methods attached to naCoeffs may be used to compute with the
74 * coefficients of our numbers, e.g., use naCoeffs->nAdd to add
75 * coefficients of our numbers */
76#define naCoeffs cf->extRing->cf
77
78/* minimal polynomial */
79#define naMinpoly naRing->qideal->m[0]
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   naGetDenom(number &a, const coeffs cf);
101number   naGetNumerator(number &a, const coeffs cf);
102number   naGcd(number a, number b, const coeffs cf);
103int      naSize(number a, const coeffs cf);
104void     naDelete(number *a, const coeffs cf);
105void     naCoeffWrite(const coeffs cf, BOOLEAN details);
106//number   naIntDiv(number a, number b, const coeffs cf);
107const char * naRead(const char *s, number *a, const coeffs cf);
108
109static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param);
110
111
112/// returns NULL if p == NULL, otherwise makes p monic by dividing
113///   by its leading coefficient (only done if this is not already 1);
114///   this assumes that we are over a ground field so that division
115///   is well-defined;
116///   modifies p
117// void      p_Monic(poly p, const ring r);
118
119///   assumes that p and q are univariate polynomials in r,
120///   mentioning the same variable;
121///   assumes a global monomial ordering in r;
122///   assumes that not both p and q are NULL;
123///   returns the gcd of p and q;
124///   leaves p and q unmodified
125// poly      p_Gcd(const poly p, const poly q, const ring r);
126
127/* returns NULL if p == NULL, otherwise makes p monic by dividing
128   by its leading coefficient (only done if this is not already 1);
129   this assumes that we are over a ground field so that division
130   is well-defined;
131   modifies p */
132static inline void p_Monic(poly p, const ring r)
133{
134  if (p == NULL) return;
135  number n = n_Init(1, r->cf);
136  if (p->next==NULL) { p_SetCoeff(p,n,r); return; }
137  poly pp = p;
138  number lc = p_GetCoeff(p, r);
139  if (n_IsOne(lc, r->cf)) return;
140  number lcInverse = n_Invers(lc, r->cf);
141  p_SetCoeff(p, n, r);   // destroys old leading coefficient!
142  pIter(p);
143  while (p != NULL)
144  {
145    number n = n_Mult(p_GetCoeff(p, r), lcInverse, r->cf);
146    n_Normalize(n,r->cf);
147    p_SetCoeff(p, n, r);   // destroys old leading coefficient!
148    pIter(p);
149  }
150  n_Delete(&lcInverse, r->cf);
151  p = pp;
152}
153
154/// see p_Gcd;
155///   additional assumption: deg(p) >= deg(q);
156///   must destroy p and q (unless one of them is returned)
157static inline poly p_GcdHelper(poly &p, poly &q, const ring r)
158{
159  while (q != NULL)
160  {
161    p_PolyDiv(p, q, FALSE, r);
162    // swap p and q:
163    poly& t = q;
164    q = p;
165    p = t;
166
167  }
168  return p;
169}
170
171/* assumes that p and q are univariate polynomials in r,
172   mentioning the same variable;
173   assumes a global monomial ordering in r;
174   assumes that not both p and q are NULL;
175   returns the gcd of p and q;
176   leaves p and q unmodified */
177static inline poly      p_Gcd(const poly p, const poly q, const ring r)
178{
179  assume((p != NULL) || (q != NULL));
180
181  poly a = p; poly b = q;
182  if (p_Deg(a, r) < p_Deg(b, r)) { a = q; b = p; }
183  a = p_Copy(a, r); b = p_Copy(b, r);
184
185  /* We have to make p monic before we return it, so that if the
186     gcd is a unit in the ground field, we will actually return 1. */
187  a = p_GcdHelper(a, b, r);
188  p_Monic(a, r);
189  return a;
190}
191
192/* see p_ExtGcd;
193   additional assumption: deg(p) >= deg(q);
194   must destroy p and q (unless one of them is returned) */
195static inline poly p_ExtGcdHelper(poly &p, poly &pFactor, poly &q, poly &qFactor,
196                                  ring r)
197{
198  if (q == NULL)
199  {
200    qFactor = NULL;
201    pFactor = p_ISet(1, r);
202    p_SetCoeff(pFactor, n_Invers(p_GetCoeff(p, r), r->cf), r);
203    p_Monic(p, r);
204    return p;
205  }
206  else
207  {
208    poly pDivQ = p_PolyDiv(p, q, TRUE, r);
209    poly ppFactor = NULL; poly qqFactor = NULL;
210    poly theGcd = p_ExtGcdHelper(q, qqFactor, p, ppFactor, r);
211    pFactor = ppFactor;
212    qFactor = p_Add_q(qqFactor,
213                      p_Neg(p_Mult_q(pDivQ, p_Copy(ppFactor, r), r), r),
214                      r);
215    return theGcd;
216  }
217}
218
219
220/* assumes that p and q are univariate polynomials in r,
221   mentioning the same variable;
222   assumes a global monomial ordering in r;
223   assumes that not both p and q are NULL;
224   returns the gcd of p and q;
225   moreover, afterwards pFactor and qFactor contain appropriate
226   factors such that gcd(p, q) = p * pFactor + q * qFactor;
227   leaves p and q unmodified */
228poly p_ExtGcd(poly p, poly &pFactor, poly q, poly &qFactor, ring r)
229{
230  assume((p != NULL) || (q != NULL));
231  poly a = p; poly b = q; BOOLEAN aCorrespondsToP = TRUE;
232  if (p_Deg(a, r) < p_Deg(b, r))
233    { a = q; b = p; aCorrespondsToP = FALSE; }
234  a = p_Copy(a, r); b = p_Copy(b, r);
235  poly aFactor = NULL; poly bFactor = NULL;
236  poly theGcd = p_ExtGcdHelper(a, aFactor, b, bFactor, r);
237  if (aCorrespondsToP) { pFactor = aFactor; qFactor = bFactor; }
238  else                 { pFactor = bFactor; qFactor = aFactor; }
239  return theGcd;
240}
241
242
243
244#ifdef LDEBUG
245BOOLEAN naDBTest(number a, const char *f, const int l, const coeffs cf)
246{
247  if (a == NULL) return TRUE;
248  p_Test((poly)a, naRing);
249  if (getCoeffType(cf)==n_algExt)
250  {
251    if((((poly)a)!=naMinpoly)
252    && p_Totaldegree((poly)a, naRing) >= p_Totaldegree(naMinpoly, naRing)
253    && (p_Totaldegree((poly)a, naRing)> 1)) // allow to output par(1)
254    {
255      dReportError("deg >= deg(minpoly) in %s:%d\n",f,l);
256      return FALSE;
257    }
258  }
259  return TRUE;
260}
261#endif
262
263void heuristicReduce(poly &p, poly reducer, const coeffs cf);
264void definiteReduce(poly &p, poly reducer, const coeffs cf);
265
266/* returns the bottom field in this field extension tower; if the tower
267   is flat, i.e., if there is no extension, then r itself is returned;
268   as a side-effect, the counter 'height' is filled with the height of
269   the extension tower (in case the tower is flat, 'height' is zero) */
270static coeffs nCoeff_bottom(const coeffs r, int &height)
271{
272  assume(r != NULL);
273  coeffs cf = r;
274  height = 0;
275  while (nCoeff_is_Extension(cf))
276  {
277    assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL);
278    cf = cf->extRing->cf;
279    height++;
280  }
281  return cf;
282}
283
284BOOLEAN naIsZero(number a, const coeffs cf)
285{
286  naTest(a);
287  return (a == NULL);
288}
289
290void naDelete(number * a, const coeffs cf)
291{
292  if (*a == NULL) return;
293  if (((poly)*a)==naMinpoly) { *a=NULL;return;}
294  poly aAsPoly = (poly)(*a);
295  p_Delete(&aAsPoly, naRing);
296  *a = NULL;
297}
298
299BOOLEAN naEqual(number a, number b, const coeffs cf)
300{
301  naTest(a); naTest(b);
302  /// simple tests
303  if (a == NULL) return (b == NULL);
304  if (b == NULL) return (a == NULL);
305  return p_EqualPolys((poly)a,(poly)b,naRing);
306}
307
308number naCopy(number a, const coeffs cf)
309{
310  naTest(a);
311  if (a == NULL) return NULL;
312  if (((poly)a)==naMinpoly) return a;
313  return (number)p_Copy((poly)a, naRing);
314}
315
316number naGetNumerator(number &a, const coeffs cf)
317{
318  return naCopy(a, cf);
319}
320
321number naGetDenom(number &a, const coeffs cf)
322{
323  naTest(a);
324  return naInit(1, cf);
325}
326
327BOOLEAN naIsOne(number a, const coeffs cf)
328{
329  naTest(a);
330  poly aAsPoly = (poly)a;
331  if ((a==NULL) || (!p_IsConstant(aAsPoly, naRing))) return FALSE;
332  return n_IsOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
333}
334
335BOOLEAN naIsMOne(number a, const coeffs cf)
336{
337  naTest(a);
338  poly aAsPoly = (poly)a;
339  if ((a==NULL) || (!p_IsConstant(aAsPoly, naRing))) return FALSE;
340  return n_IsMOne(p_GetCoeff(aAsPoly, naRing), naCoeffs);
341}
342
343/// this is in-place, modifies a
344number naNeg(number a, const coeffs cf)
345{
346  naTest(a);
347  if (a != NULL) a = (number)p_Neg((poly)a, naRing);
348  return a;
349}
350
351number naInit(long i, const coeffs cf)
352{
353  if (i == 0) return NULL;
354  else        return (number)p_ISet(i, naRing);
355}
356
357int naInt(number &a, const coeffs cf)
358{
359  naTest(a);
360  poly aAsPoly = (poly)a;
361  if(aAsPoly == NULL)
362    return 0;
363  if (!p_IsConstant(aAsPoly, naRing))
364    return 0;
365  assume( aAsPoly != NULL );
366  return n_Int(p_GetCoeff(aAsPoly, naRing), naCoeffs);
367}
368
369/* TRUE iff (a != 0 and (b == 0 or deg(a) > deg(b) or (deg(a)==deg(b) && lc(a)>lc(b))) */
370BOOLEAN naGreater(number a, number b, const coeffs cf)
371{
372  naTest(a); naTest(b);
373  if (naIsZero(a, cf))
374  {
375    if (naIsZero(b, cf)) return FALSE;
376    return !n_GreaterZero(pGetCoeff((poly)b),cf);
377  }
378  if (naIsZero(b, cf))
379  {
380    return n_GreaterZero(pGetCoeff((poly)a),cf);
381  }
382  int aDeg = p_Totaldegree((poly)a, naRing);
383  int bDeg = p_Totaldegree((poly)b, naRing);
384  if (aDeg>bDeg) return TRUE;
385  if (aDeg<bDeg) return FALSE;
386  return n_Greater(pGetCoeff((poly)a),pGetCoeff((poly)b),naCoeffs);
387}
388
389/* TRUE iff a != 0 and (LC(a) > 0 or deg(a) > 0) */
390BOOLEAN naGreaterZero(number a, const coeffs cf)
391{
392  naTest(a);
393  if (a == NULL)                                            return FALSE;
394  if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE;
395  if (p_Totaldegree((poly)a, naRing) > 0)                   return TRUE;
396  return FALSE;
397}
398
399void naCoeffWrite(const coeffs cf, BOOLEAN details)
400{
401  assume( cf != NULL );
402
403  const ring A = cf->extRing;
404
405  assume( A != NULL );
406  assume( A->cf != NULL );
407
408  n_CoeffWrite(A->cf, details);
409
410//  rWrite(A);
411
412  const int P = rVar(A);
413  assume( P > 0 );
414
415  Print("//   %d parameter    : ", P);
416
417  for (int nop=0; nop < P; nop ++)
418    Print("%s ", rRingVar(nop, A));
419
420  PrintLn();
421
422  const ideal I = A->qideal;
423
424  assume( I != NULL );
425  assume( IDELEMS(I) == 1 );
426
427
428  if ( details )
429  {
430    PrintS("//   minpoly        : (");
431    p_Write0( I->m[0], A);
432    PrintS(")");
433  }
434  else
435    PrintS("//   minpoly        : ...");
436
437  PrintLn();
438
439/*
440  char *x = rRingVar(0, A);
441
442  Print("//   Coefficients live in the extension field K[%s]/<f(%s)>\n", x, x);
443  Print("//   with the minimal polynomial f(%s) = %s\n", x,
444        p_String(A->qideal->m[0], A));
445  PrintS("//   and K: ");
446*/
447}
448
449number naAdd(number a, number b, const coeffs cf)
450{
451  naTest(a); naTest(b);
452  if (a == NULL) return naCopy(b, cf);
453  if (b == NULL) return naCopy(a, cf);
454  poly aPlusB = p_Add_q(p_Copy((poly)a, naRing),
455                        p_Copy((poly)b, naRing), naRing);
456  //definiteReduce(aPlusB, naMinpoly, cf);
457  return (number)aPlusB;
458}
459
460number naSub(number a, number b, const coeffs cf)
461{
462  naTest(a); naTest(b);
463  if (b == NULL) return naCopy(a, cf);
464  poly minusB = p_Neg(p_Copy((poly)b, naRing), naRing);
465  if (a == NULL) return (number)minusB;
466  poly aMinusB = p_Add_q(p_Copy((poly)a, naRing), minusB, naRing);
467  //definiteReduce(aMinusB, naMinpoly, cf);
468  return (number)aMinusB;
469}
470
471number naMult(number a, number b, const coeffs cf)
472{
473  naTest(a); naTest(b);
474  if ((a == NULL)||(b == NULL)) return NULL;
475  poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing),
476                          p_Copy((poly)b, naRing), naRing);
477  definiteReduce(aTimesB, naMinpoly, cf);
478  p_Normalize(aTimesB,naRing);
479  return (number)aTimesB;
480}
481
482number naDiv(number a, number b, const coeffs cf)
483{
484  naTest(a); naTest(b);
485  if (b == NULL) WerrorS(nDivBy0);
486  if (a == NULL) return NULL;
487  poly bInverse = (poly)naInvers(b, cf);
488  if(bInverse != NULL) // b is non-zero divisor!
489  {
490    poly aDivB = p_Mult_q(p_Copy((poly)a, naRing), bInverse, naRing);
491    definiteReduce(aDivB, naMinpoly, cf);
492    p_Normalize(aDivB,naRing);
493    return (number)aDivB;
494  }
495  return NULL;
496}
497
498/* 0^0 = 0;
499   for |exp| <= 7 compute power by a simple multiplication loop;
500   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
501      p^13 = p^1 * p^4 * p^8, where we utilise that
502      p^(2^(k+1)) = p^(2^k) * p^(2^k);
503   intermediate reduction modulo the minimal polynomial is controlled by
504   the in-place method heuristicReduce(poly, poly, coeffs); see there.
505*/
506void naPower(number a, int exp, number *b, const coeffs cf)
507{
508  naTest(a);
509
510  /* special cases first */
511  if (a == NULL)
512  {
513    if (exp >= 0) *b = NULL;
514    else          WerrorS(nDivBy0);
515    return;
516  }
517  else if (exp ==  0) { *b = naInit(1, cf); return; }
518  else if (exp ==  1) { *b = naCopy(a, cf); return; }
519  else if (exp == -1) { *b = naInvers(a, cf); return; }
520
521  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
522
523  /* now compute a^expAbs */
524  poly pow; poly aAsPoly = (poly)a;
525  if (expAbs <= 7)
526  {
527    pow = p_Copy(aAsPoly, naRing);
528    for (int i = 2; i <= expAbs; i++)
529    {
530      pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing);
531      heuristicReduce(pow, naMinpoly, cf);
532    }
533    definiteReduce(pow, naMinpoly, cf);
534  }
535  else
536  {
537    pow = p_ISet(1, naRing);
538    poly factor = p_Copy(aAsPoly, naRing);
539    while (expAbs != 0)
540    {
541      if (expAbs & 1)
542      {
543        pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing);
544        heuristicReduce(pow, naMinpoly, cf);
545      }
546      expAbs = expAbs / 2;
547      if (expAbs != 0)
548      {
549        factor = p_Mult_q(factor, p_Copy(factor, naRing), naRing);
550        heuristicReduce(factor, naMinpoly, cf);
551      }
552    }
553    p_Delete(&factor, naRing);
554    definiteReduce(pow, naMinpoly, cf);
555  }
556
557  /* invert if original exponent was negative */
558  number n = (number)pow;
559  if (exp < 0)
560  {
561    number m = naInvers(n, cf);
562    naDelete(&n, cf);
563    n = m;
564  }
565  *b = n;
566}
567
568/* may reduce p modulo the reducer by calling definiteReduce;
569   the decision is made based on the following heuristic
570   (which should also only be changed here in this method):
571      if (deg(p) > 10*deg(reducer) then perform reduction;
572   modifies p */
573void heuristicReduce(poly &p, poly reducer, const coeffs cf)
574{
575  #ifdef LDEBUG
576  p_Test((poly)p, naRing);
577  p_Test((poly)reducer, naRing);
578  #endif
579  if (p_Totaldegree(p, naRing) > 10 * p_Totaldegree(reducer, naRing))
580    definiteReduce(p, reducer, cf);
581}
582
583void naWriteLong(number &a, const coeffs cf)
584{
585  naTest(a);
586  if (a == NULL)
587    StringAppendS("0");
588  else
589  {
590    poly aAsPoly = (poly)a;
591    /* basically, just write aAsPoly using p_Write,
592       but use brackets around the output, if a is not
593       a constant living in naCoeffs = cf->extRing->cf */
594    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
595    if (useBrackets) StringAppendS("(");
596    p_String0Long(aAsPoly, naRing, naRing);
597    if (useBrackets) StringAppendS(")");
598  }
599}
600
601void naWriteShort(number &a, const coeffs cf)
602{
603  naTest(a);
604  if (a == NULL)
605    StringAppendS("0");
606  else
607  {
608    poly aAsPoly = (poly)a;
609    /* basically, just write aAsPoly using p_Write,
610       but use brackets around the output, if a is not
611       a constant living in naCoeffs = cf->extRing->cf */
612    BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing));
613    if (useBrackets) StringAppendS("(");
614    p_String0Short(aAsPoly, naRing, naRing);
615    if (useBrackets) StringAppendS(")");
616  }
617}
618
619const char * naRead(const char *s, number *a, const coeffs cf)
620{
621  poly aAsPoly;
622  const char * result = p_Read(s, aAsPoly, naRing);
623  if (aAsPoly!=NULL) definiteReduce(aAsPoly, naMinpoly, cf);
624  *a = (number)aAsPoly;
625  return result;
626}
627
628#if 0
629/* implemented by the rule lcm(a, b) = a * b / gcd(a, b) */
630number naLcm(number a, number b, const coeffs cf)
631{
632  naTest(a); naTest(b);
633  if (a == NULL) return NULL;
634  if (b == NULL) return NULL;
635  number theProduct = (number)p_Mult_q(p_Copy((poly)a, naRing),
636                                       p_Copy((poly)b, naRing), naRing);
637  /* note that theProduct needs not be reduced w.r.t. naMinpoly;
638     but the final division will take care of the necessary reduction */
639  number theGcd = naGcd(a, b, cf);
640  return naDiv(theProduct, theGcd, cf);
641}
642#endif
643number napNormalizeHelper(number b, const coeffs cf)
644{
645  number h=n_Init(1,naRing->cf);
646  poly bb=(poly)b;
647  number d;
648  while(bb!=NULL)
649  {
650    d=n_NormalizeHelper(h,pGetCoeff(bb), naRing->cf);
651    n_Delete(&h,naRing->cf);
652    h=d;
653    pIter(bb);
654  }
655  return h;
656}
657number naLcmContent(number a, number b, const coeffs cf)
658{
659  if (nCoeff_is_Zp(naRing->cf)) return naCopy(a,cf);
660#if 0
661  else {
662    number g = ndGcd(a, b, cf);
663    return g;
664  }
665#else
666  {
667    a=(number)p_Copy((poly)a,naRing);
668    number t=napNormalizeHelper(b,cf);
669    if(!n_IsOne(t,naRing->cf))
670    {
671      number bt, rr;
672      poly xx=(poly)a;
673      while (xx!=NULL)
674      {
675        bt = n_SubringGcd(t, pGetCoeff(xx), naRing->cf);
676        rr = n_Mult(t, pGetCoeff(xx), naRing->cf);
677        n_Delete(&pGetCoeff(xx),naRing->cf);
678        pGetCoeff(xx) = n_Div(rr, bt, naRing->cf);
679        n_Normalize(pGetCoeff(xx),naRing->cf);
680        n_Delete(&bt,naRing->cf);
681        n_Delete(&rr,naRing->cf);
682        pIter(xx);
683      }
684    }
685    n_Delete(&t,naRing->cf);
686    return (number) a;
687  }
688#endif
689}
690
691/* expects *param to be castable to AlgExtInfo */
692static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
693{
694  if (n_algExt != n) return FALSE;
695  AlgExtInfo *e = (AlgExtInfo *)param;
696  /* for extension coefficient fields we expect the underlying
697     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
698     this expectation is based on the assumption that we have properly
699     registered cf and perform reference counting rather than creating
700     multiple copies of the same coefficient field/domain/ring */
701  if (naRing == e->r)
702    return TRUE;
703  /* (Note that then also the minimal ideals will necessarily be
704     the same, as they are attached to the ring.) */
705
706  // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)...
707  if( rEqual(naRing, e->r, TRUE) ) // also checks the equality of qideals
708  {
709    const ideal mi = naRing->qideal;
710    assume( IDELEMS(mi) == 1 );
711    const ideal ii = e->r->qideal;
712    assume( IDELEMS(ii) == 1 );
713
714    // TODO: the following should be extended for 2 *equal* rings...
715    assume( p_EqualPolys(mi->m[0], ii->m[0], naRing, e->r) );
716
717    rDelete(e->r);
718
719    return TRUE;
720  }
721
722  return FALSE;
723
724}
725
726int naSize(number a, const coeffs cf)
727{
728  if (a == NULL) return -1;
729  /* this has been taken from the old implementation of field extensions,
730     where we computed the sum of the degree and the number of terms in
731     (poly)a; so we leave it at that, for the time being;
732     maybe, the number of terms alone is a better measure? */
733  poly aAsPoly = (poly)a;
734  int theDegree = 0; int noOfTerms = 0;
735  while (aAsPoly != NULL)
736  {
737    noOfTerms++;
738    int d = p_GetExp(aAsPoly, 1, naRing);
739    if (d > theDegree) theDegree = d;
740    pIter(aAsPoly);
741  }
742  return theDegree + noOfTerms;
743}
744
745/* performs polynomial division and overrides p by the remainder
746   of division of p by the reducer;
747   modifies p */
748void definiteReduce(poly &p, poly reducer, const coeffs cf)
749{
750  #ifdef LDEBUG
751  p_Test((poly)p, naRing);
752  p_Test((poly)reducer, naRing);
753  #endif
754  if ((p!=NULL) && (p_GetExp(p,1,naRing)>=p_GetExp(reducer,1,naRing)))
755  {
756    p_PolyDiv(p, reducer, FALSE, naRing);
757  }
758}
759
760void  naNormalize(number &a, const coeffs cf)
761{
762  poly aa=(poly)a;
763  if (aa!=naMinpoly)
764    definiteReduce(aa,naMinpoly,cf);
765  a=(number)aa;
766}
767
768number naConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
769{
770  if (n.isZero()) return NULL;
771  poly p=convFactoryPSingP(n,naRing);
772  return (number)p;
773}
774CanonicalForm naConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
775{
776  naTest(n);
777  if (n==NULL) return CanonicalForm(0);
778
779  return convSingPFactoryP((poly)n,naRing);
780}
781
782/* IMPORTANT NOTE: Since an algebraic field extension is again a field,
783                   the gcd of two elements is not very interesting. (It
784                   is actually any unit in the field, i.e., any non-
785                   zero element.) Note that the below method does not operate
786                   in this strong sense but rather computes the gcd of
787                   two given elements in the underlying polynomial ring. */
788number naGcd(number a, number b, const coeffs cf)
789{
790  if (a==NULL)  return naCopy(b,cf);
791  if (b==NULL)  return naCopy(a,cf);
792
793  poly ax=(poly)a;
794  poly bx=(poly)b;
795  if (pNext(ax)!=NULL)
796    return (number)p_Copy(ax, naRing);
797  else
798  {
799    if(nCoeff_is_Zp(naRing->cf))
800      return naInit(1,cf);
801    else
802    {
803      number x = n_Copy(pGetCoeff((poly)a),naRing->cf);
804      if (n_IsOne(x,naRing->cf))
805        return (number)p_NSet(x,naRing);
806      while (pNext(ax)!=NULL)
807      {
808        pIter(ax);
809        number y = n_SubringGcd(x, pGetCoeff(ax), naRing->cf);
810        n_Delete(&x,naRing->cf);
811        x = y;
812        if (n_IsOne(x,naRing->cf))
813          return (number)p_NSet(x,naRing);
814      }
815      do
816      {
817        number y = n_SubringGcd(x, pGetCoeff(bx), naRing->cf);
818        n_Delete(&x,naRing->cf);
819        x = y;
820        if (n_IsOne(x,naRing->cf))
821          return (number)p_NSet(x,naRing);
822        pIter(bx);
823      }
824      while (bx!=NULL);
825      return (number)p_NSet(x,naRing);
826    }
827  }
828#if 0
829  naTest(a); naTest(b);
830  const ring R = naRing;
831  return (number) singclap_gcd(p_Copy((poly)a, R), p_Copy((poly)b, R), R);
832#endif
833//  return (number)p_Gcd((poly)a, (poly)b, naRing);
834}
835
836number naInvers(number a, const coeffs cf)
837{
838  naTest(a);
839  if (a == NULL) WerrorS(nDivBy0);
840
841  poly aFactor = NULL; poly mFactor = NULL; poly theGcd = NULL;
842// singclap_extgcd!
843  const BOOLEAN ret = singclap_extgcd ((poly)a, naMinpoly, theGcd, aFactor, mFactor, naRing);
844
845  assume( !ret );
846
847//  if( ret ) theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
848
849  naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor);
850  p_Delete(&mFactor, naRing);
851
852  //  /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */
853  //  assume(naIsOne((number)theGcd, cf));
854
855  if( !naIsOne((number)theGcd, cf) )
856  {
857    WerrorS("zero divisor found - your minpoly is not irreducible");
858    p_Delete(&aFactor, naRing); aFactor = NULL;
859  }
860  p_Delete(&theGcd, naRing);
861
862  return (number)(aFactor);
863}
864
865/* assumes that src = Q or Z, dst = Q(a) */
866number naMap00(number a, const coeffs src, const coeffs dst)
867{
868  if (n_IsZero(a, src)) return NULL;
869  assume(src->rep == dst->extRing->cf->rep);
870  poly result = p_One(dst->extRing);
871  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
872  return (number)result;
873}
874
875/* assumes that src = Z, dst = K(a) */
876number naMapZ0(number a, const coeffs src, const coeffs dst)
877{
878  if (n_IsZero(a, src)) return NULL;
879  poly result = p_One(dst->extRing);
880  nMapFunc nMap=n_SetMap(src,dst->extRing->cf);
881  p_SetCoeff(result, nMap(a, src, dst->extRing->cf), dst->extRing);
882  if (n_IsZero(pGetCoeff(result),dst->extRing->cf))
883    p_Delete(&result,dst->extRing);
884  return (number)result;
885}
886
887/* assumes that src = Z/p, dst = Q(a) */
888number naMapP0(number a, const coeffs src, const coeffs dst)
889{
890  if (n_IsZero(a, src)) return NULL;
891  /* mapping via intermediate int: */
892  int n = n_Int(a, src);
893  number q = n_Init(n, dst->extRing->cf);
894  poly result = p_One(dst->extRing);
895  p_SetCoeff(result, q, dst->extRing);
896  return (number)result;
897}
898
899#if 0
900/* assumes that either src = Q(a), dst = Q(a), or
901                       src = Z/p(a), dst = Z/p(a)     */
902number naCopyMap(number a, const coeffs src, const coeffs dst)
903{
904  return naCopy(a, dst);
905}
906#endif
907
908number naCopyTrans2AlgExt(number a, const coeffs src, const coeffs dst)
909{
910  assume (nCoeff_is_transExt (src));
911  assume (nCoeff_is_algExt (dst));
912  fraction fa=(fraction)a;
913  poly p, q;
914  if (rSamePolyRep(src->extRing, dst->extRing))
915  {
916    p = p_Copy(NUM(fa),src->extRing);
917    if (!DENIS1(fa))
918    {
919      q = p_Copy(DEN(fa),src->extRing);
920      assume (q != NULL);
921    }
922  }
923  else
924  {
925    assume ((strcmp(rRingVar(0,src->extRing),rRingVar(0,dst->extRing))==0) && (rVar (src->extRing) == rVar (dst->extRing)));
926
927    nMapFunc nMap= n_SetMap (src->extRing->cf, dst->extRing->cf);
928
929    assume (nMap != NULL);
930    p= p_PermPoly (NUM (fa), NULL, src->extRing, dst->extRing,nMap, NULL,rVar (src->extRing));
931    if (!DENIS1(fa))
932    {
933      q= p_PermPoly (DEN (fa), NULL, src->extRing, dst->extRing,nMap, NULL,rVar (src->extRing));
934      assume (q != NULL);
935    }
936  }
937  definiteReduce(p, dst->extRing->qideal->m[0], dst);
938  p_Test (p, dst->extRing);
939  if (!DENIS1(fa))
940  {
941    definiteReduce(q, dst->extRing->qideal->m[0], dst);
942    p_Test (q, dst->extRing);
943    if (q != NULL)
944    {
945      number t= naDiv ((number)p,(number)q, dst);
946      p_Delete (&p, dst->extRing);
947      p_Delete (&q, dst->extRing);
948      return t;
949    }
950    WerrorS ("mapping denominator to zero");
951  }
952  return (number) p;
953}
954
955/* assumes that src = Q, dst = Z/p(a) */
956number naMap0P(number a, const coeffs src, const coeffs dst)
957{
958  if (n_IsZero(a, src)) return NULL;
959  // int p = rChar(dst->extRing);
960
961  number q = nlModP(a, src, dst->extRing->cf);
962
963  poly result = p_NSet(q, dst->extRing);
964
965  return (number)result;
966}
967
968/* assumes that src = Z/p, dst = Z/p(a) */
969number naMapPP(number a, const coeffs src, const coeffs dst)
970{
971  if (n_IsZero(a, src)) return NULL;
972  assume(src == dst->extRing->cf);
973  poly result = p_One(dst->extRing);
974  p_SetCoeff(result, n_Copy(a, src), dst->extRing);
975  return (number)result;
976}
977
978/* assumes that src = Z/u, dst = Z/p(a), where u != p */
979number naMapUP(number a, const coeffs src, const coeffs dst)
980{
981  if (n_IsZero(a, src)) return NULL;
982  /* mapping via intermediate int: */
983  int n = n_Int(a, src);
984  number q = n_Init(n, dst->extRing->cf);
985  poly result = p_One(dst->extRing);
986  p_SetCoeff(result, q, dst->extRing);
987  return (number)result;
988}
989
990number naGenMap(number a, const coeffs cf, const coeffs dst)
991{
992  if (a==NULL) return NULL;
993
994  const ring rSrc = cf->extRing;
995  const ring rDst = dst->extRing;
996
997  const nMapFunc nMap=n_SetMap(rSrc->cf,rDst->cf);
998  poly f = (poly)a;
999  poly g = prMapR(f, nMap, rSrc, rDst);
1000
1001  n_Test((number)g, dst);
1002  return (number)g;
1003}
1004
1005number naGenTrans2AlgExt(number a, const coeffs cf, const coeffs dst)
1006{
1007  if (a==NULL) return NULL;
1008
1009  const ring rSrc = cf->extRing;
1010  const ring rDst = dst->extRing;
1011
1012  const nMapFunc nMap=n_SetMap(rSrc->cf,rDst->cf);
1013  fraction f = (fraction)a;
1014  poly g = prMapR(NUM(f), nMap, rSrc, rDst);
1015
1016  number result=NULL;
1017  poly h = NULL;
1018
1019  if (!DENIS1(f))
1020     h = prMapR(DEN(f), nMap, rSrc, rDst);
1021
1022  if (h!=NULL)
1023  {
1024    result=naDiv((number)g,(number)h,dst);
1025    p_Delete(&g,dst->extRing);
1026    p_Delete(&h,dst->extRing);
1027  }
1028  else
1029    result=(number)g;
1030
1031  n_Test((number)result, dst);
1032  return (number)result;
1033}
1034
1035nMapFunc naSetMap(const coeffs src, const coeffs dst)
1036{
1037  /* dst is expected to be an algebraic field extension */
1038  assume(getCoeffType(dst) == ID);
1039
1040  if( src == dst ) return ndCopyMap;
1041
1042  int h = 0; /* the height of the extension tower given by dst */
1043  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
1044  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
1045
1046  /* for the time being, we only provide maps if h = 1 or 0 */
1047  if (h==0)
1048  {
1049    if ((src->rep==n_rep_gap_rat) && nCoeff_is_Q(bDst))
1050      return naMap00;                            /// Q or Z     -->  Q(a)
1051    if ((src->rep==n_rep_gap_gmp) && nCoeff_is_Q(bDst))
1052      return naMapZ0;                            /// Z     -->  Q(a)
1053    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
1054      return naMapP0;                            /// Z/p   -->  Q(a)
1055    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
1056      return naMap0P;                            /// Q      --> Z/p(a)
1057    if ((src->rep==n_rep_gap_gmp) && nCoeff_is_Zp(bDst))
1058      return naMapZ0;                            /// Z     -->  Z/p(a)
1059    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
1060    {
1061      if (src->ch == dst->ch) return naMapPP;    /// Z/p    --> Z/p(a)
1062      else return naMapUP;                       /// Z/u    --> Z/p(a)
1063    }
1064  }
1065  if (h != 1) return NULL;
1066  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
1067  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
1068
1069  nMapFunc nMap=n_SetMap(src->extRing->cf,dst->extRing->cf);
1070  if (rSamePolyRep(src->extRing, dst->extRing) && (strcmp(rRingVar(0, src->extRing), rRingVar(0, dst->extRing)) == 0))
1071  {
1072    if (src->type==n_algExt)
1073       return ndCopyMap; // naCopyMap;         /// K(a)   --> K(a)
1074    else
1075       return naCopyTrans2AlgExt;
1076  }
1077  else if ((nMap!=NULL) && (strcmp(rRingVar(0,src->extRing),rRingVar(0,dst->extRing))==0) && (rVar (src->extRing) == rVar (dst->extRing)))
1078  {
1079    if (src->type==n_algExt)
1080       return naGenMap; // naCopyMap;         /// K(a)   --> K'(a)
1081    else
1082       return naGenTrans2AlgExt;
1083  }
1084
1085  return NULL;                                           /// default
1086}
1087
1088static int naParDeg(number a, const coeffs cf)
1089{
1090  if (a == NULL) return -1;
1091  poly aa=(poly)a;
1092  return cf->extRing->pFDeg(aa,cf->extRing);
1093}
1094
1095/// return the specified parameter as a number in the given alg. field
1096static number naParameter(const int iParameter, const coeffs cf)
1097{
1098  assume(getCoeffType(cf) == ID);
1099
1100  const ring R = cf->extRing;
1101  assume( R != NULL );
1102  assume( 0 < iParameter && iParameter <= rVar(R) );
1103
1104  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
1105
1106  return (number) p;
1107}
1108
1109
1110/// if m == var(i)/1 => return i,
1111int naIsParam(number m, const coeffs cf)
1112{
1113  assume(getCoeffType(cf) == ID);
1114
1115  const ring R = cf->extRing;
1116  assume( R != NULL );
1117
1118  return p_Var( (poly)m, R );
1119}
1120
1121
1122static void naClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1123{
1124  assume(cf != NULL);
1125  assume(getCoeffType(cf) == ID);
1126  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
1127
1128  const ring   R = cf->extRing;
1129  assume(R != NULL);
1130  const coeffs Q = R->cf;
1131  assume(Q != NULL);
1132  assume(nCoeff_is_Q(Q));
1133
1134  numberCollectionEnumerator.Reset();
1135
1136  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
1137  {
1138    c = n_Init(1, cf);
1139    return;
1140  }
1141
1142  naTest(numberCollectionEnumerator.Current());
1143
1144  // part 1, find a small candidate for gcd
1145  int s1; int s=2147483647; // max. int
1146
1147  const BOOLEAN lc_is_pos=naGreaterZero(numberCollectionEnumerator.Current(),cf);
1148
1149  int normalcount = 0;
1150
1151  poly cand1, cand;
1152
1153  do
1154  {
1155    number& n = numberCollectionEnumerator.Current();
1156    naNormalize(n, cf); ++normalcount;
1157
1158    naTest(n);
1159
1160    cand1 = (poly)n;
1161
1162    s1 = p_Deg(cand1, R); // naSize?
1163    if (s>s1)
1164    {
1165      cand = cand1;
1166      s = s1;
1167    }
1168  } while (numberCollectionEnumerator.MoveNext() );
1169
1170//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
1171
1172  cand = p_Copy(cand, R);
1173  // part 2: compute gcd(cand,all coeffs)
1174
1175  numberCollectionEnumerator.Reset();
1176
1177  int length = 0;
1178  while (numberCollectionEnumerator.MoveNext() )
1179  {
1180    number& n = numberCollectionEnumerator.Current();
1181    ++length;
1182
1183    if( (--normalcount) <= 0)
1184      naNormalize(n, cf);
1185
1186    naTest(n);
1187
1188//    p_InpGcd(cand, (poly)n, R);
1189
1190    cand = singclap_gcd(cand, p_Copy((poly)n, R), R);
1191
1192//    cand1 = p_Gcd(cand,(poly)n, R); p_Delete(&cand, R); cand = cand1;
1193
1194    assume( naGreaterZero((number)cand, cf) ); // ???
1195/*
1196    if(p_IsConstant(cand,R))
1197    {
1198      c = cand;
1199
1200      if(!lc_is_pos)
1201      {
1202        // make the leading coeff positive
1203        c = nlNeg(c, cf);
1204        numberCollectionEnumerator.Reset();
1205
1206        while (numberCollectionEnumerator.MoveNext() )
1207        {
1208          number& nn = numberCollectionEnumerator.Current();
1209          nn = nlNeg(nn, cf);
1210        }
1211      }
1212      return;
1213    }
1214*/
1215
1216  }
1217
1218
1219  // part3: all coeffs = all coeffs / cand
1220  if (!lc_is_pos)
1221    cand = p_Neg(cand, R);
1222
1223  c = (number)cand; naTest(c);
1224
1225  poly cInverse = (poly)naInvers(c, cf);
1226  assume(cInverse != NULL); // c is non-zero divisor!?
1227
1228
1229  numberCollectionEnumerator.Reset();
1230
1231
1232  while (numberCollectionEnumerator.MoveNext() )
1233  {
1234    number& n = numberCollectionEnumerator.Current();
1235
1236    assume( length > 0 );
1237
1238    if( --length > 0 )
1239    {
1240      assume( cInverse != NULL );
1241      n = (number) p_Mult_q(p_Copy(cInverse, R), (poly)n, R);
1242    }
1243    else
1244    {
1245      n = (number) p_Mult_q(cInverse, (poly)n, R);
1246      cInverse = NULL;
1247      assume(length == 0);
1248    }
1249
1250    definiteReduce((poly &)n, naMinpoly, cf);
1251  }
1252
1253  assume(length == 0);
1254  assume(cInverse == NULL); //   p_Delete(&cInverse, R);
1255
1256  // Quick and dirty fix for constant content clearing... !?
1257  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
1258
1259  number cc;
1260
1261  n_ClearContent(itr, cc, Q); // TODO: get rid of (-LC) normalization!?
1262
1263  // over alg. ext. of Q // takes over the input number
1264  c = (number) p_Mult_nn( (poly)c, cc, R);
1265//      p_Mult_q(p_NSet(cc, R), , R);
1266
1267  n_Delete(&cc, Q);
1268
1269  // TODO: the above is not enough! need GCD's of polynomial coeffs...!
1270/*
1271  // old and wrong part of p_Content
1272    if (rField_is_Q_a(r) && !CLEARENUMERATORS) // should not be used anymore if CLEARENUMERATORS is 1
1273    {
1274      // we only need special handling for alg. ext.
1275      if (getCoeffType(r->cf)==n_algExt)
1276      {
1277        number hzz = n_Init(1, r->cf->extRing->cf);
1278        p=ph;
1279        while (p!=NULL)
1280        { // each monom: coeff in Q_a
1281          poly c_n_n=(poly)pGetCoeff(p);
1282          poly c_n=c_n_n;
1283          while (c_n!=NULL)
1284          { // each monom: coeff in Q
1285            d=n_NormalizeHelper(hzz,pGetCoeff(c_n),r->cf->extRing->cf);
1286            n_Delete(&hzz,r->cf->extRing->cf);
1287            hzz=d;
1288            pIter(c_n);
1289          }
1290          pIter(p);
1291        }
1292        // hzz contains the 1/lcm of all denominators in c_n_n
1293        h=n_Invers(hzz,r->cf->extRing->cf);
1294        n_Delete(&hzz,r->cf->extRing->cf);
1295        n_Normalize(h,r->cf->extRing->cf);
1296        if(!n_IsOne(h,r->cf->extRing->cf))
1297        {
1298          p=ph;
1299          while (p!=NULL)
1300          { // each monom: coeff in Q_a
1301            poly c_n=(poly)pGetCoeff(p);
1302            while (c_n!=NULL)
1303            { // each monom: coeff in Q
1304              d=n_Mult(h,pGetCoeff(c_n),r->cf->extRing->cf);
1305              n_Normalize(d,r->cf->extRing->cf);
1306              n_Delete(&pGetCoeff(c_n),r->cf->extRing->cf);
1307              pGetCoeff(c_n)=d;
1308              pIter(c_n);
1309            }
1310            pIter(p);
1311          }
1312        }
1313        n_Delete(&h,r->cf->extRing->cf);
1314      }
1315    }
1316*/
1317
1318
1319//  c = n_Init(1, cf); assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
1320}
1321
1322
1323static void naClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1324{
1325  assume(cf != NULL);
1326  assume(getCoeffType(cf) == ID);
1327  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
1328
1329  assume(cf->extRing != NULL);
1330  const coeffs Q = cf->extRing->cf;
1331  assume(Q != NULL);
1332  assume(nCoeff_is_Q(Q));
1333  number n;
1334  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
1335  n_ClearDenominators(itr, n, Q); // this should probably be fine...
1336  c = (number)p_NSet(n, cf->extRing); // over alg. ext. of Q // takes over the input number
1337}
1338
1339void naKillChar(coeffs cf)
1340{
1341   if ((--cf->extRing->ref) == 0)
1342     rDelete(cf->extRing);
1343}
1344
1345char* naCoeffString(const coeffs r) // currently also for tranext.
1346{
1347  const char* const* p=n_ParameterNames(r);
1348  int l=0;
1349  int i;
1350  for(i=0; i<n_NumberOfParameters(r);i++)
1351  {
1352    l+=(strlen(p[i])+1);
1353  }
1354  char *s=(char *)omAlloc(l+10+1);
1355  s[0]='\0';
1356  snprintf(s,10+1,"%d",r->ch); /* Fp(a) or Q(a) */
1357  char tt[2];
1358  tt[0]=',';
1359  tt[1]='\0';
1360  for(i=0; i<n_NumberOfParameters(r);i++)
1361  {
1362    strcat(s,tt);
1363    strcat(s,p[i]);
1364  }
1365  return s;
1366}
1367
1368number  naChineseRemainder(number *x, number *q,int rl, BOOLEAN sym,const coeffs cf)
1369{
1370  poly *P=(poly*)omAlloc(rl*sizeof(poly*));
1371  number *X=(number *)omAlloc(rl*sizeof(number));
1372  int i;
1373  for(i=0;i<rl;i++) P[i]=p_Copy((poly)(x[i]),cf->extRing);
1374  poly result=p_ChineseRemainder(P,X,q,rl,cf->extRing);
1375  omFreeSize(X,rl*sizeof(number));
1376  omFreeSize(P,rl*sizeof(poly*));
1377  return ((number)result);
1378}
1379
1380number  naFarey(number p, number n, const coeffs cf)
1381{
1382  // n is really a bigint
1383  poly result=p_Farey(p_Copy((poly)p,cf->extRing),n,cf->extRing);
1384  return ((number)result);
1385}
1386
1387
1388BOOLEAN naInitChar(coeffs cf, void * infoStruct)
1389{
1390  assume( infoStruct != NULL );
1391
1392  AlgExtInfo *e = (AlgExtInfo *)infoStruct;
1393  /// first check whether cf->extRing != NULL and delete old ring???
1394
1395  assume(e->r                     != NULL);      // extRing;
1396  assume(e->r->cf                 != NULL);      // extRing->cf;
1397
1398  assume((e->r->qideal            != NULL) &&    // minideal has one
1399         (IDELEMS(e->r->qideal)   == 1)    &&    // non-zero generator
1400         (e->r->qideal->m[0]      != NULL)    ); // at m[0];
1401
1402  assume( cf != NULL );
1403  assume(getCoeffType(cf) == ID);                     // coeff type;
1404
1405  e->r->ref ++; // increase the ref.counter for the ground poly. ring!
1406  const ring R = e->r; // no copy!
1407  cf->extRing  = R;
1408
1409  /* propagate characteristic up so that it becomes
1410     directly accessible in cf: */
1411  cf->ch = R->cf->ch;
1412
1413  cf->is_field=TRUE;
1414  cf->is_domain=TRUE;
1415  cf->rep=n_rep_poly;
1416
1417  #ifdef LDEBUG
1418  p_Test((poly)naMinpoly, naRing);
1419  #endif
1420
1421  cf->cfCoeffString = naCoeffString;
1422
1423  cf->cfGreaterZero  = naGreaterZero;
1424  cf->cfGreater      = naGreater;
1425  cf->cfEqual        = naEqual;
1426  cf->cfIsZero       = naIsZero;
1427  cf->cfIsOne        = naIsOne;
1428  cf->cfIsMOne       = naIsMOne;
1429  cf->cfInit         = naInit;
1430  cf->cfFarey        = naFarey;
1431  cf->cfChineseRemainder= naChineseRemainder;
1432  cf->cfInt          = naInt;
1433  cf->cfInpNeg       = naNeg;
1434  cf->cfAdd          = naAdd;
1435  cf->cfSub          = naSub;
1436  cf->cfMult         = naMult;
1437  cf->cfDiv          = naDiv;
1438  cf->cfExactDiv     = naDiv;
1439  cf->cfPower        = naPower;
1440  cf->cfCopy         = naCopy;
1441
1442  cf->cfWriteLong        = naWriteLong;
1443
1444  if( rCanShortOut(naRing) )
1445    cf->cfWriteShort = naWriteShort;
1446  else
1447    cf->cfWriteShort = naWriteLong;
1448
1449  cf->cfRead         = naRead;
1450  cf->cfDelete       = naDelete;
1451  cf->cfSetMap       = naSetMap;
1452  cf->cfGetDenom     = naGetDenom;
1453  cf->cfGetNumerator = naGetNumerator;
1454  cf->cfRePart       = naCopy;
1455  cf->cfCoeffWrite   = naCoeffWrite;
1456  cf->cfNormalize    = naNormalize;
1457  cf->cfKillChar     = naKillChar;
1458#ifdef LDEBUG
1459  cf->cfDBTest       = naDBTest;
1460#endif
1461  cf->cfGcd          = naGcd;
1462  cf->cfNormalizeHelper          = naLcmContent;
1463  cf->cfSize         = naSize;
1464  cf->nCoeffIsEqual  = naCoeffIsEqual;
1465  cf->cfInvers       = naInvers;
1466  cf->convFactoryNSingN=naConvFactoryNSingN;
1467  cf->convSingNFactoryN=naConvSingNFactoryN;
1468  cf->cfParDeg = naParDeg;
1469
1470  cf->iNumberOfParameters = rVar(R);
1471  cf->pParameterNames = (const char**)R->names;
1472  cf->cfParameter = naParameter;
1473  cf->has_simple_Inverse= R->cf->has_simple_Inverse;
1474  /* cf->has_simple_Alloc= FALSE; */
1475
1476  if( nCoeff_is_Q(R->cf) )
1477  {
1478    cf->cfClearContent = naClearContent;
1479    cf->cfClearDenominators = naClearDenominators;
1480  }
1481
1482  return FALSE;
1483}
1484
1485template class CRecursivePolyCoeffsEnumerator<NAConverter>;
1486
1487template class IAccessor<snumber*>;
1488
1489/* --------------------------------------------------------------------*/
1490#if 0
1491void npolyCoeffWrite(const coeffs cf, BOOLEAN details)
1492{
1493  assume( cf != NULL );
1494
1495  const ring A = cf->extRing;
1496
1497  assume( A != NULL );
1498  Print("// polynomial ring as coefficient ring :\n");
1499  rWrite(A);
1500  PrintLn();
1501}
1502number npolyMult(number a, number b, const coeffs cf)
1503{
1504  naTest(a); naTest(b);
1505  if ((a == NULL)||(b == NULL)) return NULL;
1506  poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing),
1507                          p_Copy((poly)b, naRing), naRing);
1508  return (number)aTimesB;
1509}
1510
1511void npolyPower(number a, int exp, number *b, const coeffs cf)
1512{
1513  naTest(a);
1514
1515  /* special cases first */
1516  if (a == NULL)
1517  {
1518    if (exp >= 0) *b = NULL;
1519    else          WerrorS(nDivBy0);
1520    return;
1521  }
1522  else if (exp ==  0) { *b = naInit(1, cf); return; }
1523  else if (exp ==  1) { *b = naCopy(a, cf); return; }
1524  else if (exp == -1) { *b = naInvers(a, cf); return; }
1525
1526  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
1527
1528  /* now compute a^expAbs */
1529  poly pow; poly aAsPoly = (poly)a;
1530  if (expAbs <= 7)
1531  {
1532    pow = p_Copy(aAsPoly, naRing);
1533    for (int i = 2; i <= expAbs; i++)
1534    {
1535      pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing);
1536    }
1537  }
1538  else
1539  {
1540    pow = p_ISet(1, naRing);
1541    poly factor = p_Copy(aAsPoly, naRing);
1542    while (expAbs != 0)
1543    {
1544      if (expAbs & 1)
1545      {
1546        pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing);
1547      }
1548      expAbs = expAbs / 2;
1549      if (expAbs != 0)
1550      {
1551        factor = p_Mult_q(factor, p_Copy(factor, naRing), naRing);
1552      }
1553    }
1554    p_Delete(&factor, naRing);
1555  }
1556
1557  /* invert if original exponent was negative */
1558  number n = (number)pow;
1559  if (exp < 0)
1560  {
1561    number m = npolyInvers(n, cf);
1562    naDelete(&n, cf);
1563    n = m;
1564  }
1565  *b = n;
1566}
1567
1568number npolyDiv(number a, number b, const coeffs cf)
1569{
1570  naTest(a); naTest(b);
1571  if (b == NULL) WerrorS(nDivBy0);
1572  if (a == NULL) return NULL;
1573  poly p=singclap_pdivide((poly)a,(poly)b,naRing);
1574  return (number)p;
1575}
1576
1577
1578BOOLEAN npolyInitChar(coeffs cf, void * infoStruct)
1579{
1580  assume( infoStruct != NULL );
1581
1582  AlgExtInfo *e = (AlgExtInfo *)infoStruct;
1583  /// first check whether cf->extRing != NULL and delete old ring???
1584
1585  assume(e->r                     != NULL);      // extRing;
1586  assume(e->r->cf                 != NULL);      // extRing->cf;
1587
1588  assume( cf != NULL );
1589
1590  e->r->ref ++; // increase the ref.counter for the ground poly. ring!
1591  const ring R = e->r; // no copy!
1592  cf->extRing  = R;
1593
1594  /* propagate characteristic up so that it becomes
1595     directly accessible in cf: */
1596  cf->ch = R->cf->ch;
1597  cf->is_field=FALSE;
1598  cf->is_domain=TRUE;
1599
1600  cf->cfCoeffString = naCoeffString;
1601
1602  cf->cfGreaterZero  = naGreaterZero;
1603  cf->cfGreater      = naGreater;
1604  cf->cfEqual        = naEqual;
1605  cf->cfIsZero       = naIsZero;
1606  cf->cfIsOne        = naIsOne;
1607  cf->cfIsMOne       = naIsMOne;
1608  cf->cfInit         = naInit;
1609  cf->cfFarey        = naFarey;
1610  cf->cfChineseRemainder= naChineseRemainder;
1611  cf->cfInt          = naInt;
1612  cf->cfInpNeg       = naNeg;
1613  cf->cfAdd          = naAdd;
1614  cf->cfSub          = naSub;
1615  cf->cfMult         = npolyMult;
1616  cf->cfDiv          = npolyDiv;
1617  cf->cfPower        = naPower;
1618  cf->cfCopy         = naCopy;
1619
1620  cf->cfWriteLong        = naWriteLong;
1621
1622  if( rCanShortOut(naRing) )
1623    cf->cfWriteShort = naWriteShort;
1624  else
1625    cf->cfWriteShort = naWriteLong;
1626
1627  cf->cfRead         = naRead;
1628  cf->cfDelete       = naDelete;
1629  cf->cfSetMap       = naSetMap;
1630  cf->cfGetDenom     = naGetDenom;
1631  cf->cfGetNumerator = naGetNumerator;
1632  cf->cfRePart       = naCopy;
1633  cf->cfCoeffWrite   = npolyCoeffWrite;
1634  cf->cfNormalize    = npolyNormalize;
1635  cf->cfKillChar     = naKillChar;
1636#ifdef LDEBUG
1637  cf->cfDBTest       = naDBTest;
1638#endif
1639  cf->cfGcd          = naGcd;
1640  cf->cfNormalizeHelper          = naLcmContent;
1641  cf->cfSize         = naSize;
1642  cf->nCoeffIsEqual  = naCoeffIsEqual;
1643  cf->cfInvers       = npolyInvers;
1644  cf->convFactoryNSingN=naConvFactoryNSingN;
1645  cf->convSingNFactoryN=naConvSingNFactoryN;
1646  cf->cfParDeg = naParDeg;
1647
1648  cf->iNumberOfParameters = rVar(R);
1649  cf->pParameterNames = (const char**)R->names;
1650  cf->cfParameter = naParameter;
1651  cf->has_simple_Inverse=FALSE;
1652  /* cf->has_simple_Alloc= FALSE; */
1653
1654  if( nCoeff_is_Q(R->cf) )
1655  {
1656    cf->cfClearContent = naClearContent;
1657    cf->cfClearDenominators = naClearDenominators;
1658  }
1659
1660  return FALSE;
1661}
1662#endif
Note: See TracBrowser for help on using the repository browser.