source: git/libpolys/polys/ext_fields/algext.cc @ 8d1432e

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