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

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