source: git/libpolys/polys/ext_fields/transext.cc @ 603ac7a

jengelh-datetimespielwiese
Last change on this file since 603ac7a was 603ac7a, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: add P_Normalize to definiteGcdCancellation/ntNormalize
  • Property mode set to 100644
File size: 56.3 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: numbers in a rational function field K(t_1, .., t_s) with
6*           transcendental variables t_1, ..., t_s, where s >= 1.
7*           Denoting the implemented coeffs object by cf, then these numbers
8*           are represented as quotients of polynomials living in the
9*           polynomial ring K[t_1, .., t_s] represented by cf->extring.
10*
11*           An element of K(t_1, .., t_s) may have numerous representations,
12*           due to the possibility of common polynomial factors in the
13*           numerator and denominator. This problem is handled by a
14*           cancellation heuristic: Each number "knows" its complexity
15*           which is 0 if and only if common factors have definitely been
16*           cancelled, and some positive integer otherwise.
17*           Each arithmetic operation of two numbers with complexities c1
18*           and c2 will result in a number of complexity c1 + c2 + some
19*           penalty (specific for each arithmetic operation; see constants
20*           in the *.h file). Whenever the resulting complexity exceeds a
21*           certain threshold (see constant in the *.h file), then the
22*           cancellation heuristic will call 'factory' to compute the gcd
23*           and cancel it out in the given number. (This definite cancel-
24*           lation will also be performed at the beginning of ntWrite,
25*           ensuring that any output is free of common factors.
26*           For the special case of K = Q (i.e., when computing over the
27*           rationals), this definite cancellation procedure will also take
28*           care of nested fractions: If there are fractional coefficients
29*           in the numerator or denominator of a number, then this number
30*           is being replaced by a quotient of two polynomials over Z, or
31*           - if the denominator is a constant - by a polynomial over Q.
32*
33*           TODO: the description above needs a major update!!!
34*/
35#define TRANSEXT_PRIVATES
36
37#include "config.h"
38#include <misc/auxiliary.h>
39
40#include <omalloc/omalloc.h>
41
42#include <reporter/reporter.h>
43
44#include <coeffs/coeffs.h>
45#include <coeffs/numbers.h>
46#include <coeffs/longrat.h>
47
48#include <polys/monomials/ring.h>
49#include <polys/monomials/p_polys.h>
50#include <polys/simpleideals.h>
51
52#ifdef HAVE_FACTORY
53#include <polys/clapsing.h>
54#include <polys/clapconv.h>
55#include <factory/factory.h>
56#endif
57
58#include <polys/ext_fields/transext.h>
59#include <polys/prCopy.h>
60
61#include <polys/PolyEnumerator.h>
62
63
64/* constants for controlling the complexity of numbers */
65#define ADD_COMPLEXITY 1   /**< complexity increase due to + and - */
66#define MULT_COMPLEXITY 2   /**< complexity increase due to * and / */
67#define BOUND_COMPLEXITY 10   /**< maximum complexity of a number */
68
69/// TRUE iff num. represents 1
70#define NUMIS1(f) (p_IsOne(NUM(f), cf->extRing))
71
72#define COM(f) f->complexity
73
74
75#ifdef LDEBUG
76#define ntTest(a) assume(ntDBTest(a,__FILE__,__LINE__,cf))
77BOOLEAN  ntDBTest(number a, const char *f, const int l, const coeffs r);
78#else
79#define ntTest(a) do {} while (0)
80#endif
81
82/// Our own type!
83static const n_coeffType ID = n_transExt;
84
85/* polynomial ring in which the numerators and denominators of our
86   numbers live */
87#define ntRing cf->extRing
88
89/* coeffs object in which the coefficients of our numbers live;
90 * methods attached to ntCoeffs may be used to compute with the
91 * coefficients of our numbers, e.g., use ntCoeffs->nAdd to add
92 * coefficients of our numbers */
93#define ntCoeffs cf->extRing->cf
94
95
96omBin fractionObjectBin = omGetSpecBin(sizeof(fractionObject));
97
98/// forward declarations
99BOOLEAN  ntGreaterZero(number a, const coeffs cf);
100BOOLEAN  ntGreater(number a, number b, const coeffs cf);
101BOOLEAN  ntEqual(number a, number b, const coeffs cf);
102BOOLEAN  ntIsOne(number a, const coeffs cf);
103BOOLEAN  ntIsMOne(number a, const coeffs cf);
104BOOLEAN  ntIsZero(number a, const coeffs cf);
105number   ntInit(long i, const coeffs cf);
106int      ntInt(number &a, const coeffs cf);
107number   ntNeg(number a, const coeffs cf);
108number   ntInvers(number a, const coeffs cf);
109number   ntAdd(number a, number b, const coeffs cf);
110number   ntSub(number a, number b, const coeffs cf);
111number   ntMult(number a, number b, const coeffs cf);
112number   ntDiv(number a, number b, const coeffs cf);
113void     ntPower(number a, int exp, number *b, const coeffs cf);
114number   ntCopy(number a, const coeffs cf);
115void     ntWriteLong(number &a, const coeffs cf);
116void     ntWriteShort(number &a, const coeffs cf);
117number   ntRePart(number a, const coeffs cf);
118number   ntImPart(number a, const coeffs cf);
119number   ntGetDenom(number &a, const coeffs cf);
120number   ntGetNumerator(number &a, const coeffs cf);
121number   ntGcd(number a, number b, const coeffs cf);
122number   ntLcm(number a, number b, const coeffs cf);
123int      ntSize(number a, const coeffs cf);
124void     ntDelete(number * a, const coeffs cf);
125void     ntCoeffWrite(const coeffs cf, BOOLEAN details);
126number   ntIntDiv(number a, number b, const coeffs cf);
127const char * ntRead(const char *s, number *a, const coeffs cf);
128static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param);
129
130void heuristicGcdCancellation(number a, const coeffs cf);
131void definiteGcdCancellation(number a, const coeffs cf,
132                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed);
133void handleNestedFractionsOverQ(fraction f, const coeffs cf);
134
135#ifdef LDEBUG
136BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs cf)
137{
138  assume(getCoeffType(cf) == ID);
139
140  if (IS0(a)) return TRUE;
141
142  const fraction t = (fraction)a;
143
144  const poly num = NUM(t);
145  assume(num != NULL);   /**< t != 0 ==> numerator(t) != 0 */
146  assume( _p_Test(num, ntRing,1) );
147
148  const poly den = DEN(t);
149
150  if (den != NULL) // !DENIS1(f)
151  {
152    assume( _p_Test(den, ntRing,1) );
153
154    if(p_IsConstant(den, ntRing) && (n_IsOne(pGetCoeff(den), ntRing->cf)))
155    {
156      Print("?/1 in %s:%d\n",f,l);
157      return FALSE;
158    }
159
160    if( !n_GreaterZero(pGetCoeff(den), ntRing->cf) )
161    {
162      Print("negative sign of DEN. of a fraction in %s:%d\n",f,l);
163      return FALSE;
164    }
165
166    // test that den is over integers!?
167
168  } else
169  {  // num != NULL // den == NULL
170
171//    if( COM(t) != 0 )
172//    {
173//      Print("?//NULL with non-zero complexity: %d in %s:%d\n", COM(t), f, l);
174//      return FALSE;
175//    }
176    // test that nume is over integers!?
177  }
178  return TRUE;
179}
180#endif
181
182/* returns the bottom field in this field extension tower; if the tower
183   is flat, i.e., if there is no extension, then r itself is returned;
184   as a side-effect, the counter 'height' is filled with the height of
185   the extension tower (in case the tower is flat, 'height' is zero) */
186static coeffs nCoeff_bottom(const coeffs r, int &height)
187{
188  assume(r != NULL);
189  coeffs cf = r;
190  height = 0;
191  while (nCoeff_is_Extension(cf))
192  {
193    assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL);
194    cf = cf->extRing->cf;
195    height++;
196  }
197  return cf;
198}
199
200BOOLEAN ntIsZero(number a, const coeffs cf)
201{
202  ntTest(a); // !!!
203  return (IS0(a));
204}
205
206void ntDelete(number * a, const coeffs cf)
207{
208  ntTest(*a); // !!!
209  fraction f = (fraction)(*a);
210  if (IS0(f)) return;
211  p_Delete(&NUM(f), ntRing);
212  if (!DENIS1(f)) p_Delete(&DEN(f), ntRing);
213  omFreeBin((ADDRESS)f, fractionObjectBin);
214  *a = NULL;
215}
216
217BOOLEAN ntEqual(number a, number b, const coeffs cf)
218{
219  ntTest(a);
220  ntTest(b);
221
222  /// simple tests
223  if (a == b) return TRUE;
224  if ((IS0(a)) && (!IS0(b))) return FALSE;
225  if ((IS0(b)) && (!IS0(a))) return FALSE;
226
227  /// cheap test if gcd's have been cancelled in both numbers
228  fraction fa = (fraction)a;
229  fraction fb = (fraction)b;
230  if ((COM(fa) == 1) && (COM(fb) == 1))
231  {
232    poly f = p_Add_q(p_Copy(NUM(fa), ntRing),
233                     p_Neg(p_Copy(NUM(fb), ntRing), ntRing),
234                     ntRing);
235    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
236    if (DENIS1(fa) && DENIS1(fb))  return TRUE;
237    if (DENIS1(fa) && !DENIS1(fb)) return FALSE;
238    if (!DENIS1(fa) && DENIS1(fb)) return FALSE;
239    f = p_Add_q(p_Copy(DEN(fa), ntRing),
240                p_Neg(p_Copy(DEN(fb), ntRing), ntRing),
241                ntRing);
242    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
243    return TRUE;
244  }
245
246  /* default: the more expensive multiplication test
247              a/b = c/d  <==>  a*d = b*c */
248  poly f = p_Copy(NUM(fa), ntRing);
249  if (!DENIS1(fb)) f = p_Mult_q(f, p_Copy(DEN(fb), ntRing), ntRing);
250  poly g = p_Copy(NUM(fb), ntRing);
251  if (!DENIS1(fa)) g = p_Mult_q(g, p_Copy(DEN(fa), ntRing), ntRing);
252  poly h = p_Add_q(f, p_Neg(g, ntRing), ntRing);
253  if (h == NULL) return TRUE;
254  else
255  {
256    p_Delete(&h, ntRing);
257    return FALSE;
258  }
259}
260
261number ntCopy(number a, const coeffs cf)
262{
263  ntTest(a); // !!!
264  if (IS0(a)) return NULL;
265  fraction f = (fraction)a;
266  poly g = p_Copy(NUM(f), ntRing);
267  poly h = NULL; if (!DENIS1(f)) h = p_Copy(DEN(f), ntRing);
268  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
269  NUM(result) = g;
270  DEN(result) = h;
271  COM(result) = COM(f);
272  ntTest((number)result);
273  return (number)result;
274}
275
276/// TODO: normalization of a!?
277number ntGetNumerator(number &a, const coeffs cf)
278{
279  ntTest(a);
280  definiteGcdCancellation(a, cf, FALSE);
281
282  if (IS0(a)) return NULL;
283
284  fraction f = (fraction)a;
285  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
286
287  const BOOLEAN denis1= DENIS1 (f);
288
289  if (getCoeffType (ntCoeffs) == n_Q && !denis1)
290    handleNestedFractionsOverQ (f, cf);
291
292  if (getCoeffType (ntCoeffs) == n_Q && denis1)
293  {
294    assume( DEN (f) == NULL );
295
296    number g;
297    // TODO/NOTE: the following should not be necessary (due to
298    // Hannes!) as NUM (f) should be over Z!!!
299    CPolyCoeffsEnumerator itr(NUM(f));
300
301
302    n_ClearDenominators(itr, g, ntRing->cf);
303
304    if( !n_GreaterZero(g, ntRing->cf) )
305    {
306      NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
307      g = n_Neg(g, ntRing->cf);
308    }
309
310    // g should be a positive integer now!
311    assume( n_GreaterZero(g, ntRing->cf) );
312
313    if( !n_IsOne(g, ntRing->cf) )
314    {
315      DEN (f) = p_NSet(g, ntRing); // update COM(f)???
316      COM (f) ++;
317      assume( DEN (f) != NULL );
318    }
319    else
320      n_Delete(&g, ntRing->cf);
321
322    ntTest(a);
323  }
324
325  // Call ntNormalize instead of above?!?
326
327  NUM (result) = p_Copy (NUM (f), ntRing); // ???
328  DEN (result) = NULL;
329  COM (result) = 0;
330
331  ntTest((number)result);
332  return (number)result;
333}
334
335/// TODO: normalization of a!?
336number ntGetDenom(number &a, const coeffs cf)
337{
338  ntTest(a);
339  definiteGcdCancellation(a, cf, FALSE);
340  fraction f = (fraction)a;
341
342  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
343  DEN (result)= NULL;
344  COM (result)= 0;
345
346
347  const BOOLEAN denis1 = DENIS1 (f);
348
349  if( IS0(f) || (denis1 && getCoeffType (ntCoeffs) != n_Q) ) // */1 or 0
350  {
351    NUM (result)= p_One(ntRing);
352    ntTest((number)result);
353    return (number)result;
354  }
355
356  if (!denis1) // */* / Q
357  {
358    assume( DEN (f) != NULL );
359
360    if (getCoeffType (ntCoeffs) == n_Q)
361      handleNestedFractionsOverQ (f, cf);
362
363    ntTest(a);
364
365    if( DEN (f) != NULL ) // is it ?? // 1 now???
366    {
367      assume( !p_IsOne(DEN (f), ntRing) );
368
369      NUM (result) = p_Copy (DEN (f), ntRing);
370      ntTest((number)result);
371      return (number)result;
372    }
373//    NUM (result) = p_One(ntRing); // NOTE: just in order to be sure...
374  }
375
376  // */1 / Q
377  assume( getCoeffType (ntCoeffs) == n_Q );
378  assume( DEN (f) == NULL );
379
380  number g;
381//    poly num= p_Copy (NUM (f), ntRing); // ???
382
383
384  // TODO/NOTE: the following should not be necessary (due to
385  // Hannes!) as NUM (f) should be over Z!!!
386  CPolyCoeffsEnumerator itr(NUM(f));
387
388  n_ClearDenominators(itr, g, ntRing->cf); // may return -1 :(((
389
390  if( !n_GreaterZero(g, ntRing->cf) )
391  {
392//     NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
393//     g = n_Neg(g, ntRing->cf);
394    NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
395    g = n_Neg(g, ntRing->cf);
396  }
397
398  // g should be a positive integer now!
399  assume( n_GreaterZero(g, ntRing->cf) );
400
401  if( !n_IsOne(g, ntRing->cf) )
402  {
403    assume( n_GreaterZero(g, ntRing->cf) );
404    assume( !n_IsOne(g, ntRing->cf) );
405
406    DEN (f) = p_NSet(g, ntRing); // update COM(f)???
407    assume( DEN (f) != NULL );
408    COM (f) ++;
409
410    NUM (result)= p_Copy (DEN (f), ntRing);
411  }
412  else
413  { // common denom == 1?
414    NUM (result)= p_NSet(g, ntRing); // p_Copy (DEN (f), ntRing);
415//  n_Delete(&g, ntRing->cf);
416  }
417
418//    if (!p_IsConstant (num, ntRing) && pNext(num) != NULL)
419//    else
420//      g= p_GetAllDenom (num, ntRing);
421//    result= (fraction) ntSetMap (ntRing->cf, cf) (g, ntRing->cf, cf);
422
423  ntTest((number)result);
424  return (number)result;
425}
426
427BOOLEAN ntIsOne(number a, const coeffs cf)
428{
429  ntTest(a); // !!!
430  definiteGcdCancellation(a, cf, FALSE);
431  fraction f = (fraction)a;
432  return (f!=NULL) && DENIS1(f) && NUMIS1(f);
433}
434
435BOOLEAN ntIsMOne(number a, const coeffs cf)
436{
437  ntTest(a);
438  definiteGcdCancellation(a, cf, FALSE);
439  fraction f = (fraction)a;
440  if ((f==NULL) || (!DENIS1(f))) return FALSE;
441  poly g = NUM(f);
442  if (!p_IsConstant(g, ntRing)) return FALSE;
443  return n_IsMOne(p_GetCoeff(g, ntRing), ntCoeffs);
444}
445
446/// this is in-place, modifies a
447number ntNeg(number a, const coeffs cf)
448{
449  ntTest(a);
450  if (!IS0(a))
451  {
452    fraction f = (fraction)a;
453    NUM(f) = p_Neg(NUM(f), ntRing);
454  }
455  ntTest(a);
456  return a;
457}
458
459number ntImPart(number a, const coeffs cf)
460{
461  ntTest(a);
462  return NULL;
463}
464
465number ntInit_bigint(number longratBigIntNumber, const coeffs src, const coeffs cf)
466{
467  assume( cf != NULL );
468
469  const ring A = cf->extRing;
470
471  assume( A != NULL );
472
473  const coeffs C = A->cf;
474
475  assume( C != NULL );
476
477  number n = n_Init_bigint(longratBigIntNumber, src, C);
478
479  if ( n_IsZero(n, C) )
480  {
481    n_Delete(&n, C);
482    return NULL;
483  }
484
485  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
486
487  number den = n_GetDenom(n, C);
488
489  assume( n_GreaterZero(den, C) );
490
491  if( n_IsOne(den, C) )
492  {
493     NUM(result) = p_NSet(n, A);
494     DEN(result) = NULL;
495     n_Delete(&den, C);
496  } else
497  {
498     DEN(result) = p_NSet(den, A);
499     NUM(result) = p_NSet(n_GetNumerator(n, C), A);
500     n_Delete(&n, C);
501  }
502
503  COM(result) = 0;
504
505  ntTest((number)result);
506
507  return (number)result;
508}
509
510
511number ntInit(long i, const coeffs cf)
512{
513  if (i != 0)
514  {
515    poly p=p_ISet(i, ntRing);
516    if (p!=NULL)
517    {
518      fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
519      NUM(result) = p;
520      //DEN(result) = NULL; // done by omAlloc0Bin
521      //COM(result) = 0; // done by omAlloc0Bin
522      ntTest((number)result);
523      return (number)result;
524    }
525  }
526  return NULL;
527}
528
529
530/// takes over p!
531number ntInit(poly p, const coeffs cf)
532{
533  if (p == NULL) return NULL;
534
535    number g;
536    // TODO/NOTE: the following should not be necessary (due to
537    // Hannes!) as NUM (f) should be over Z!!!
538    CPolyCoeffsEnumerator itr(p);
539
540    n_ClearDenominators(itr, g, ntRing->cf);
541
542    if( !n_GreaterZero(g, ntRing->cf) )
543    {
544      p = p_Neg(p, ntRing); // Ugly :(((
545      g = n_Neg(g, ntRing->cf);
546    }
547
548    // g should be a positive integer now!
549    assume( n_GreaterZero(g, ntRing->cf) );
550
551    fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
552
553    if( !n_IsOne(g, ntRing->cf) )
554    {
555      DEN (f) = p_NSet(g, ntRing);
556//      COM (f) ++; // update COM(f)???
557      assume( DEN (f) != NULL );
558    }
559    else
560    {
561      DEN(f) = NULL;
562      n_Delete(&g, ntRing->cf);
563    }
564
565    NUM(f) = p;
566    COM(f) = 0;
567
568    ntTest((number)f);
569    return (number)f;
570}
571
572int ntInt(number &a, const coeffs cf)
573{
574  ntTest(a);
575  if (IS0(a)) return 0;
576  definiteGcdCancellation(a, cf, FALSE);
577  fraction f = (fraction)a;
578  if (!DENIS1(f)) return 0;
579
580  const poly aAsPoly = NUM(f);
581
582  if(aAsPoly == NULL)
583    return 0;
584
585  if (!p_IsConstant(aAsPoly, ntRing))
586    return 0;
587
588  assume( aAsPoly != NULL );
589
590  return n_Int(p_GetCoeff(aAsPoly, ntRing), ntCoeffs);
591}
592
593/* This method will only consider the numerators of a and b, without
594   cancelling gcd's before.
595   Moreover it may return TRUE only if one or both numerators
596   are zero or if their degrees are equal. Then TRUE is returned iff
597   coeff(numerator(a)) > coeff(numerator(b));
598   In all other cases, FALSE will be returned. */
599BOOLEAN ntGreater(number a, number b, const coeffs cf)
600{
601  ntTest(a);
602  ntTest(b);
603  number aNumCoeff = NULL; int aNumDeg = -1;
604  number aDenCoeff = NULL; int aDenDeg = -1;
605  number bNumCoeff = NULL; int bNumDeg = -1;
606  number bDenCoeff = NULL; int bDenDeg = -1;
607  if (!IS0(a))
608  {
609    fraction fa = (fraction)a;
610    aNumDeg = p_Totaldegree(NUM(fa), ntRing);
611    aNumCoeff = p_GetCoeff(NUM(fa), ntRing);
612    if (DEN(fa)!=NULL)
613    {
614      aDenDeg = p_Totaldegree(DEN(fa), ntRing);
615      aDenCoeff=p_GetCoeff(DEN(fa),ntRing);
616    }
617  }
618  else return !(ntGreaterZero (b,cf));
619  if (!IS0(b))
620  {
621    fraction fb = (fraction)b;
622    bNumDeg = p_Totaldegree(NUM(fb), ntRing);
623    bNumCoeff = p_GetCoeff(NUM(fb), ntRing);
624    if (DEN(fb)!=NULL)
625    {
626      bDenDeg = p_Totaldegree(DEN(fb), ntRing);
627      bDenCoeff=p_GetCoeff(DEN(fb),ntRing);
628    }
629  }
630  else return ntGreaterZero(a,cf);
631  if (aNumDeg-aDenDeg > bNumDeg-bDenDeg) return TRUE;
632  if (aNumDeg-aDenDeg < bNumDeg-bDenDeg) return FALSE;
633  number aa;
634  number bb;
635  if (bDenCoeff==NULL) aa=n_Copy(aNumCoeff,ntRing->cf);
636  else                 aa=n_Mult(aNumCoeff,bDenCoeff,ntRing->cf);
637  if (aDenCoeff==NULL) bb=n_Copy(bNumCoeff,ntRing->cf);
638  else                 bb=n_Mult(bNumCoeff,aDenCoeff,ntRing->cf);
639  BOOLEAN rr= n_Greater(aa, bb, ntCoeffs);
640  n_Delete(&aa,ntRing->cf);
641  n_Delete(&bb,ntRing->cf);
642  return rr;
643}
644
645/* this method will only consider the numerator of a, without cancelling
646   the gcd before;
647   returns TRUE iff the leading coefficient of the numerator of a is > 0
648                    or the leading term of the numerator of a is not a
649                    constant */
650BOOLEAN ntGreaterZero(number a, const coeffs cf)
651{
652  ntTest(a);
653  if (IS0(a)) return FALSE;
654  fraction f = (fraction)a;
655  poly g = NUM(f);
656  return (n_GreaterZero(p_GetCoeff(g, ntRing), ntCoeffs) ||
657          (!p_LmIsConstant(g, ntRing)));
658}
659
660void ntCoeffWrite(const coeffs cf, BOOLEAN details)
661{
662  assume( cf != NULL );
663
664  const ring A = cf->extRing;
665
666  assume( A != NULL );
667  assume( A->cf != NULL );
668
669  n_CoeffWrite(A->cf, details);
670
671//  rWrite(A);
672
673  const int P = rVar(A);
674  assume( P > 0 );
675
676  Print("//   %d parameter    : ", P);
677
678  for (int nop=0; nop < P; nop ++)
679    Print("%s ", rRingVar(nop, A));
680
681  assume( A->qideal == NULL );
682
683  PrintS("\n//   minpoly        : 0\n");
684
685/*
686  PrintS("//   Coefficients live in the rational function field\n");
687  Print("//   K(");
688  for (int i = 0; i < rVar(ntRing); i++)
689  {
690    if (i > 0) PrintS(" ");
691    Print("%s", rRingVar(i, ntRing));
692  }
693  PrintS(") with\n");
694  PrintS("//   K: "); n_CoeffWrite(cf->extRing->cf);
695*/
696}
697
698number ntAdd(number a, number b, const coeffs cf)
699{
700  ntTest(a);
701  ntTest(b);
702  if (IS0(a)) return ntCopy(b, cf);
703  if (IS0(b)) return ntCopy(a, cf);
704
705  fraction fa = (fraction)a;
706  fraction fb = (fraction)b;
707
708  poly g = p_Copy(NUM(fa), ntRing);
709  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
710  poly h = p_Copy(NUM(fb), ntRing);
711  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
712  g = p_Add_q(g, h, ntRing);
713
714  if (g == NULL) return NULL;
715
716  poly f;
717  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
718  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
719  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
720  else /* both denom's are != 1 */    f = p_Mult_q(p_Copy(DEN(fa), ntRing),
721                                                   p_Copy(DEN(fb), ntRing),
722                                                   ntRing);
723
724  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
725  NUM(result) = g;
726  DEN(result) = f;
727  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
728  heuristicGcdCancellation((number)result, cf);
729
730//  ntTest((number)result);
731
732  return (number)result;
733}
734
735number ntSub(number a, number b, const coeffs cf)
736{
737  ntTest(a);
738  ntTest(b);
739  if (IS0(a)) return ntNeg(ntCopy(b, cf), cf);
740  if (IS0(b)) return ntCopy(a, cf);
741
742  fraction fa = (fraction)a;
743  fraction fb = (fraction)b;
744
745  poly g = p_Copy(NUM(fa), ntRing);
746  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
747  poly h = p_Copy(NUM(fb), ntRing);
748  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
749  g = p_Add_q(g, p_Neg(h, ntRing), ntRing);
750
751  if (g == NULL) return NULL;
752
753  poly f;
754  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
755  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
756  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
757  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
758                                                   p_Copy(DEN(fb), ntRing),
759                                                   ntRing);
760
761  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
762  NUM(result) = g;
763  DEN(result) = f;
764  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
765  heuristicGcdCancellation((number)result, cf);
766//  ntTest((number)result);
767  return (number)result;
768}
769
770number ntMult(number a, number b, const coeffs cf)
771{
772  ntTest(a); // !!!?
773  ntTest(b); // !!!?
774
775  if (IS0(a) || IS0(b)) return NULL;
776
777  fraction fa = (fraction)a;
778  fraction fb = (fraction)b;
779
780  const poly g = pp_Mult_qq(NUM(fa), NUM(fb), ntRing);
781
782  if (g == NULL) return NULL; // may happen due to zero divisors???
783
784  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
785
786  NUM(result) = g;
787
788  const poly da = DEN(fa);
789  const poly db = DEN(fb);
790
791
792  if (db == NULL)
793  {
794    // b = ? // NULL
795
796    if(da == NULL)
797    { // both fa && fb are ?? // NULL!
798      assume (da == NULL && db == NULL);
799      DEN(result) = NULL;
800      COM(result) = 0;
801    }
802    else
803    {
804      assume (da != NULL && db == NULL);
805      DEN(result) = p_Copy(da, ntRing);
806      COM(result) = COM(fa) + MULT_COMPLEXITY;
807      heuristicGcdCancellation((number)result, cf);
808    }
809  } else
810  { // b = ?? / ??
811    if (da == NULL)
812    { // a == ? // NULL
813      assume( db != NULL && da == NULL);
814      DEN(result) = p_Copy(db, ntRing);
815      COM(result) = COM(fb) + MULT_COMPLEXITY;
816      heuristicGcdCancellation((number)result, cf);
817    }
818    else /* both den's are != 1 */
819    {
820      assume (da != NULL && db != NULL);
821      DEN(result) = pp_Mult_qq(da, db, ntRing);
822      COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
823      heuristicGcdCancellation((number)result, cf);
824    }
825  }
826
827//  ntTest((number)result);
828
829  return (number)result;
830}
831
832number ntDiv(number a, number b, const coeffs cf)
833{
834  ntTest(a);
835  ntTest(b);
836  if (IS0(a)) return NULL;
837  if (IS0(b)) WerrorS(nDivBy0);
838
839  fraction fa = (fraction)a;
840  fraction fb = (fraction)b;
841
842  poly g = p_Copy(NUM(fa), ntRing);
843  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
844
845  if (g == NULL) return NULL;   /* may happen due to zero divisors */
846
847  poly f = p_Copy(NUM(fb), ntRing);
848  if (!DENIS1(fa)) f = p_Mult_q(f, p_Copy(DEN(fa), ntRing), ntRing);
849
850  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
851  NUM(result) = g;
852  if (!n_GreaterZero(pGetCoeff(f),ntRing->cf))
853  {
854    g=p_Neg(g,ntRing);
855    f=p_Neg(f,ntRing);
856    NUM(result) = g;
857  }
858  if (!p_IsConstant(f,ntRing) || !n_IsOne(pGetCoeff(f),ntRing->cf))
859  {
860    DEN(result) = f;
861  }
862  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
863  heuristicGcdCancellation((number)result, cf);
864//  ntTest((number)result);
865  return (number)result;
866}
867
868/* 0^0 = 0;
869   for |exp| <= 7 compute power by a simple multiplication loop;
870   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
871      p^13 = p^1 * p^4 * p^8, where we utilise that
872      p^(2^(k+1)) = p^(2^k) * p^(2^k);
873   intermediate cancellation is controlled by the in-place method
874   heuristicGcdCancellation; see there.
875*/
876void ntPower(number a, int exp, number *b, const coeffs cf)
877{
878  ntTest(a);
879
880  /* special cases first */
881  if (IS0(a))
882  {
883    if (exp >= 0) *b = NULL;
884    else          WerrorS(nDivBy0);
885  }
886  else if (exp ==  0) { *b = ntInit(1, cf); return;}
887  else if (exp ==  1) { *b = ntCopy(a, cf); return;}
888  else if (exp == -1) { *b = ntInvers(a, cf); return;}
889
890  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
891
892  /* now compute a^expAbs */
893  number pow; number t;
894  if (expAbs <= 7)
895  {
896    pow = ntCopy(a, cf);
897    for (int i = 2; i <= expAbs; i++)
898    {
899      t = ntMult(pow, a, cf);
900      ntDelete(&pow, cf);
901      pow = t;
902      heuristicGcdCancellation(pow, cf);
903    }
904  }
905  else
906  {
907    pow = ntInit(1, cf);
908    number factor = ntCopy(a, cf);
909    while (expAbs != 0)
910    {
911      if (expAbs & 1)
912      {
913        t = ntMult(pow, factor, cf);
914        ntDelete(&pow, cf);
915        pow = t;
916        heuristicGcdCancellation(pow, cf);
917      }
918      expAbs = expAbs / 2;
919      if (expAbs != 0)
920      {
921        t = ntMult(factor, factor, cf);
922        ntDelete(&factor, cf);
923        factor = t;
924        heuristicGcdCancellation(factor, cf);
925      }
926    }
927    ntDelete(&factor, cf);
928  }
929
930  /* invert if original exponent was negative */
931  if (exp < 0)
932  {
933    t = ntInvers(pow, cf);
934    ntDelete(&pow, cf);
935    pow = t;
936  }
937  *b = pow;
938  ntTest(*b);
939}
940
941/* assumes that cf represents the rationals, i.e. Q, and will only
942   be called in that case;
943   assumes furthermore that f != NULL and that the denominator of f != 1;
944   generally speaking, this method removes denominators in the rational
945   coefficients of the numerator and denominator of 'a';
946   more concretely, the following normalizations will be performed,
947   where t^alpha denotes a monomial in the transcendental variables t_k
948   (1) if 'a' is of the form
949          (sum_alpha a_alpha/b_alpha * t^alpha)
950          -------------------------------------
951            (sum_beta c_beta/d_beta * t^beta)
952       with integers a_alpha, b_alpha, c_beta, d_beta, then both the
953       numerator and the denominator will be multiplied by the LCM of
954       the b_alpha's and the d_beta's (if this LCM is != 1),
955   (2) if 'a' is - e.g. after having performed step (1) - of the form
956          (sum_alpha a_alpha * t^alpha)
957          -----------------------------
958            (sum_beta c_beta * t^beta)
959       with integers a_alpha, c_beta, and with a non-constant denominator,
960       then both the numerator and the denominator will be divided by the
961       GCD of the a_alpha's and the c_beta's (if this GCD is != 1),
962   (3) if 'a' is - e.g. after having performed steps (1) and (2) - of the
963       form
964          (sum_alpha a_alpha * t^alpha)
965          -----------------------------
966                        c
967       with integers a_alpha, and c != 1, then 'a' will be replaced by
968       (sum_alpha a_alpha/c * t^alpha);
969   this procedure does not alter COM(f) (this has to be done by the
970   calling procedure);
971   modifies f */
972void handleNestedFractionsOverQ(fraction f, const coeffs cf)
973{
974  assume(nCoeff_is_Q(ntCoeffs));
975  assume(!IS0(f));
976  assume(!DENIS1(f));
977
978  if (!p_IsConstant(DEN(f), ntRing))
979  { /* step (1); see documentation of this procedure above */
980    p_Normalize(NUM(f), ntRing);
981    p_Normalize(DEN(f), ntRing);
982    number lcmOfDenominators = n_Init(1, ntCoeffs);
983    number c; number tmp;
984    poly p = NUM(f);
985    /* careful when using n_Lcm!!! It computes the lcm of the numerator
986       of the 1st argument and the denominator of the 2nd!!! */
987    while (p != NULL)
988    {
989      c = p_GetCoeff(p, ntRing);
990      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
991      n_Delete(&lcmOfDenominators, ntCoeffs);
992      lcmOfDenominators = tmp;
993      pIter(p);
994    }
995    p = DEN(f);
996    while (p != NULL)
997    {
998      c = p_GetCoeff(p, ntRing);
999      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
1000      n_Delete(&lcmOfDenominators, ntCoeffs);
1001      lcmOfDenominators = tmp;
1002      pIter(p);
1003    }
1004    if (!n_IsOne(lcmOfDenominators, ntCoeffs))
1005    { /* multiply NUM(f) and DEN(f) with lcmOfDenominators */
1006      NUM(f) = p_Mult_nn(NUM(f), lcmOfDenominators, ntRing);
1007      p_Normalize(NUM(f), ntRing);
1008      DEN(f) = p_Mult_nn(DEN(f), lcmOfDenominators, ntRing);
1009      p_Normalize(DEN(f), ntRing);
1010    }
1011    n_Delete(&lcmOfDenominators, ntCoeffs);
1012    if (!p_IsConstant(DEN(f), ntRing))
1013    { /* step (2); see documentation of this procedure above */
1014      p = NUM(f);
1015      number gcdOfCoefficients = n_Copy(p_GetCoeff(p, ntRing), ntCoeffs);
1016      pIter(p);
1017      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
1018      {
1019        c = p_GetCoeff(p, ntRing);
1020        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
1021        n_Delete(&gcdOfCoefficients, ntCoeffs);
1022        gcdOfCoefficients = tmp;
1023        pIter(p);
1024      }
1025      p = DEN(f);
1026      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
1027      {
1028        c = p_GetCoeff(p, ntRing);
1029        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
1030        n_Delete(&gcdOfCoefficients, ntCoeffs);
1031        gcdOfCoefficients = tmp;
1032        pIter(p);
1033      }
1034      if (!n_IsOne(gcdOfCoefficients, ntCoeffs))
1035      { /* divide NUM(f) and DEN(f) by gcdOfCoefficients */
1036        number inverseOfGcdOfCoefficients = n_Invers(gcdOfCoefficients,
1037                                                     ntCoeffs);
1038        NUM(f) = p_Mult_nn(NUM(f), inverseOfGcdOfCoefficients, ntRing);
1039        p_Normalize(NUM(f), ntRing);
1040        DEN(f) = p_Mult_nn(DEN(f), inverseOfGcdOfCoefficients, ntRing);
1041        p_Normalize(DEN(f), ntRing);
1042        n_Delete(&inverseOfGcdOfCoefficients, ntCoeffs);
1043      }
1044      n_Delete(&gcdOfCoefficients, ntCoeffs);
1045    }
1046  }
1047  if (p_IsConstant(DEN(f), ntRing) &&
1048      (!n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs)))
1049  { /* step (3); see documentation of this procedure above */
1050    number inverseOfDen = n_Invers(p_GetCoeff(DEN(f), ntRing), ntCoeffs);
1051    NUM(f) = p_Mult_nn(NUM(f), inverseOfDen, ntRing);
1052    n_Delete(&inverseOfDen, ntCoeffs);
1053    p_Delete(&DEN(f), ntRing);
1054    DEN(f) = NULL;
1055  }
1056
1057  /* Now, due to the above computations, DEN(f) may have become the
1058     1-polynomial which needs to be represented by NULL: */
1059  if ((DEN(f) != NULL) &&
1060      p_IsConstant(DEN(f), ntRing) &&
1061      n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
1062  {
1063    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1064  }
1065
1066  if( DEN(f) != NULL )
1067    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1068    {
1069      NUM(f) = p_Neg(NUM(f), ntRing);
1070      DEN(f) = p_Neg(DEN(f), ntRing);
1071    }
1072
1073  ntTest((number)f); // TODO!
1074}
1075
1076/* modifies a */
1077void heuristicGcdCancellation(number a, const coeffs cf)
1078{
1079//  ntTest(a); // !!!!????
1080  if (IS0(a)) return;
1081
1082  fraction f = (fraction)a;
1083  if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; ntTest(a); return; }
1084
1085  assume( DEN(f) != NULL );
1086
1087  /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
1088  if (p_EqualPolys(NUM(f), DEN(f), ntRing))
1089  { /* numerator and denominator are both != 1 */
1090    p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
1091    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1092    COM(f) = 0;
1093  }
1094  else
1095  {
1096    if (COM(f) > BOUND_COMPLEXITY)
1097      definiteGcdCancellation(a, cf, TRUE);
1098
1099  // TODO: check if it is enough to put the following into definiteGcdCancellation?!
1100  if( DEN(f) != NULL )
1101    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1102    {
1103      NUM(f) = p_Neg(NUM(f), ntRing);
1104      DEN(f) = p_Neg(DEN(f), ntRing);
1105    }
1106  }
1107
1108  ntTest(a);
1109}
1110
1111/// modifies a
1112void definiteGcdCancellation(number a, const coeffs cf,
1113                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed)
1114{
1115  ntTest(a); // !!!!
1116
1117  fraction f = (fraction)a;
1118
1119  if (IS0(a)) return;
1120  if (NUM(f)!=NULL) p_Normalize(NUM(f), ntRing);
1121  if (DEN(f)!=NULL) p_Normalize(DEN(f), ntRing);
1122  if (!simpleTestsHaveAlreadyBeenPerformed)
1123  {
1124    if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
1125
1126    /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
1127    if (p_EqualPolys(NUM(f), DEN(f), ntRing))
1128    { /* numerator and denominator are both != 1 */
1129      p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
1130      p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1131      COM(f) = 0;
1132      ntTest(a); // !!!!
1133      return;
1134    }
1135  }
1136
1137#ifdef HAVE_FACTORY
1138  /* Note that, over Q, singclap_gcd will remove the denominators in all
1139     rational coefficients of pNum and pDen, before starting to compute
1140     the gcd. Thus, we do not need to ensure that the coefficients of
1141     pNum and pDen live in Z; they may well be elements of Q\Z. */
1142  /* singclap_gcd destroys its arguments; we hence need copies: */
1143  poly pGcd = singclap_gcd(p_Copy(NUM(f), ntRing), p_Copy(DEN(f), ntRing), cf->extRing);
1144  if (p_IsConstant(pGcd, ntRing) &&
1145      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
1146  { /* gcd = 1; nothing to cancel;
1147       Suppose the given rational function field is over Q. Although the
1148       gcd is 1, we may have produced fractional coefficients in NUM(f),
1149       DEN(f), or both, due to previous arithmetics. The next call will
1150       remove those nested fractions, in case there are any. */
1151    if (nCoeff_is_Zp(ntCoeffs) && p_IsConstant (DEN (f), ntRing))
1152    {
1153      NUM (f) = p_Div_nn (NUM (f), p_GetCoeff (DEN(f),ntRing), ntRing);
1154      //poly newNum= singclap_pdivide (NUM(f), DEN (f), ntRing);
1155      //p_Delete(&NUM (f), ntRing);
1156      //NUM (f)= newNum;
1157      p_Delete(&DEN (f), ntRing);
1158      DEN (f) = NULL;
1159      COM (f) = 0;
1160    } else if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
1161  }
1162  else
1163  { /* We divide both NUM(f) and DEN(f) by the gcd which is known
1164       to be != 1. */
1165    poly newNum = singclap_pdivide(NUM(f), pGcd, ntRing);
1166    p_Delete(&NUM(f), ntRing);
1167    NUM(f) = newNum;
1168    poly newDen = singclap_pdivide(DEN(f), pGcd, ntRing);
1169    p_Delete(&DEN(f), ntRing);
1170    DEN(f) = newDen;
1171    if (p_IsConstant(DEN(f), ntRing) &&
1172        n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
1173    {
1174      /* DEN(f) = 1 needs to be represented by NULL! */
1175      p_Delete(&DEN(f), ntRing);
1176      newDen = NULL;
1177    }
1178    else
1179    { /* Note that over Q, by cancelling the gcd, we may have produced
1180         fractional coefficients in NUM(f), DEN(f), or both. The next
1181         call will remove those nested fractions, in case there are
1182         any. */
1183      if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
1184    }
1185  }
1186  COM(f) = 0;
1187  p_Delete(&pGcd, ntRing);
1188
1189  if( DEN(f) != NULL )
1190    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1191    {
1192      NUM(f) = p_Neg(NUM(f), ntRing);
1193      DEN(f) = p_Neg(DEN(f), ntRing);
1194    }
1195#endif /* HAVE_FACTORY */
1196
1197  ntTest(a); // !!!!
1198}
1199
1200// NOTE: modifies a
1201void ntWriteLong(number &a, const coeffs cf)
1202{
1203  ntTest(a);
1204  definiteGcdCancellation(a, cf, FALSE);
1205  if (IS0(a))
1206    StringAppendS("0");
1207  else
1208  {
1209    fraction f = (fraction)a;
1210    // stole logic from napWrite from kernel/longtrans.cc of legacy singular
1211    BOOLEAN omitBrackets = p_IsConstant(NUM(f), ntRing);
1212    if (!omitBrackets) StringAppendS("(");
1213    p_String0Long(NUM(f), ntRing, ntRing);
1214    if (!omitBrackets) StringAppendS(")");
1215    if (!DENIS1(f))
1216    {
1217      StringAppendS("/");
1218      omitBrackets = p_IsConstant(DEN(f), ntRing);
1219      if (!omitBrackets) StringAppendS("(");
1220      p_String0Long(DEN(f), ntRing, ntRing);
1221      if (!omitBrackets) StringAppendS(")");
1222    }
1223  }
1224  ntTest(a); // !!!!
1225}
1226
1227// NOTE: modifies a
1228void ntWriteShort(number &a, const coeffs cf)
1229{
1230  ntTest(a);
1231  definiteGcdCancellation(a, cf, FALSE);
1232  if (IS0(a))
1233    StringAppendS("0");
1234  else
1235  {
1236    fraction f = (fraction)a;
1237    // stole logic from napWrite from kernel/longtrans.cc of legacy singular
1238    BOOLEAN omitBrackets = p_IsConstant(NUM(f), ntRing);
1239    if (!omitBrackets) StringAppendS("(");
1240    p_String0Short(NUM(f), ntRing, ntRing);
1241    if (!omitBrackets) StringAppendS(")");
1242    if (!DENIS1(f))
1243    {
1244      StringAppendS("/");
1245      omitBrackets = p_IsConstant(DEN(f), ntRing);
1246      if (!omitBrackets) StringAppendS("(");
1247      p_String0Short(DEN(f), ntRing, ntRing);
1248      if (!omitBrackets) StringAppendS(")");
1249    }
1250  }
1251  ntTest(a);
1252}
1253
1254const char * ntRead(const char *s, number *a, const coeffs cf)
1255{
1256  poly p;
1257  const char * result = p_Read(s, p, ntRing);
1258  if (p == NULL) *a = NULL;
1259  else *a = ntInit(p, cf);
1260  return result;
1261}
1262
1263void ntNormalize (number &a, const coeffs cf)
1264{
1265  if ((a!=NULL))
1266  {
1267    definiteGcdCancellation(a, cf, FALSE);
1268    if ((DEN(a)!=NULL)
1269    &&(!n_GreaterZero(pGetCoeff(DEN(a)),ntRing->cf)))
1270    {
1271      NUM(a)=p_Neg(NUM(a),ntRing);
1272      DEN(a)=p_Neg(DEN(a),ntRing);
1273    }
1274  }
1275  ntTest(a); // !!!!
1276}
1277
1278/* expects *param to be castable to TransExtInfo */
1279static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
1280{
1281  if (ID != n) return FALSE;
1282  TransExtInfo *e = (TransExtInfo *)param;
1283  /* for rational function fields we expect the underlying
1284     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
1285     this expectation is based on the assumption that we have properly
1286     registered cf and perform reference counting rather than creating
1287     multiple copies of the same coefficient field/domain/ring */
1288  if (ntRing == e->r)
1289    return TRUE;
1290
1291  // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)...
1292  if( rEqual(ntRing, e->r, TRUE) )
1293  {
1294    rDelete(e->r);
1295    return TRUE;
1296  }
1297
1298  return FALSE;
1299}
1300
1301number ntLcm(number a, number b, const coeffs cf)
1302{
1303  ntTest(a);
1304  ntTest(b);
1305  fraction fb = (fraction)b;
1306  if ((b==NULL)||(DEN(fb)==NULL)) return ntCopy(a,cf);
1307#ifdef HAVE_FACTORY
1308  fraction fa = (fraction)a;
1309  /* singclap_gcd destroys its arguments; we hence need copies: */
1310  poly pa = p_Copy(NUM(fa), ntRing);
1311  poly pb = p_Copy(DEN(fb), ntRing);
1312
1313  /* Note that, over Q, singclap_gcd will remove the denominators in all
1314     rational coefficients of pa and pb, before starting to compute
1315     the gcd. Thus, we do not need to ensure that the coefficients of
1316     pa and pb live in Z; they may well be elements of Q\Z. */
1317  poly pGcd = singclap_gcd(pa, pb, cf->extRing);
1318  if (p_IsConstant(pGcd, ntRing) &&
1319      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
1320  { /* gcd = 1; return pa*pb*/
1321    p_Delete(&pGcd,ntRing);
1322    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1323    NUM(result) = pp_Mult_qq(NUM(fa),DEN(fb),ntRing);
1324
1325    ntTest((number)result); // !!!!
1326
1327    return (number)result;
1328  }
1329
1330
1331  /* return pa*pb/gcd */
1332    poly newNum = singclap_pdivide(NUM(fa), pGcd, ntRing);
1333    p_Delete(&pGcd,ntRing);
1334    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1335    NUM(result) = p_Mult_q(p_Copy(DEN(fb),ntRing),newNum,ntRing);
1336    ntTest((number)result); // !!!!
1337    return (number)result;
1338
1339#else
1340  Print("// factory needed: transext.cc:ntLcm\n");
1341  return NULL;
1342#endif /* HAVE_FACTORY */
1343  return NULL;
1344}
1345
1346number ntGcd(number a, number b, const coeffs cf)
1347{
1348  ntTest(a);
1349  ntTest(b);
1350  if (a==NULL) return ntCopy(b,cf);
1351  if (b==NULL) return ntCopy(a,cf);
1352#ifdef HAVE_FACTORY
1353  fraction fa = (fraction)a;
1354  fraction fb = (fraction)b;
1355  /* singclap_gcd destroys its arguments; we hence need copies: */
1356  poly pa = p_Copy(NUM(fa), ntRing);
1357  poly pb = p_Copy(NUM(fb), ntRing);
1358
1359  /* Note that, over Q, singclap_gcd will remove the denominators in all
1360     rational coefficients of pa and pb, before starting to compute
1361     the gcd. Thus, we do not need to ensure that the coefficients of
1362     pa and pb live in Z; they may well be elements of Q\Z. */
1363  poly pGcd = singclap_gcd(pa, pb, cf->extRing);
1364  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1365  NUM(result) = pGcd;
1366  ntTest((number)result); // !!!!
1367  return (number)result;
1368#else
1369  Print("// factory needed: transext.cc:ntGcd\n");
1370  return NULL;
1371#endif /* HAVE_FACTORY */
1372}
1373
1374int ntSize(number a, const coeffs cf)
1375{
1376  ntTest(a);
1377  if (IS0(a)) return -1;
1378  /* this has been taken from the old implementation of field extensions,
1379     where we computed the sum of the degrees and the numbers of terms in
1380     the numerator and denominator of a; so we leave it at that, for the
1381     time being */
1382  fraction f = (fraction)a;
1383  poly p = NUM(f);
1384  int noOfTerms = 0;
1385  int numDegree = 0;
1386  while (p != NULL)
1387  {
1388    noOfTerms++;
1389    int d = 0;
1390    for (int i = 1; i <= rVar(ntRing); i++)
1391      d += p_GetExp(p, i, ntRing);
1392    if (d > numDegree) numDegree = d;
1393    pIter(p);
1394  }
1395  int denDegree = 0;
1396  if (!DENIS1(f))
1397  {
1398    p = DEN(f);
1399    while (p != NULL)
1400    {
1401      noOfTerms++;
1402      int d = 0;
1403      for (int i = 1; i <= rVar(ntRing); i++)
1404        d += p_GetExp(p, i, ntRing);
1405      if (d > denDegree) denDegree = d;
1406      pIter(p);
1407    }
1408  }
1409  ntTest(a); // !!!!
1410  return numDegree + denDegree + noOfTerms;
1411}
1412
1413number ntInvers(number a, const coeffs cf)
1414{
1415  ntTest(a);
1416  if (IS0(a))
1417  {
1418    WerrorS(nDivBy0);
1419    return NULL;
1420  }
1421  fraction f = (fraction)a;
1422  assume( f != NULL );
1423
1424  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1425
1426  assume( NUM(f) != NULL );
1427  const poly den = DEN(f);
1428
1429  if (den == NULL)
1430    NUM(result) = p_One(ntRing);
1431  else
1432    NUM(result) = p_Copy(den, ntRing);
1433
1434  if( !NUMIS1(f) )
1435  {
1436    poly num_f=NUM(f);
1437    BOOLEAN neg= !n_GreaterZero(pGetCoeff(num_f),ntRing->cf);
1438    if (neg)
1439    {
1440      num_f=p_Neg(p_Copy(num_f, ntRing), ntRing);
1441      NUM(result)=p_Neg(NUM(result), ntRing);
1442    }
1443    else
1444    {
1445      num_f=p_Copy(num_f, ntRing);
1446    }
1447    DEN(result) = num_f;
1448    COM(result) = COM(f);
1449    if (neg)
1450    {
1451      if (p_IsOne(num_f, ntRing))
1452      {
1453        DEN(result)=NULL;
1454        COM(result) = 0;
1455        p_Delete(&num_f,ntRing);
1456      }
1457    }
1458  }
1459  else
1460  {
1461    DEN(result) = NULL;
1462    COM(result) = 0;
1463  }
1464  ntTest((number)result); // !!!!
1465  return (number)result;
1466}
1467
1468/* assumes that src = Q, dst = Q(t_1, ..., t_s) */
1469number ntMap00(number a, const coeffs src, const coeffs dst)
1470{
1471  if (n_IsZero(a, src)) return NULL;
1472  assume(n_Test(a, src));
1473  assume(src == dst->extRing->cf);
1474  return ntInit(p_NSet(n_Copy(a, src), dst->extRing), dst);
1475}
1476
1477/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
1478number ntMapP0(number a, const coeffs src, const coeffs dst)
1479{
1480  if (n_IsZero(a, src)) return NULL;
1481  assume(n_Test(a, src));
1482  /* mapping via intermediate int: */
1483  int n = n_Int(a, src);
1484  number q = n_Init(n, dst->extRing->cf);
1485  if (n_IsZero(q, dst->extRing->cf))
1486  {
1487    n_Delete(&q, dst->extRing->cf);
1488    return NULL;
1489  }
1490  return ntInit(p_NSet(q, dst->extRing), dst);
1491}
1492
1493/* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
1494                       src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
1495number ntCopyMap(number a, const coeffs cf, const coeffs dst)
1496{
1497//  if (n_IsZero(a, cf)) return NULL;
1498
1499  ntTest(a);
1500
1501  if (IS0(a)) return NULL;
1502
1503  const ring rSrc = cf->extRing;
1504  const ring rDst = dst->extRing;
1505
1506  if( rSrc == rDst )
1507    return ntCopy(a, dst); // USUALLY WRONG!
1508
1509  fraction f = (fraction)a;
1510  poly g = prCopyR(NUM(f), rSrc, rDst);
1511
1512  poly h = NULL;
1513
1514  if (!DENIS1(f))
1515     h = prCopyR(DEN(f), rSrc, rDst);
1516
1517  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1518
1519  NUM(result) = g;
1520  DEN(result) = h;
1521  COM(result) = COM(f);
1522  assume(n_Test((number)result, dst));
1523  return (number)result;
1524}
1525
1526number ntCopyAlg(number a, const coeffs cf, const coeffs dst)
1527{
1528  assume( n_Test(a, cf) );
1529  if (n_IsZero(a, cf)) return NULL;
1530
1531  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1532  // DEN(f) = NULL; COM(f) = 0;
1533  NUM(f) = prCopyR((poly)a, cf->extRing, dst->extRing);
1534  assume(n_Test((number)f, dst));
1535  return (number)f;
1536}
1537
1538/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
1539number ntMap0P(number a, const coeffs src, const coeffs dst)
1540{
1541  assume( n_Test(a, src) );
1542  if (n_IsZero(a, src)) return NULL;
1543  // int p = rChar(dst->extRing);
1544  number q = nlModP(a, src, dst->extRing->cf);
1545
1546  if (n_IsZero(q, dst->extRing->cf))
1547  {
1548    n_Delete(&q, dst->extRing->cf);
1549    return NULL;
1550  }
1551
1552  poly g = p_NSet(q, dst->extRing);
1553
1554  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1555  NUM(f) = g; // DEN(f) = NULL; COM(f) = 0;
1556  assume(n_Test((number)f, dst));
1557  return (number)f;
1558}
1559
1560/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
1561number ntMapPP(number a, const coeffs src, const coeffs dst)
1562{
1563  assume( n_Test(a, src) );
1564  if (n_IsZero(a, src)) return NULL;
1565  assume(src == dst->extRing->cf);
1566  poly p = p_One(dst->extRing);
1567  p_SetCoeff(p, n_Copy(a, src), dst->extRing);
1568  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1569  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1570  assume(n_Test((number)f, dst));
1571  return (number)f;
1572}
1573
1574/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
1575number ntMapUP(number a, const coeffs src, const coeffs dst)
1576{
1577  assume( n_Test(a, src) );
1578  if (n_IsZero(a, src)) return NULL;
1579  /* mapping via intermediate int: */
1580  int n = n_Int(a, src);
1581  number q = n_Init(n, dst->extRing->cf);
1582  poly p;
1583  if (n_IsZero(q, dst->extRing->cf))
1584  {
1585    n_Delete(&q, dst->extRing->cf);
1586    return NULL;
1587  }
1588  p = p_One(dst->extRing);
1589  p_SetCoeff(p, q, dst->extRing);
1590  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1591  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1592  assume(n_Test((number)f, dst));
1593  return (number)f;
1594}
1595
1596nMapFunc ntSetMap(const coeffs src, const coeffs dst)
1597{
1598  /* dst is expected to be a rational function field */
1599  assume(getCoeffType(dst) == ID);
1600
1601  if( src == dst ) return ndCopyMap;
1602
1603  int h = 0; /* the height of the extension tower given by dst */
1604  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
1605  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
1606
1607  /* for the time being, we only provide maps if h = 1 and if b is Q or
1608     some field Z/pZ: */
1609  if (h==0)
1610  {
1611    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
1612      return ntMap00;                                 /// Q       -->  Q(T)
1613    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
1614      return ntMapP0;                                 /// Z/p     -->  Q(T)
1615    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
1616      return ntMap0P;                                 /// Q       --> Z/p(T)
1617    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
1618    {
1619      if (src->ch == dst->ch) return ntMapPP;         /// Z/p     --> Z/p(T)
1620      else return ntMapUP;                            /// Z/u     --> Z/p(T)
1621    }
1622  }
1623  if (h != 1) return NULL;
1624  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
1625
1626  /* Let T denote the sequence of transcendental extension variables, i.e.,
1627     K[t_1, ..., t_s] =: K[T];
1628     Let moreover, for any such sequence T, T' denote any subsequence of T
1629     of the form t_1, ..., t_w with w <= s. */
1630
1631  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
1632
1633  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
1634  {
1635    if (rVar(src->extRing) > rVar(dst->extRing))
1636       return NULL;
1637
1638    for (int i = 0; i < rVar(src->extRing); i++)
1639      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0)
1640         return NULL;
1641
1642    if (src->type==n_transExt)
1643       return ntCopyMap;          /// Q(T')   --> Q(T)
1644    else
1645       return ntCopyAlg;
1646  }
1647
1648  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
1649  {
1650    if (rVar(src->extRing) > rVar(dst->extRing))
1651       return NULL;
1652
1653    for (int i = 0; i < rVar(src->extRing); i++)
1654      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0)
1655         return NULL;
1656
1657    if (src->type==n_transExt)
1658       return ntCopyMap;         /// Z/p(T') --> Z/p(T)
1659    else
1660       return ntCopyAlg;
1661  }
1662
1663  return NULL;                                 /// default
1664}
1665#if 0
1666nMapFunc ntSetMap_T(const coeffs src, const coeffs dst)
1667{
1668  nMapFunc n=ntSetMap(src,dst);
1669  if (n==ntCopyAlg) printf("n=ntCopyAlg\n");
1670  else if (n==ntCopyMap) printf("n=ntCopyMap\n");
1671  else if (n==ntMapUP) printf("n=ntMapUP\n");
1672  else if (n==ntMap0P) printf("n=ntMap0P\n");
1673  else if (n==ntMapP0) printf("n=ntMapP0\n");
1674  else if (n==ntMap00) printf("n=ntMap00\n");
1675  else if (n==NULL) printf("n=NULL\n");
1676  else printf("n=?\n");
1677  return n;
1678}
1679#endif
1680
1681void ntKillChar(coeffs cf)
1682{
1683  if ((--cf->extRing->ref) == 0)
1684    rDelete(cf->extRing);
1685}
1686#ifdef HAVE_FACTORY
1687number ntConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
1688{
1689  if (n.isZero()) return NULL;
1690  poly p=convFactoryPSingP(n,ntRing);
1691  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1692  NUM(result) = p;
1693  //DEN(result) = NULL; // done by omAlloc0Bin
1694  //COM(result) = 0; // done by omAlloc0Bin
1695  ntTest((number)result);
1696  return (number)result;
1697}
1698CanonicalForm ntConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
1699{
1700  ntTest(n);
1701  if (IS0(n)) return CanonicalForm(0);
1702
1703  fraction f = (fraction)n;
1704  return convSingPFactoryP(NUM(f),ntRing);
1705}
1706#endif
1707
1708static int ntParDeg(number a, const coeffs cf)
1709{
1710  ntTest(a);
1711  if (IS0(a)) return -1;
1712  fraction fa = (fraction)a;
1713  return cf->extRing->pFDeg(NUM(fa),cf->extRing);
1714}
1715
1716/// return the specified parameter as a number in the given trans.ext.
1717static number ntParameter(const int iParameter, const coeffs cf)
1718{
1719  assume(getCoeffType(cf) == ID);
1720
1721  const ring R = cf->extRing;
1722  assume( R != NULL );
1723  assume( 0 < iParameter && iParameter <= rVar(R) );
1724
1725  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
1726
1727//  return (number) p;
1728
1729  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1730  NUM(f) = p;
1731  DEN(f) = NULL;
1732  COM(f) = 0;
1733
1734  ntTest((number)f);
1735
1736  return (number)f;
1737}
1738
1739/// if m == var(i)/1 => return i,
1740int ntIsParam(number m, const coeffs cf)
1741{
1742  ntTest(m);
1743  assume(getCoeffType(cf) == ID);
1744
1745  const ring R = cf->extRing;
1746  assume( R != NULL );
1747
1748  fraction f = (fraction)m;
1749
1750  if( DEN(f) != NULL )
1751    return 0;
1752
1753  return p_Var( NUM(f), R );
1754}
1755
1756struct NTNumConverter
1757{
1758  static inline poly convert(const number& n)
1759  {
1760    // suitable for trans. ext. numbers that are fractions of polys
1761    return NUM((fraction)n); // return the numerator
1762  }
1763};
1764
1765
1766static void ntClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1767{
1768  assume(cf != NULL);
1769  assume(getCoeffType(cf) == ID);
1770  // all coeffs are given by fractions of polynomails over integers!!!
1771  // without denominators!!!
1772
1773  const ring   R = cf->extRing;
1774  assume(R != NULL);
1775  const coeffs Q = R->cf;
1776  assume(Q != NULL);
1777  assume(nCoeff_is_Q(Q));
1778
1779
1780  numberCollectionEnumerator.Reset();
1781
1782  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
1783  {
1784    c = ntInit(1, cf);
1785    return;
1786  }
1787
1788  // all coeffs are given by integers after returning from this routine
1789
1790  // part 1, collect product of all denominators /gcds
1791  poly cand = NULL;
1792
1793  do
1794  {
1795    number &n = numberCollectionEnumerator.Current();
1796
1797    ntNormalize(n, cf);
1798
1799    fraction f = (fraction)n;
1800
1801    assume( f != NULL );
1802
1803    const poly den = DEN(f);
1804
1805    assume( den == NULL ); // ?? / 1 ?
1806
1807    const poly num = NUM(f);
1808
1809    if( cand == NULL )
1810      cand = p_Copy(num, R);
1811    else
1812      cand = singclap_gcd(cand, p_Copy(num, R), R); // gcd(cand, num)
1813
1814    if( p_IsConstant(cand, R) )
1815      break;
1816  }
1817  while( numberCollectionEnumerator.MoveNext() ) ;
1818
1819
1820  // part2: all coeffs = all coeffs * cand
1821  if( cand != NULL )
1822  {
1823  if( !p_IsConstant(cand, R) )
1824  {
1825    c = ntInit(cand, cf);
1826    numberCollectionEnumerator.Reset();
1827    while (numberCollectionEnumerator.MoveNext() )
1828    {
1829      number &n = numberCollectionEnumerator.Current();
1830      const number t = ntDiv(n, c, cf); // TODO: rewrite!?
1831      ntDelete(&n, cf);
1832      n = t;
1833    }
1834  } // else NUM (result) = p_One(R);
1835  else { p_Delete(&cand, R); cand = NULL; }
1836  }
1837
1838  // Quick and dirty fix for constant content clearing: consider numerators???
1839  CRecursivePolyCoeffsEnumerator<NTNumConverter> itr(numberCollectionEnumerator); // recursively treat the NUM(numbers) as polys!
1840  number cc;
1841
1842  n_ClearContent(itr, cc, Q);
1843  number g = ntInit(p_NSet(cc, R), cf);
1844
1845  if( cand != NULL )
1846  {
1847    number gg = ntMult(g, c, cf);
1848    ntDelete(&g, cf);
1849    ntDelete(&c, cf); c = gg;
1850  } else
1851    c = g;
1852  ntTest(c);
1853}
1854
1855static void ntClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1856{
1857  assume(cf != NULL);
1858  assume(getCoeffType(cf) == ID); // both over Q(a) and Zp(a)!
1859  // all coeffs are given by fractions of polynomails over integers!!!
1860
1861  numberCollectionEnumerator.Reset();
1862
1863  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
1864  {
1865    c = ntInit(1, cf);
1866    return;
1867  }
1868
1869  // all coeffs are given by integers after returning from this routine
1870
1871  // part 1, collect product of all denominators /gcds
1872  poly cand = NULL;
1873
1874  const ring R = cf->extRing;
1875  assume(R != NULL);
1876
1877  const coeffs Q = R->cf;
1878  assume(Q != NULL);
1879//  assume(nCoeff_is_Q(Q));
1880
1881  do
1882  {
1883    number &n = numberCollectionEnumerator.Current();
1884
1885    ntNormalize(n, cf);
1886
1887    fraction f = (fraction)ntGetDenom (n, cf);
1888
1889    assume( f != NULL );
1890
1891    const poly den = NUM(f);
1892
1893    if( den == NULL ) // ?? / 1 ?
1894      continue;
1895
1896    if( cand == NULL )
1897      cand = p_Copy(den, R);
1898    else
1899    {
1900      // cand === LCM( cand, den )!!!!
1901      // NOTE: maybe it's better to make the product and clearcontent afterwards!?
1902      // TODO: move the following to factory?
1903      poly gcd = singclap_gcd(p_Copy(cand, R), p_Copy(den, R), R); // gcd(cand, den) is monic no mater leading coeffs! :((((
1904      if (nCoeff_is_Q (Q))
1905      {
1906        number LcGcd= n_Gcd (p_GetCoeff (cand, R), p_GetCoeff(den, R), Q);
1907        gcd = p_Mult_nn(gcd, LcGcd, R);
1908        n_Delete(&LcGcd,Q);
1909      }
1910//      assume( n_IsOne(pGetCoeff(gcd), Q) ); // TODO: this may be wrong...
1911      cand = p_Mult_q(cand, p_Copy(den, R), R); // cand *= den
1912      const poly t = singclap_pdivide( cand, gcd, R ); // cand' * den / gcd(cand', den)
1913      p_Delete(&cand, R);
1914      p_Delete(&gcd, R);
1915      cand = t;
1916    }
1917  }
1918  while( numberCollectionEnumerator.MoveNext() );
1919
1920  if( cand == NULL )
1921  {
1922    c = ntInit(1, cf);
1923    return;
1924  }
1925
1926  c = ntInit(cand, cf);
1927
1928  numberCollectionEnumerator.Reset();
1929
1930  number d = NULL;
1931
1932  while (numberCollectionEnumerator.MoveNext() )
1933  {
1934    number &n = numberCollectionEnumerator.Current();
1935    number t = ntMult(n, c, cf); // TODO: rewrite!?
1936    ntDelete(&n, cf);
1937
1938    ntNormalize(t, cf); // TODO: needed?
1939    n = t;
1940
1941    fraction f = (fraction)t;
1942    assume( f != NULL );
1943
1944    const poly den = DEN(f);
1945
1946    if( den != NULL ) // ?? / ?? ?
1947    {
1948      assume( p_IsConstant(den, R) );
1949      assume( pNext(den) == NULL );
1950
1951      if( d == NULL )
1952        d = n_Copy(pGetCoeff(den), Q);
1953      else
1954      {
1955        number g = n_Lcm(d, pGetCoeff(den), Q);
1956        n_Delete(&d, Q); d = g;
1957      }
1958    }
1959  }
1960
1961  if( d != NULL )
1962  {
1963    numberCollectionEnumerator.Reset();
1964    while (numberCollectionEnumerator.MoveNext() )
1965    {
1966      number &n = numberCollectionEnumerator.Current();
1967      fraction f = (fraction)n;
1968
1969      assume( f != NULL );
1970
1971      const poly den = DEN(f);
1972
1973      if( den == NULL ) // ?? / 1 ?
1974        NUM(f) = p_Mult_nn(NUM(f), d, R);
1975      else
1976      {
1977        assume( p_IsConstant(den, R) );
1978        assume( pNext(den) == NULL );
1979
1980        number ddd = n_Div(d, pGetCoeff(den), Q); // but be an integer now!!!
1981        NUM(f) = p_Mult_nn(NUM(f), ddd, R);
1982        n_Delete(&ddd, Q);
1983
1984        p_Delete(&DEN(f), R);
1985        DEN(f) = NULL; // TODO: check if this is needed!?
1986      }
1987
1988      assume( DEN(f) == NULL );
1989    }
1990
1991    NUM(c) = p_Mult_nn(NUM(c), d, R);
1992    n_Delete(&d, Q);
1993  }
1994
1995
1996  ntTest(c);
1997}
1998
1999BOOLEAN ntInitChar(coeffs cf, void * infoStruct)
2000{
2001
2002  assume( infoStruct != NULL );
2003
2004  TransExtInfo *e = (TransExtInfo *)infoStruct;
2005
2006  assume( e->r                != NULL);      // extRing;
2007  assume( e->r->cf            != NULL);      // extRing->cf;
2008  assume( e->r->qideal == NULL );
2009
2010  assume( cf != NULL );
2011  assume(getCoeffType(cf) == ID);                // coeff type;
2012
2013  ring R = e->r;
2014  assume(R != NULL);
2015
2016  R->ref ++; // increase the ref.counter for the ground poly. ring!
2017
2018  cf->extRing           = R;
2019  /* propagate characteristic up so that it becomes
2020     directly accessible in cf: */
2021  cf->ch = R->cf->ch;
2022  cf->factoryVarOffset = R->cf->factoryVarOffset + rVar(R);
2023
2024  cf->cfGreaterZero  = ntGreaterZero;
2025  cf->cfGreater      = ntGreater;
2026  cf->cfEqual        = ntEqual;
2027  cf->cfIsZero       = ntIsZero;
2028  cf->cfIsOne        = ntIsOne;
2029  cf->cfIsMOne       = ntIsMOne;
2030  cf->cfInit         = ntInit;
2031  cf->cfInit_bigint  = ntInit_bigint;
2032  cf->cfInt          = ntInt;
2033  cf->cfNeg          = ntNeg;
2034  cf->cfAdd          = ntAdd;
2035  cf->cfSub          = ntSub;
2036  cf->cfMult         = ntMult;
2037  cf->cfDiv          = ntDiv;
2038  cf->cfExactDiv     = ntDiv;
2039  cf->cfPower        = ntPower;
2040  cf->cfCopy         = ntCopy;
2041  cf->cfWriteLong    = ntWriteLong;
2042  cf->cfRead         = ntRead;
2043  cf->cfNormalize    = ntNormalize;
2044  cf->cfDelete       = ntDelete;
2045  cf->cfSetMap       = ntSetMap;
2046  cf->cfGetDenom     = ntGetDenom;
2047  cf->cfGetNumerator = ntGetNumerator;
2048  cf->cfRePart       = ntCopy;
2049  cf->cfImPart       = ntImPart;
2050  cf->cfCoeffWrite   = ntCoeffWrite;
2051#ifdef LDEBUG
2052  cf->cfDBTest       = ntDBTest;
2053#endif
2054  cf->cfGcd          = ntGcd;
2055  cf->cfLcm          = ntLcm;
2056  cf->cfSize         = ntSize;
2057  cf->nCoeffIsEqual  = ntCoeffIsEqual;
2058  cf->cfInvers       = ntInvers;
2059  cf->cfIntDiv       = ntDiv;
2060  cf->cfKillChar     = ntKillChar;
2061
2062  if( rCanShortOut(ntRing) )
2063    cf->cfWriteShort = ntWriteShort;
2064  else
2065    cf->cfWriteShort = ntWriteLong;
2066
2067#ifndef HAVE_FACTORY
2068  PrintS("// Warning: The 'factory' module is not available.\n");
2069  PrintS("//          Hence gcd's cannot be cancelled in any\n");
2070  PrintS("//          computed fraction!\n");
2071#else
2072  cf->convFactoryNSingN =ntConvFactoryNSingN;
2073  cf->convSingNFactoryN =ntConvSingNFactoryN;
2074#endif
2075  cf->cfParDeg = ntParDeg;
2076
2077  cf->iNumberOfParameters = rVar(R);
2078  cf->pParameterNames = R->names;
2079  cf->cfParameter = ntParameter;
2080
2081  if( nCoeff_is_Q(R->cf) )
2082    cf->cfClearContent = ntClearContent;
2083
2084  cf->cfClearDenominators = ntClearDenominators;
2085
2086  return FALSE;
2087}
Note: See TracBrowser for help on using the repository browser.