source: git/libpolys/polys/ext_fields/transext.cc @ e5c9e5

jengelh-datetimespielwiese
Last change on this file since e5c9e5 was e5c9e5, checked in by Hans Schoenemann <hannes@…>, 10 years ago
fix: normalization for trans. ext.
  • Property mode set to 100644
File size: 57.0 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), ntCoeffs)))
155    {
156      Print("?/1 in %s:%d\n",f,l);
157      return FALSE;
158    }
159
160    if( !n_GreaterZero(pGetCoeff(den), ntCoeffs) )
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  if (IS0(a)) return NULL;
281
282  definiteGcdCancellation(a, cf, FALSE);
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, ntCoeffs);
303
304    if( !n_GreaterZero(g, ntCoeffs) )
305    {
306      NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
307      g = n_Neg(g, ntCoeffs);
308    }
309
310    // g should be a positive integer now!
311    assume( n_GreaterZero(g, ntCoeffs) );
312
313    if( !n_IsOne(g, ntCoeffs) )
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, ntCoeffs);
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, ntCoeffs); // may return -1 :(((
389
390  if( !n_GreaterZero(g, ntCoeffs) )
391  {
392//     NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
393//     g = n_Neg(g, ntCoeffs);
394    NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
395    g = n_Neg(g, ntCoeffs);
396  }
397
398  // g should be a positive integer now!
399  assume( n_GreaterZero(g, ntCoeffs) );
400
401  if( !n_IsOne(g, ntCoeffs) )
402  {
403    assume( n_GreaterZero(g, ntCoeffs) );
404    assume( !n_IsOne(g, ntCoeffs) );
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, ntCoeffs);
416  }
417
418//    if (!p_IsConstant (num, ntRing) && pNext(num) != NULL)
419//    else
420//      g= p_GetAllDenom (num, ntRing);
421//    result= (fraction) ntSetMap (ntCoeffs, cf) (g, ntCoeffs, 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  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
536
537  if (nCoeff_is_Q(ntCoeffs))
538  {
539    number g;
540    // TODO/NOTE: the following should not be necessary (due to
541    // Hannes!) as NUM (f) should be over Z!!!
542    // but it is not: normalize it
543    CPolyCoeffsEnumerator itr(p);
544
545    n_ClearDenominators(itr, g, ntCoeffs);
546
547    if( !n_GreaterZero(g, ntCoeffs) )
548    {
549      p = p_Neg(p, ntRing);
550      g = n_Neg(g, ntCoeffs);
551    }
552
553    // g should be a positive integer now!
554    assume( n_GreaterZero(g, ntCoeffs) );
555
556    if( !n_IsOne(g, ntCoeffs) )
557    {
558      DEN (f) = p_NSet(g, ntRing);
559      assume( DEN (f) != NULL );
560    }
561    else
562    {
563      //DEN(f) = NULL; // done by omAlloc0
564      n_Delete(&g, ntCoeffs);
565    }
566  }
567
568  NUM(f) = p;
569  COM(f) = 0;
570
571  ntTest((number)f);
572  return (number)f;
573}
574
575int ntInt(number &a, const coeffs cf)
576{
577  ntTest(a);
578  if (IS0(a)) return 0;
579  definiteGcdCancellation(a, cf, FALSE);
580  fraction f = (fraction)a;
581  if (!DENIS1(f)) return 0;
582
583  const poly aAsPoly = NUM(f);
584
585  if(aAsPoly == NULL)
586    return 0;
587
588  if (!p_IsConstant(aAsPoly, ntRing))
589    return 0;
590
591  assume( aAsPoly != NULL );
592
593  return n_Int(p_GetCoeff(aAsPoly, ntRing), ntCoeffs);
594}
595
596/* This method will only consider the numerators of a and b, without
597   cancelling gcd's before.
598   Moreover it may return TRUE only if one or both numerators
599   are zero or if their degrees are equal. Then TRUE is returned iff
600   coeff(numerator(a)) > coeff(numerator(b));
601   In all other cases, FALSE will be returned. */
602BOOLEAN ntGreater(number a, number b, const coeffs cf)
603{
604  ntTest(a);
605  ntTest(b);
606  number aNumCoeff = NULL; int aNumDeg = -1;
607  number aDenCoeff = NULL; int aDenDeg = -1;
608  number bNumCoeff = NULL; int bNumDeg = -1;
609  number bDenCoeff = NULL; int bDenDeg = -1;
610  if (!IS0(a))
611  {
612    fraction fa = (fraction)a;
613    aNumDeg = p_Totaldegree(NUM(fa), ntRing);
614    aNumCoeff = p_GetCoeff(NUM(fa), ntRing);
615    if (DEN(fa)!=NULL)
616    {
617      aDenDeg = p_Totaldegree(DEN(fa), ntRing);
618      aDenCoeff=p_GetCoeff(DEN(fa),ntRing);
619    }
620  }
621  else return !(ntGreaterZero (b,cf));
622  if (!IS0(b))
623  {
624    fraction fb = (fraction)b;
625    bNumDeg = p_Totaldegree(NUM(fb), ntRing);
626    bNumCoeff = p_GetCoeff(NUM(fb), ntRing);
627    if (DEN(fb)!=NULL)
628    {
629      bDenDeg = p_Totaldegree(DEN(fb), ntRing);
630      bDenCoeff=p_GetCoeff(DEN(fb),ntRing);
631    }
632  }
633  else return ntGreaterZero(a,cf);
634  if (aNumDeg-aDenDeg > bNumDeg-bDenDeg) return TRUE;
635  if (aNumDeg-aDenDeg < bNumDeg-bDenDeg) return FALSE;
636  number aa;
637  number bb;
638  if (bDenCoeff==NULL) aa=n_Copy(aNumCoeff,ntCoeffs);
639  else                 aa=n_Mult(aNumCoeff,bDenCoeff,ntCoeffs);
640  if (aDenCoeff==NULL) bb=n_Copy(bNumCoeff,ntCoeffs);
641  else                 bb=n_Mult(bNumCoeff,aDenCoeff,ntCoeffs);
642  BOOLEAN rr= n_Greater(aa, bb, ntCoeffs);
643  n_Delete(&aa,ntCoeffs);
644  n_Delete(&bb,ntCoeffs);
645  return rr;
646}
647
648/* this method will only consider the numerator of a, without cancelling
649   the gcd before;
650   returns TRUE iff the leading coefficient of the numerator of a is > 0
651                    or the leading term of the numerator of a is not a
652                    constant */
653BOOLEAN ntGreaterZero(number a, const coeffs cf)
654{
655  ntTest(a);
656  if (IS0(a)) return FALSE;
657  fraction f = (fraction)a;
658  poly g = NUM(f);
659  return (!p_LmIsConstant(g,ntRing)|| n_GreaterZero(pGetCoeff(g), ntCoeffs));
660}
661
662void ntCoeffWrite(const coeffs cf, BOOLEAN details)
663{
664  assume( cf != NULL );
665
666  const ring A = cf->extRing;
667
668  assume( A != NULL );
669  assume( A->cf != NULL );
670
671  n_CoeffWrite(A->cf, details);
672
673//  rWrite(A);
674
675  const int P = rVar(A);
676  assume( P > 0 );
677
678  Print("//   %d parameter    : ", P);
679
680  for (int nop=0; nop < P; nop ++)
681    Print("%s ", rRingVar(nop, A));
682
683  assume( A->qideal == NULL );
684
685  PrintS("\n//   minpoly        : 0\n");
686
687/*
688  PrintS("//   Coefficients live in the rational function field\n");
689  Print("//   K(");
690  for (int i = 0; i < rVar(ntRing); i++)
691  {
692    if (i > 0) PrintS(" ");
693    Print("%s", rRingVar(i, ntRing));
694  }
695  PrintS(") with\n");
696  PrintS("//   K: "); n_CoeffWrite(cf->extRing->cf);
697*/
698}
699
700number ntAdd(number a, number b, const coeffs cf)
701{
702  ntTest(a);
703  ntTest(b);
704  if (IS0(a)) return ntCopy(b, cf);
705  if (IS0(b)) return ntCopy(a, cf);
706
707  fraction fa = (fraction)a;
708  fraction fb = (fraction)b;
709
710  poly g = p_Copy(NUM(fa), ntRing);
711  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
712  poly h = p_Copy(NUM(fb), ntRing);
713  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
714  g = p_Add_q(g, h, ntRing);
715
716  if (g == NULL) return NULL;
717
718  poly f;
719  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
720  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
721  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
722  else /* both denom's are != 1 */    f = p_Mult_q(p_Copy(DEN(fa), ntRing),
723                                                   p_Copy(DEN(fb), ntRing),
724                                                   ntRing);
725
726  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
727  NUM(result) = g;
728  DEN(result) = f;
729  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
730  heuristicGcdCancellation((number)result, cf);
731
732//  ntTest((number)result);
733
734  return (number)result;
735}
736
737number ntSub(number a, number b, const coeffs cf)
738{
739  ntTest(a);
740  ntTest(b);
741  if (IS0(a)) return ntNeg(ntCopy(b, cf), cf);
742  if (IS0(b)) return ntCopy(a, cf);
743
744  fraction fa = (fraction)a;
745  fraction fb = (fraction)b;
746
747  poly g = p_Copy(NUM(fa), ntRing);
748  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
749  poly h = p_Copy(NUM(fb), ntRing);
750  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
751  g = p_Add_q(g, p_Neg(h, ntRing), ntRing);
752
753  if (g == NULL) return NULL;
754
755  poly f;
756  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
757  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
758  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
759  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
760                                                   p_Copy(DEN(fb), ntRing),
761                                                   ntRing);
762
763  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
764  NUM(result) = g;
765  DEN(result) = f;
766  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
767  heuristicGcdCancellation((number)result, cf);
768//  ntTest((number)result);
769  return (number)result;
770}
771
772number ntMult(number a, number b, const coeffs cf)
773{
774  ntTest(a); // !!!?
775  ntTest(b); // !!!?
776
777  if (IS0(a) || IS0(b)) return NULL;
778
779  fraction fa = (fraction)a;
780  fraction fb = (fraction)b;
781
782  const poly g = pp_Mult_qq(NUM(fa), NUM(fb), ntRing);
783
784  if (g == NULL) return NULL; // may happen due to zero divisors???
785
786  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
787
788  NUM(result) = g;
789
790  const poly da = DEN(fa);
791  const poly db = DEN(fb);
792
793
794  if (db == NULL)
795  {
796    // b = ? // NULL
797
798    if(da == NULL)
799    { // both fa && fb are ?? // NULL!
800      assume (da == NULL && db == NULL);
801      DEN(result) = NULL;
802      COM(result) = 0;
803    }
804    else
805    {
806      assume (da != NULL && db == NULL);
807      DEN(result) = p_Copy(da, ntRing);
808      COM(result) = COM(fa) + MULT_COMPLEXITY;
809      heuristicGcdCancellation((number)result, cf);
810    }
811  } else
812  { // b = ?? / ??
813    if (da == NULL)
814    { // a == ? // NULL
815      assume( db != NULL && da == NULL);
816      DEN(result) = p_Copy(db, ntRing);
817      COM(result) = COM(fb) + MULT_COMPLEXITY;
818      heuristicGcdCancellation((number)result, cf);
819    }
820    else /* both den's are != 1 */
821    {
822      assume (da != NULL && db != NULL);
823      DEN(result) = pp_Mult_qq(da, db, ntRing);
824      COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
825      heuristicGcdCancellation((number)result, cf);
826    }
827  }
828
829//  ntTest((number)result);
830
831  return (number)result;
832}
833
834number ntDiv(number a, number b, const coeffs cf)
835{
836  ntTest(a);
837  ntTest(b);
838  if (IS0(a)) return NULL;
839  if (IS0(b)) WerrorS(nDivBy0);
840
841  fraction fa = (fraction)a;
842  fraction fb = (fraction)b;
843
844  poly g = p_Copy(NUM(fa), ntRing);
845  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
846
847  if (g == NULL) return NULL;   /* may happen due to zero divisors */
848
849  poly f = p_Copy(NUM(fb), ntRing);
850  if (!DENIS1(fa)) f = p_Mult_q(f, p_Copy(DEN(fa), ntRing), ntRing);
851
852  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
853  NUM(result) = g;
854  if (!n_GreaterZero(pGetCoeff(f),ntCoeffs))
855  {
856    g=p_Neg(g,ntRing);
857    f=p_Neg(f,ntRing);
858    NUM(result) = g;
859  }
860  if (!p_IsConstant(f,ntRing) || !n_IsOne(pGetCoeff(f),ntCoeffs))
861  {
862    DEN(result) = f;
863  }
864  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
865  heuristicGcdCancellation((number)result, cf);
866//  ntTest((number)result);
867  return (number)result;
868}
869
870/* 0^0 = 0;
871   for |exp| <= 7 compute power by a simple multiplication loop;
872   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
873      p^13 = p^1 * p^4 * p^8, where we utilise that
874      p^(2^(k+1)) = p^(2^k) * p^(2^k);
875   intermediate cancellation is controlled by the in-place method
876   heuristicGcdCancellation; see there.
877*/
878void ntPower(number a, int exp, number *b, const coeffs cf)
879{
880  ntTest(a);
881
882  /* special cases first */
883  if (IS0(a))
884  {
885    if (exp >= 0) *b = NULL;
886    else          WerrorS(nDivBy0);
887  }
888  else if (exp ==  0) { *b = ntInit(1, cf); return;}
889  else if (exp ==  1) { *b = ntCopy(a, cf); return;}
890  else if (exp == -1) { *b = ntInvers(a, cf); return;}
891
892  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
893
894  /* now compute a^expAbs */
895  number pow; number t;
896  if (expAbs <= 7)
897  {
898    pow = ntCopy(a, cf);
899    for (int i = 2; i <= expAbs; i++)
900    {
901      t = ntMult(pow, a, cf);
902      ntDelete(&pow, cf);
903      pow = t;
904      heuristicGcdCancellation(pow, cf);
905    }
906  }
907  else
908  {
909    pow = ntInit(1, cf);
910    number factor = ntCopy(a, cf);
911    while (expAbs != 0)
912    {
913      if (expAbs & 1)
914      {
915        t = ntMult(pow, factor, cf);
916        ntDelete(&pow, cf);
917        pow = t;
918        heuristicGcdCancellation(pow, cf);
919      }
920      expAbs = expAbs / 2;
921      if (expAbs != 0)
922      {
923        t = ntMult(factor, factor, cf);
924        ntDelete(&factor, cf);
925        factor = t;
926        heuristicGcdCancellation(factor, cf);
927      }
928    }
929    ntDelete(&factor, cf);
930  }
931
932  /* invert if original exponent was negative */
933  if (exp < 0)
934  {
935    t = ntInvers(pow, cf);
936    ntDelete(&pow, cf);
937    pow = t;
938  }
939  *b = pow;
940  ntTest(*b);
941}
942
943/* assumes that cf represents the rationals, i.e. Q, and will only
944   be called in that case;
945   assumes furthermore that f != NULL and that the denominator of f != 1;
946   generally speaking, this method removes denominators in the rational
947   coefficients of the numerator and denominator of 'a';
948   more concretely, the following normalizations will be performed,
949   where t^alpha denotes a monomial in the transcendental variables t_k
950   (1) if 'a' is of the form
951          (sum_alpha a_alpha/b_alpha * t^alpha)
952          -------------------------------------
953            (sum_beta c_beta/d_beta * t^beta)
954       with integers a_alpha, b_alpha, c_beta, d_beta, then both the
955       numerator and the denominator will be multiplied by the LCM of
956       the b_alpha's and the d_beta's (if this LCM is != 1),
957   (2) if 'a' is - e.g. after having performed step (1) - of the form
958          (sum_alpha a_alpha * t^alpha)
959          -----------------------------
960            (sum_beta c_beta * t^beta)
961       with integers a_alpha, c_beta, and with a non-constant denominator,
962       then both the numerator and the denominator will be divided by the
963       GCD of the a_alpha's and the c_beta's (if this GCD is != 1),
964   (3) if 'a' is - e.g. after having performed steps (1) and (2) - of the
965       form
966          (sum_alpha a_alpha * t^alpha)
967          -----------------------------
968                        c
969       with integers a_alpha, and c != 1, then 'a' will be replaced by
970       (sum_alpha a_alpha/c * t^alpha);
971   this procedure does not alter COM(f) (this has to be done by the
972   calling procedure);
973   modifies f */
974void handleNestedFractionsOverQ(fraction f, const coeffs cf)
975{
976  assume(nCoeff_is_Q(ntCoeffs));
977  assume(!IS0(f));
978  assume(!DENIS1(f));
979
980  if (!p_IsConstant(DEN(f), ntRing))
981  { /* step (1); see documentation of this procedure above */
982    p_Normalize(NUM(f), ntRing);
983    p_Normalize(DEN(f), ntRing);
984    number lcmOfDenominators = n_Init(1, ntCoeffs);
985    number c; number tmp;
986    poly p = NUM(f);
987    /* careful when using n_Lcm!!! It computes the lcm of the numerator
988       of the 1st argument and the denominator of the 2nd!!! */
989    while (p != NULL)
990    {
991      c = p_GetCoeff(p, ntRing);
992      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
993      n_Delete(&lcmOfDenominators, ntCoeffs);
994      lcmOfDenominators = tmp;
995      pIter(p);
996    }
997    p = DEN(f);
998    while (p != NULL)
999    {
1000      c = p_GetCoeff(p, ntRing);
1001      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
1002      n_Delete(&lcmOfDenominators, ntCoeffs);
1003      lcmOfDenominators = tmp;
1004      pIter(p);
1005    }
1006    if (!n_IsOne(lcmOfDenominators, ntCoeffs))
1007    { /* multiply NUM(f) and DEN(f) with lcmOfDenominators */
1008      NUM(f) = p_Mult_nn(NUM(f), lcmOfDenominators, ntRing);
1009      p_Normalize(NUM(f), ntRing);
1010      DEN(f) = p_Mult_nn(DEN(f), lcmOfDenominators, ntRing);
1011      p_Normalize(DEN(f), ntRing);
1012    }
1013    n_Delete(&lcmOfDenominators, ntCoeffs);
1014    if (!p_IsConstant(DEN(f), ntRing))
1015    { /* step (2); see documentation of this procedure above */
1016      p = NUM(f);
1017      number gcdOfCoefficients = n_Copy(p_GetCoeff(p, ntRing), ntCoeffs);
1018      pIter(p);
1019      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
1020      {
1021        c = p_GetCoeff(p, ntRing);
1022        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
1023        n_Delete(&gcdOfCoefficients, ntCoeffs);
1024        gcdOfCoefficients = tmp;
1025        pIter(p);
1026      }
1027      p = DEN(f);
1028      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
1029      {
1030        c = p_GetCoeff(p, ntRing);
1031        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
1032        n_Delete(&gcdOfCoefficients, ntCoeffs);
1033        gcdOfCoefficients = tmp;
1034        pIter(p);
1035      }
1036      if (!n_IsOne(gcdOfCoefficients, ntCoeffs))
1037      { /* divide NUM(f) and DEN(f) by gcdOfCoefficients */
1038        number inverseOfGcdOfCoefficients = n_Invers(gcdOfCoefficients,
1039                                                     ntCoeffs);
1040        NUM(f) = p_Mult_nn(NUM(f), inverseOfGcdOfCoefficients, ntRing);
1041        p_Normalize(NUM(f), ntRing);
1042        DEN(f) = p_Mult_nn(DEN(f), inverseOfGcdOfCoefficients, ntRing);
1043        p_Normalize(DEN(f), ntRing);
1044        n_Delete(&inverseOfGcdOfCoefficients, ntCoeffs);
1045      }
1046      n_Delete(&gcdOfCoefficients, ntCoeffs);
1047    }
1048  }
1049  if (p_IsConstant(DEN(f), ntRing) &&
1050      (!n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs)))
1051  { /* step (3); see documentation of this procedure above */
1052    number inverseOfDen = n_Invers(p_GetCoeff(DEN(f), ntRing), ntCoeffs);
1053    NUM(f) = p_Mult_nn(NUM(f), inverseOfDen, ntRing);
1054    n_Delete(&inverseOfDen, ntCoeffs);
1055    p_Delete(&DEN(f), ntRing);
1056    DEN(f) = NULL;
1057  }
1058
1059  /* Now, due to the above computations, DEN(f) may have become the
1060     1-polynomial which needs to be represented by NULL: */
1061  if ((DEN(f) != NULL) &&
1062      p_IsConstant(DEN(f), ntRing) &&
1063      n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
1064  {
1065    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1066  }
1067
1068  if( DEN(f) != NULL )
1069    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1070    {
1071      NUM(f) = p_Neg(NUM(f), ntRing);
1072      DEN(f) = p_Neg(DEN(f), ntRing);
1073    }
1074
1075  ntTest((number)f); // TODO!
1076}
1077
1078/* modifies a */
1079void heuristicGcdCancellation(number a, const coeffs cf)
1080{
1081//  ntTest(a); // !!!!????
1082  if (IS0(a)) return;
1083
1084  fraction f = (fraction)a;
1085  if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; ntTest(a); return; }
1086
1087  assume( DEN(f) != NULL );
1088
1089  /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
1090  if (p_EqualPolys(NUM(f), DEN(f), ntRing))
1091  { /* numerator and denominator are both != 1 */
1092    p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
1093    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1094    COM(f) = 0;
1095  }
1096  else
1097  {
1098    if (COM(f) > BOUND_COMPLEXITY)
1099      definiteGcdCancellation(a, cf, TRUE);
1100
1101    // TODO: check if it is enough to put the following into definiteGcdCancellation?!
1102    if( DEN(f) != NULL )
1103    {
1104      if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1105      {
1106        NUM(f) = p_Neg(NUM(f), ntRing);
1107        DEN(f) = p_Neg(DEN(f), ntRing);
1108      }
1109      if (ntCoeffs->has_simple_Inverse)
1110      {
1111        if (!n_IsOne(pGetCoeff(DEN(f)),ntCoeffs))
1112        {
1113          number inv=n_Invers(pGetCoeff(DEN(f)),ntCoeffs);
1114          DEN(f)=p_Mult_nn(DEN(f),inv,ntRing);
1115          NUM(f)=p_Mult_nn(NUM(f),inv,ntRing);
1116        }
1117        if(p_LmIsConstant(DEN(f),ntRing))
1118        {
1119          p_Delete(&DEN(f),ntRing);
1120          COM(f)=0;
1121        }
1122      }
1123    }
1124  }
1125
1126  ntTest(a);
1127}
1128
1129/// modifies a
1130void definiteGcdCancellation(number a, const coeffs cf,
1131                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed)
1132{
1133  ntTest(a); // !!!!
1134
1135  fraction f = (fraction)a;
1136
1137  if (IS0(a)) return;
1138  if (NUM(f)!=NULL) p_Normalize(NUM(f), ntRing);
1139  if (DEN(f)!=NULL) p_Normalize(DEN(f), ntRing);
1140  if (!simpleTestsHaveAlreadyBeenPerformed)
1141  {
1142    if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
1143
1144    /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
1145    if (p_EqualPolys(NUM(f), DEN(f), ntRing))
1146    { /* numerator and denominator are both != 1 */
1147      p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
1148      p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1149      COM(f) = 0;
1150      ntTest(a); // !!!!
1151      return;
1152    }
1153  }
1154
1155#ifdef HAVE_FACTORY
1156  /* Note that, over Q, singclap_gcd will remove the denominators in all
1157     rational coefficients of pNum and pDen, before starting to compute
1158     the gcd. Thus, we do not need to ensure that the coefficients of
1159     pNum and pDen live in Z; they may well be elements of Q\Z. */
1160  /* singclap_gcd destroys its arguments; we hence need copies: */
1161  poly pGcd = singclap_gcd(p_Copy(NUM(f), ntRing), p_Copy(DEN(f), ntRing), cf->extRing);
1162  if (p_IsConstant(pGcd, ntRing) &&
1163      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
1164  { /* gcd = 1; nothing to cancel;
1165       Suppose the given rational function field is over Q. Although the
1166       gcd is 1, we may have produced fractional coefficients in NUM(f),
1167       DEN(f), or both, due to previous arithmetics. The next call will
1168       remove those nested fractions, in case there are any. */
1169    if (nCoeff_is_Zp(ntCoeffs) && p_IsConstant (DEN (f), ntRing))
1170    {
1171      NUM (f) = p_Div_nn (NUM (f), p_GetCoeff (DEN(f),ntRing), ntRing);
1172      //poly newNum= singclap_pdivide (NUM(f), DEN (f), ntRing);
1173      //p_Delete(&NUM (f), ntRing);
1174      //NUM (f)= newNum;
1175      p_Delete(&DEN (f), ntRing);
1176      DEN (f) = NULL;
1177      COM (f) = 0;
1178    } else if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
1179  }
1180  else
1181  { /* We divide both NUM(f) and DEN(f) by the gcd which is known
1182       to be != 1. */
1183    poly newNum = singclap_pdivide(NUM(f), pGcd, ntRing);
1184    p_Delete(&NUM(f), ntRing);
1185    NUM(f) = newNum;
1186    poly newDen = singclap_pdivide(DEN(f), pGcd, ntRing);
1187    p_Delete(&DEN(f), ntRing);
1188    DEN(f) = newDen;
1189    if (p_IsConstant(DEN(f), ntRing) &&
1190        n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
1191    {
1192      /* DEN(f) = 1 needs to be represented by NULL! */
1193      p_Delete(&DEN(f), ntRing);
1194      newDen = NULL;
1195    }
1196    else
1197    { /* Note that over Q, by cancelling the gcd, we may have produced
1198         fractional coefficients in NUM(f), DEN(f), or both. The next
1199         call will remove those nested fractions, in case there are
1200         any. */
1201      if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
1202    }
1203  }
1204  COM(f) = 0;
1205  p_Delete(&pGcd, ntRing);
1206
1207  if( DEN(f) != NULL )
1208    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1209    {
1210      NUM(f) = p_Neg(NUM(f), ntRing);
1211      DEN(f) = p_Neg(DEN(f), ntRing);
1212    }
1213#endif /* HAVE_FACTORY */
1214
1215  ntTest(a); // !!!!
1216}
1217
1218// NOTE: modifies a
1219void ntWriteLong(number &a, const coeffs cf)
1220{
1221  ntTest(a);
1222  definiteGcdCancellation(a, cf, FALSE);
1223  if (IS0(a))
1224    StringAppendS("0");
1225  else
1226  {
1227    fraction f = (fraction)a;
1228    // stole logic from napWrite from kernel/longtrans.cc of legacy singular
1229    BOOLEAN omitBrackets = p_IsConstant(NUM(f), ntRing);
1230    if (!omitBrackets) StringAppendS("(");
1231    p_String0Long(NUM(f), ntRing, ntRing);
1232    if (!omitBrackets) StringAppendS(")");
1233    if (!DENIS1(f))
1234    {
1235      StringAppendS("/");
1236      omitBrackets = p_IsConstant(DEN(f), ntRing);
1237      if (!omitBrackets) StringAppendS("(");
1238      p_String0Long(DEN(f), ntRing, ntRing);
1239      if (!omitBrackets) StringAppendS(")");
1240    }
1241  }
1242  ntTest(a); // !!!!
1243}
1244
1245// NOTE: modifies a
1246void ntWriteShort(number &a, const coeffs cf)
1247{
1248  ntTest(a);
1249  definiteGcdCancellation(a, cf, FALSE);
1250  if (IS0(a))
1251    StringAppendS("0");
1252  else
1253  {
1254    fraction f = (fraction)a;
1255    // stole logic from napWrite from kernel/longtrans.cc of legacy singular
1256    BOOLEAN omitBrackets = p_IsConstant(NUM(f), ntRing);
1257    if (!omitBrackets) StringAppendS("(");
1258    p_String0Short(NUM(f), ntRing, ntRing);
1259    if (!omitBrackets) StringAppendS(")");
1260    if (!DENIS1(f))
1261    {
1262      StringAppendS("/");
1263      omitBrackets = p_IsConstant(DEN(f), ntRing);
1264      if (!omitBrackets) StringAppendS("(");
1265      p_String0Short(DEN(f), ntRing, ntRing);
1266      if (!omitBrackets) StringAppendS(")");
1267    }
1268  }
1269  ntTest(a);
1270}
1271
1272const char * ntRead(const char *s, number *a, const coeffs cf)
1273{
1274  poly p;
1275  const char * result = p_Read(s, p, ntRing);
1276  if (p == NULL) *a = NULL;
1277  else *a = ntInit(p, cf);
1278  return result;
1279}
1280
1281void ntNormalize (number &a, const coeffs cf)
1282{
1283  if ((a!=NULL))
1284  {
1285    definiteGcdCancellation(a, cf, FALSE);
1286    if ((DEN(a)!=NULL)
1287    &&(!n_GreaterZero(pGetCoeff(DEN(a)),ntCoeffs)))
1288    {
1289      NUM(a)=p_Neg(NUM(a),ntRing);
1290      DEN(a)=p_Neg(DEN(a),ntRing);
1291    }
1292  }
1293  ntTest(a); // !!!!
1294}
1295
1296/* expects *param to be castable to TransExtInfo */
1297static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
1298{
1299  if (ID != n) return FALSE;
1300  TransExtInfo *e = (TransExtInfo *)param;
1301  /* for rational function fields we expect the underlying
1302     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
1303     this expectation is based on the assumption that we have properly
1304     registered cf and perform reference counting rather than creating
1305     multiple copies of the same coefficient field/domain/ring */
1306  if (ntRing == e->r)
1307    return TRUE;
1308
1309  // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)...
1310  if( rEqual(ntRing, e->r, TRUE) )
1311  {
1312    rDelete(e->r);
1313    return TRUE;
1314  }
1315
1316  return FALSE;
1317}
1318
1319number ntLcm(number a, number b, const coeffs cf)
1320{
1321  ntTest(a);
1322  ntTest(b);
1323  fraction fb = (fraction)b;
1324  if ((b==NULL)||(DEN(fb)==NULL)) return ntCopy(a,cf);
1325#ifdef HAVE_FACTORY
1326  fraction fa = (fraction)a;
1327  /* singclap_gcd destroys its arguments; we hence need copies: */
1328  poly pa = p_Copy(NUM(fa), ntRing);
1329  poly pb = p_Copy(DEN(fb), ntRing);
1330
1331  /* Note that, over Q, singclap_gcd will remove the denominators in all
1332     rational coefficients of pa and pb, before starting to compute
1333     the gcd. Thus, we do not need to ensure that the coefficients of
1334     pa and pb live in Z; they may well be elements of Q\Z. */
1335  poly pGcd = singclap_gcd(pa, pb, cf->extRing);
1336  if (p_IsConstant(pGcd, ntRing) &&
1337      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
1338  { /* gcd = 1; return pa*pb*/
1339    p_Delete(&pGcd,ntRing);
1340    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1341    NUM(result) = pp_Mult_qq(NUM(fa),DEN(fb),ntRing);
1342
1343    ntTest((number)result); // !!!!
1344
1345    return (number)result;
1346  }
1347
1348
1349  /* return pa*pb/gcd */
1350    poly newNum = singclap_pdivide(NUM(fa), pGcd, ntRing);
1351    p_Delete(&pGcd,ntRing);
1352    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1353    NUM(result) = p_Mult_q(p_Copy(DEN(fb),ntRing),newNum,ntRing);
1354    ntTest((number)result); // !!!!
1355    return (number)result;
1356
1357#else
1358  Print("// factory needed: transext.cc:ntLcm\n");
1359  return NULL;
1360#endif /* HAVE_FACTORY */
1361  return NULL;
1362}
1363
1364number ntGcd(number a, number b, const coeffs cf)
1365{
1366  ntTest(a);
1367  ntTest(b);
1368  if (a==NULL) return ntCopy(b,cf);
1369  if (b==NULL) return ntCopy(a,cf);
1370#ifdef HAVE_FACTORY
1371  fraction fa = (fraction)a;
1372  fraction fb = (fraction)b;
1373  /* singclap_gcd destroys its arguments; we hence need copies: */
1374  poly pa = p_Copy(NUM(fa), ntRing);
1375  poly pb = p_Copy(NUM(fb), ntRing);
1376
1377  /* Note that, over Q, singclap_gcd will remove the denominators in all
1378     rational coefficients of pa and pb, before starting to compute
1379     the gcd. Thus, we do not need to ensure that the coefficients of
1380     pa and pb live in Z; they may well be elements of Q\Z. */
1381  poly pGcd = singclap_gcd(pa, pb, cf->extRing);
1382  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1383  NUM(result) = pGcd;
1384  ntTest((number)result); // !!!!
1385  return (number)result;
1386#else
1387  Print("// factory needed: transext.cc:ntGcd\n");
1388  return NULL;
1389#endif /* HAVE_FACTORY */
1390}
1391
1392int ntSize(number a, const coeffs cf)
1393{
1394  ntTest(a);
1395  if (IS0(a)) return -1;
1396  /* this has been taken from the old implementation of field extensions,
1397     where we computed the sum of the degrees and the numbers of terms in
1398     the numerator and denominator of a; so we leave it at that, for the
1399     time being */
1400  fraction f = (fraction)a;
1401  poly p = NUM(f);
1402  int noOfTerms = 0;
1403  int numDegree = 0;
1404  while (p != NULL)
1405  {
1406    noOfTerms++;
1407    int d = 0;
1408    for (int i = 1; i <= rVar(ntRing); i++)
1409      d += p_GetExp(p, i, ntRing);
1410    if (d > numDegree) numDegree = d;
1411    pIter(p);
1412  }
1413  int denDegree = 0;
1414  if (!DENIS1(f))
1415  {
1416    p = DEN(f);
1417    while (p != NULL)
1418    {
1419      noOfTerms++;
1420      int d = 0;
1421      for (int i = 1; i <= rVar(ntRing); i++)
1422        d += p_GetExp(p, i, ntRing);
1423      if (d > denDegree) denDegree = d;
1424      pIter(p);
1425    }
1426  }
1427  ntTest(a); // !!!!
1428  return numDegree + denDegree + noOfTerms;
1429}
1430
1431number ntInvers(number a, const coeffs cf)
1432{
1433  ntTest(a);
1434  if (IS0(a))
1435  {
1436    WerrorS(nDivBy0);
1437    return NULL;
1438  }
1439  fraction f = (fraction)a;
1440  assume( f != NULL );
1441
1442  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1443
1444  assume( NUM(f) != NULL );
1445  const poly den = DEN(f);
1446
1447  if (den == NULL)
1448    NUM(result) = p_One(ntRing);
1449  else
1450    NUM(result) = p_Copy(den, ntRing);
1451
1452  if( !NUMIS1(f) )
1453  {
1454    poly num_f=NUM(f);
1455    BOOLEAN neg= !n_GreaterZero(pGetCoeff(num_f),ntCoeffs);
1456    if (neg)
1457    {
1458      num_f=p_Neg(p_Copy(num_f, ntRing), ntRing);
1459      NUM(result)=p_Neg(NUM(result), ntRing);
1460    }
1461    else
1462    {
1463      num_f=p_Copy(num_f, ntRing);
1464    }
1465    DEN(result) = num_f;
1466    COM(result) = COM(f);
1467    if (neg)
1468    {
1469      if (p_IsOne(num_f, ntRing))
1470      {
1471        DEN(result)=NULL;
1472        COM(result) = 0;
1473        p_Delete(&num_f,ntRing);
1474      }
1475    }
1476  }
1477  else
1478  {
1479    DEN(result) = NULL;
1480    COM(result) = 0;
1481  }
1482  ntTest((number)result); // !!!!
1483  return (number)result;
1484}
1485
1486/* assumes that src = Q, dst = Q(t_1, ..., t_s) */
1487number ntMap00(number a, const coeffs src, const coeffs dst)
1488{
1489  if (n_IsZero(a, src)) return NULL;
1490  assume(n_Test(a, src));
1491  assume(src == dst->extRing->cf);
1492  if ((SR_HDL(a) & SR_INT) || (a->s==3))
1493  {
1494    number res=ntInit(p_NSet(n_Copy(a, src), dst->extRing), dst);
1495    n_Test(res,dst);
1496    return res;
1497  }
1498  number nn=nlGetDenom(a,src);
1499  number zz=nlGetNumerator(a,src);
1500  number res=ntInit(p_NSet(zz,dst->extRing), dst);
1501  fraction ff=(fraction)res;
1502  if (n_IsOne(nn,src)) DEN(ff)=NULL;
1503  else                 DEN(ff)=p_NSet(nn,dst->extRing);
1504  n_Test((number)ff,dst);
1505  return (number)ff;
1506}
1507
1508/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
1509number ntMapP0(number a, const coeffs src, const coeffs dst)
1510{
1511  if (n_IsZero(a, src)) return NULL;
1512  assume(n_Test(a, src));
1513  /* mapping via intermediate int: */
1514  int n = n_Int(a, src);
1515  number q = n_Init(n, dst->extRing->cf);
1516  if (n_IsZero(q, dst->extRing->cf))
1517  {
1518    n_Delete(&q, dst->extRing->cf);
1519    return NULL;
1520  }
1521  return ntInit(p_NSet(q, dst->extRing), dst);
1522}
1523
1524/* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
1525                       src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
1526number ntCopyMap(number a, const coeffs cf, const coeffs dst)
1527{
1528//  if (n_IsZero(a, cf)) return NULL;
1529
1530  ntTest(a);
1531
1532  if (IS0(a)) return NULL;
1533
1534  const ring rSrc = cf->extRing;
1535  const ring rDst = dst->extRing;
1536
1537  if( rSrc == rDst )
1538    return ntCopy(a, dst); // USUALLY WRONG!
1539
1540  fraction f = (fraction)a;
1541  poly g = prCopyR(NUM(f), rSrc, rDst);
1542
1543  poly h = NULL;
1544
1545  if (!DENIS1(f))
1546     h = prCopyR(DEN(f), rSrc, rDst);
1547
1548  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1549
1550  NUM(result) = g;
1551  DEN(result) = h;
1552  COM(result) = COM(f);
1553  assume(n_Test((number)result, dst));
1554  return (number)result;
1555}
1556
1557number ntCopyAlg(number a, const coeffs cf, const coeffs dst)
1558{
1559  assume( n_Test(a, cf) );
1560  if (n_IsZero(a, cf)) return NULL;
1561
1562  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1563  // DEN(f) = NULL; COM(f) = 0;
1564  NUM(f) = prCopyR((poly)a, cf->extRing, dst->extRing);
1565  assume(n_Test((number)f, dst));
1566  return (number)f;
1567}
1568
1569/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
1570number ntMap0P(number a, const coeffs src, const coeffs dst)
1571{
1572  assume( n_Test(a, src) );
1573  if (n_IsZero(a, src)) return NULL;
1574  // int p = rChar(dst->extRing);
1575  number q = nlModP(a, src, dst->extRing->cf);
1576
1577  if (n_IsZero(q, dst->extRing->cf))
1578  {
1579    n_Delete(&q, dst->extRing->cf);
1580    return NULL;
1581  }
1582
1583  poly g = p_NSet(q, dst->extRing);
1584
1585  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1586  NUM(f) = g; // DEN(f) = NULL; COM(f) = 0;
1587  assume(n_Test((number)f, dst));
1588  return (number)f;
1589}
1590
1591/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
1592number ntMapPP(number a, const coeffs src, const coeffs dst)
1593{
1594  assume( n_Test(a, src) );
1595  if (n_IsZero(a, src)) return NULL;
1596  assume(src == dst->extRing->cf);
1597  poly p = p_One(dst->extRing);
1598  p_SetCoeff(p, n_Copy(a, src), dst->extRing);
1599  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1600  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1601  assume(n_Test((number)f, dst));
1602  return (number)f;
1603}
1604
1605/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
1606number ntMapUP(number a, const coeffs src, const coeffs dst)
1607{
1608  assume( n_Test(a, src) );
1609  if (n_IsZero(a, src)) return NULL;
1610  /* mapping via intermediate int: */
1611  int n = n_Int(a, src);
1612  number q = n_Init(n, dst->extRing->cf);
1613  poly p;
1614  if (n_IsZero(q, dst->extRing->cf))
1615  {
1616    n_Delete(&q, dst->extRing->cf);
1617    return NULL;
1618  }
1619  p = p_One(dst->extRing);
1620  p_SetCoeff(p, q, dst->extRing);
1621  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1622  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1623  assume(n_Test((number)f, dst));
1624  return (number)f;
1625}
1626
1627nMapFunc ntSetMap(const coeffs src, const coeffs dst)
1628{
1629  /* dst is expected to be a rational function field */
1630  assume(getCoeffType(dst) == ID);
1631
1632  if( src == dst ) return ndCopyMap;
1633
1634  int h = 0; /* the height of the extension tower given by dst */
1635  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
1636  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
1637
1638  /* for the time being, we only provide maps if h = 1 and if b is Q or
1639     some field Z/pZ: */
1640  if (h==0)
1641  {
1642    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
1643      return ntMap00;                                 /// Q       -->  Q(T)
1644    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
1645      return ntMapP0;                                 /// Z/p     -->  Q(T)
1646    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
1647      return ntMap0P;                                 /// Q       --> Z/p(T)
1648    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
1649    {
1650      if (src->ch == dst->ch) return ntMapPP;         /// Z/p     --> Z/p(T)
1651      else return ntMapUP;                            /// Z/u     --> Z/p(T)
1652    }
1653  }
1654  if (h != 1) return NULL;
1655  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
1656
1657  /* Let T denote the sequence of transcendental extension variables, i.e.,
1658     K[t_1, ..., t_s] =: K[T];
1659     Let moreover, for any such sequence T, T' denote any subsequence of T
1660     of the form t_1, ..., t_w with w <= s. */
1661
1662  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
1663
1664  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
1665  {
1666    if (rVar(src->extRing) > rVar(dst->extRing))
1667       return NULL;
1668
1669    for (int i = 0; i < rVar(src->extRing); i++)
1670      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0)
1671         return NULL;
1672
1673    if (src->type==n_transExt)
1674       return ntCopyMap;          /// Q(T')   --> Q(T)
1675    else
1676       return ntCopyAlg;
1677  }
1678
1679  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
1680  {
1681    if (rVar(src->extRing) > rVar(dst->extRing))
1682       return NULL;
1683
1684    for (int i = 0; i < rVar(src->extRing); i++)
1685      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0)
1686         return NULL;
1687
1688    if (src->type==n_transExt)
1689       return ntCopyMap;         /// Z/p(T') --> Z/p(T)
1690    else
1691       return ntCopyAlg;
1692  }
1693
1694  return NULL;                                 /// default
1695}
1696#if 0
1697nMapFunc ntSetMap_T(const coeffs src, const coeffs dst)
1698{
1699  nMapFunc n=ntSetMap(src,dst);
1700  if (n==ntCopyAlg) printf("n=ntCopyAlg\n");
1701  else if (n==ntCopyMap) printf("n=ntCopyMap\n");
1702  else if (n==ntMapUP) printf("n=ntMapUP\n");
1703  else if (n==ntMap0P) printf("n=ntMap0P\n");
1704  else if (n==ntMapP0) printf("n=ntMapP0\n");
1705  else if (n==ntMap00) printf("n=ntMap00\n");
1706  else if (n==NULL) printf("n=NULL\n");
1707  else printf("n=?\n");
1708  return n;
1709}
1710#endif
1711
1712void ntKillChar(coeffs cf)
1713{
1714  if ((--cf->extRing->ref) == 0)
1715    rDelete(cf->extRing);
1716}
1717#ifdef HAVE_FACTORY
1718number ntConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
1719{
1720  if (n.isZero()) return NULL;
1721  poly p=convFactoryPSingP(n,ntRing);
1722  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1723  NUM(result) = p;
1724  //DEN(result) = NULL; // done by omAlloc0Bin
1725  //COM(result) = 0; // done by omAlloc0Bin
1726  ntTest((number)result);
1727  return (number)result;
1728}
1729CanonicalForm ntConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
1730{
1731  ntTest(n);
1732  if (IS0(n)) return CanonicalForm(0);
1733
1734  fraction f = (fraction)n;
1735  return convSingPFactoryP(NUM(f),ntRing);
1736}
1737#endif
1738
1739static int ntParDeg(number a, const coeffs cf)
1740{
1741  ntTest(a);
1742  if (IS0(a)) return -1;
1743  fraction fa = (fraction)a;
1744  return cf->extRing->pFDeg(NUM(fa),cf->extRing);
1745}
1746
1747/// return the specified parameter as a number in the given trans.ext.
1748static number ntParameter(const int iParameter, const coeffs cf)
1749{
1750  assume(getCoeffType(cf) == ID);
1751
1752  const ring R = cf->extRing;
1753  assume( R != NULL );
1754  assume( 0 < iParameter && iParameter <= rVar(R) );
1755
1756  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
1757  p_Test(p,R);
1758
1759  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1760  NUM(f) = p;
1761  DEN(f) = NULL;
1762  COM(f) = 0;
1763
1764  ntTest((number)f);
1765
1766  return (number)f;
1767}
1768
1769/// if m == var(i)/1 => return i,
1770int ntIsParam(number m, const coeffs cf)
1771{
1772  ntTest(m);
1773  assume(getCoeffType(cf) == ID);
1774
1775  const ring R = cf->extRing;
1776  assume( R != NULL );
1777
1778  fraction f = (fraction)m;
1779
1780  if( DEN(f) != NULL )
1781    return 0;
1782
1783  return p_Var( NUM(f), R );
1784}
1785
1786struct NTNumConverter
1787{
1788  static inline poly convert(const number& n)
1789  {
1790    // suitable for trans. ext. numbers that are fractions of polys
1791    return NUM((fraction)n); // return the numerator
1792  }
1793};
1794
1795
1796static void ntClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1797{
1798  assume(cf != NULL);
1799  assume(getCoeffType(cf) == ID);
1800  // all coeffs are given by fractions of polynomails over integers!!!
1801  // without denominators!!!
1802
1803  const ring   R = cf->extRing;
1804  assume(R != NULL);
1805  const coeffs Q = R->cf;
1806  assume(Q != NULL);
1807  assume(nCoeff_is_Q(Q));
1808
1809
1810  numberCollectionEnumerator.Reset();
1811
1812  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
1813  {
1814    c = ntInit(1, cf);
1815    return;
1816  }
1817
1818  // all coeffs are given by integers after returning from this routine
1819
1820  // part 1, collect product of all denominators /gcds
1821  poly cand = NULL;
1822
1823  do
1824  {
1825    number &n = numberCollectionEnumerator.Current();
1826
1827    ntNormalize(n, cf);
1828
1829    fraction f = (fraction)n;
1830
1831    assume( f != NULL );
1832
1833    const poly den = DEN(f);
1834
1835    assume( den == NULL ); // ?? / 1 ?
1836
1837    const poly num = NUM(f);
1838
1839    if( cand == NULL )
1840      cand = p_Copy(num, R);
1841    else
1842      cand = singclap_gcd(cand, p_Copy(num, R), R); // gcd(cand, num)
1843
1844    if( p_IsConstant(cand, R) )
1845      break;
1846  }
1847  while( numberCollectionEnumerator.MoveNext() ) ;
1848
1849
1850  // part2: all coeffs = all coeffs * cand
1851  if( cand != NULL )
1852  {
1853  if( !p_IsConstant(cand, R) )
1854  {
1855    c = ntInit(cand, cf);
1856    numberCollectionEnumerator.Reset();
1857    while (numberCollectionEnumerator.MoveNext() )
1858    {
1859      number &n = numberCollectionEnumerator.Current();
1860      const number t = ntDiv(n, c, cf); // TODO: rewrite!?
1861      ntDelete(&n, cf);
1862      n = t;
1863    }
1864  } // else NUM (result) = p_One(R);
1865  else { p_Delete(&cand, R); cand = NULL; }
1866  }
1867
1868  // Quick and dirty fix for constant content clearing: consider numerators???
1869  CRecursivePolyCoeffsEnumerator<NTNumConverter> itr(numberCollectionEnumerator); // recursively treat the NUM(numbers) as polys!
1870  number cc;
1871
1872  n_ClearContent(itr, cc, Q);
1873  number g = ntInit(p_NSet(cc, R), cf);
1874
1875  if( cand != NULL )
1876  {
1877    number gg = ntMult(g, c, cf);
1878    ntDelete(&g, cf);
1879    ntDelete(&c, cf); c = gg;
1880  } else
1881    c = g;
1882  ntTest(c);
1883}
1884
1885static void ntClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1886{
1887  assume(cf != NULL);
1888  assume(getCoeffType(cf) == ID); // both over Q(a) and Zp(a)!
1889  // all coeffs are given by fractions of polynomails over integers!!!
1890
1891  numberCollectionEnumerator.Reset();
1892
1893  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
1894  {
1895    c = ntInit(1, cf);
1896    return;
1897  }
1898
1899  // all coeffs are given by integers after returning from this routine
1900
1901  // part 1, collect product of all denominators /gcds
1902  poly cand = NULL;
1903
1904  const ring R = cf->extRing;
1905  assume(R != NULL);
1906
1907  const coeffs Q = R->cf;
1908  assume(Q != NULL);
1909//  assume(nCoeff_is_Q(Q));
1910
1911  do
1912  {
1913    number &n = numberCollectionEnumerator.Current();
1914
1915    ntNormalize(n, cf);
1916
1917    fraction f = (fraction)ntGetDenom (n, cf);
1918
1919    assume( f != NULL );
1920
1921    const poly den = NUM(f);
1922
1923    if( den == NULL ) // ?? / 1 ?
1924      continue;
1925
1926    if( cand == NULL )
1927      cand = p_Copy(den, R);
1928    else
1929    {
1930      // cand === LCM( cand, den )!!!!
1931      // NOTE: maybe it's better to make the product and clearcontent afterwards!?
1932      // TODO: move the following to factory?
1933      poly gcd = singclap_gcd(p_Copy(cand, R), p_Copy(den, R), R); // gcd(cand, den) is monic no mater leading coeffs! :((((
1934      if (nCoeff_is_Q (Q))
1935      {
1936        number LcGcd= n_Gcd (p_GetCoeff (cand, R), p_GetCoeff(den, R), Q);
1937        gcd = p_Mult_nn(gcd, LcGcd, R);
1938        n_Delete(&LcGcd,Q);
1939      }
1940//      assume( n_IsOne(pGetCoeff(gcd), Q) ); // TODO: this may be wrong...
1941      cand = p_Mult_q(cand, p_Copy(den, R), R); // cand *= den
1942      const poly t = singclap_pdivide( cand, gcd, R ); // cand' * den / gcd(cand', den)
1943      p_Delete(&cand, R);
1944      p_Delete(&gcd, R);
1945      cand = t;
1946    }
1947  }
1948  while( numberCollectionEnumerator.MoveNext() );
1949
1950  if( cand == NULL )
1951  {
1952    c = ntInit(1, cf);
1953    return;
1954  }
1955
1956  c = ntInit(cand, cf);
1957
1958  numberCollectionEnumerator.Reset();
1959
1960  number d = NULL;
1961
1962  while (numberCollectionEnumerator.MoveNext() )
1963  {
1964    number &n = numberCollectionEnumerator.Current();
1965    number t = ntMult(n, c, cf); // TODO: rewrite!?
1966    ntDelete(&n, cf);
1967
1968    ntNormalize(t, cf); // TODO: needed?
1969    n = t;
1970
1971    fraction f = (fraction)t;
1972    assume( f != NULL );
1973
1974    const poly den = DEN(f);
1975
1976    if( den != NULL ) // ?? / ?? ?
1977    {
1978      assume( p_IsConstant(den, R) );
1979      assume( pNext(den) == NULL );
1980
1981      if( d == NULL )
1982        d = n_Copy(pGetCoeff(den), Q);
1983      else
1984      {
1985        number g = n_Lcm(d, pGetCoeff(den), Q);
1986        n_Delete(&d, Q); d = g;
1987      }
1988    }
1989  }
1990
1991  if( d != NULL )
1992  {
1993    numberCollectionEnumerator.Reset();
1994    while (numberCollectionEnumerator.MoveNext() )
1995    {
1996      number &n = numberCollectionEnumerator.Current();
1997      fraction f = (fraction)n;
1998
1999      assume( f != NULL );
2000
2001      const poly den = DEN(f);
2002
2003      if( den == NULL ) // ?? / 1 ?
2004        NUM(f) = p_Mult_nn(NUM(f), d, R);
2005      else
2006      {
2007        assume( p_IsConstant(den, R) );
2008        assume( pNext(den) == NULL );
2009
2010        number ddd = n_Div(d, pGetCoeff(den), Q); // but be an integer now!!!
2011        NUM(f) = p_Mult_nn(NUM(f), ddd, R);
2012        n_Delete(&ddd, Q);
2013
2014        p_Delete(&DEN(f), R);
2015        DEN(f) = NULL; // TODO: check if this is needed!?
2016      }
2017
2018      assume( DEN(f) == NULL );
2019    }
2020
2021    NUM(c) = p_Mult_nn(NUM(c), d, R);
2022    n_Delete(&d, Q);
2023  }
2024
2025
2026  ntTest(c);
2027}
2028
2029BOOLEAN ntInitChar(coeffs cf, void * infoStruct)
2030{
2031
2032  assume( infoStruct != NULL );
2033
2034  TransExtInfo *e = (TransExtInfo *)infoStruct;
2035
2036  assume( e->r                != NULL);      // extRing;
2037  assume( e->r->cf            != NULL);      // extRing->cf;
2038  assume( e->r->qideal == NULL );
2039
2040  assume( cf != NULL );
2041  assume(getCoeffType(cf) == ID);                // coeff type;
2042
2043  ring R = e->r;
2044  assume(R != NULL);
2045
2046  R->ref ++; // increase the ref.counter for the ground poly. ring!
2047
2048  cf->extRing           = R;
2049  /* propagate characteristic up so that it becomes
2050     directly accessible in cf: */
2051  cf->ch = R->cf->ch;
2052  cf->factoryVarOffset = R->cf->factoryVarOffset + rVar(R);
2053
2054  cf->cfGreaterZero  = ntGreaterZero;
2055  cf->cfGreater      = ntGreater;
2056  cf->cfEqual        = ntEqual;
2057  cf->cfIsZero       = ntIsZero;
2058  cf->cfIsOne        = ntIsOne;
2059  cf->cfIsMOne       = ntIsMOne;
2060  cf->cfInit         = ntInit;
2061  cf->cfInit_bigint  = ntInit_bigint;
2062  cf->cfInt          = ntInt;
2063  cf->cfNeg          = ntNeg;
2064  cf->cfAdd          = ntAdd;
2065  cf->cfSub          = ntSub;
2066  cf->cfMult         = ntMult;
2067  cf->cfDiv          = ntDiv;
2068  cf->cfExactDiv     = ntDiv;
2069  cf->cfPower        = ntPower;
2070  cf->cfCopy         = ntCopy;
2071  cf->cfWriteLong    = ntWriteLong;
2072  cf->cfRead         = ntRead;
2073  cf->cfNormalize    = ntNormalize;
2074  cf->cfDelete       = ntDelete;
2075  cf->cfSetMap       = ntSetMap;
2076  cf->cfGetDenom     = ntGetDenom;
2077  cf->cfGetNumerator = ntGetNumerator;
2078  cf->cfRePart       = ntCopy;
2079  cf->cfImPart       = ntImPart;
2080  cf->cfCoeffWrite   = ntCoeffWrite;
2081#ifdef LDEBUG
2082  cf->cfDBTest       = ntDBTest;
2083#endif
2084  cf->cfGcd          = ntGcd;
2085  cf->cfLcm          = ntLcm;
2086  cf->cfSize         = ntSize;
2087  cf->nCoeffIsEqual  = ntCoeffIsEqual;
2088  cf->cfInvers       = ntInvers;
2089  cf->cfIntDiv       = ntDiv;
2090  cf->cfKillChar     = ntKillChar;
2091
2092  if( rCanShortOut(ntRing) )
2093    cf->cfWriteShort = ntWriteShort;
2094  else
2095    cf->cfWriteShort = ntWriteLong;
2096
2097#ifndef HAVE_FACTORY
2098  PrintS("// Warning: The 'factory' module is not available.\n");
2099  PrintS("//          Hence gcd's cannot be cancelled in any\n");
2100  PrintS("//          computed fraction!\n");
2101#else
2102  cf->convFactoryNSingN =ntConvFactoryNSingN;
2103  cf->convSingNFactoryN =ntConvSingNFactoryN;
2104#endif
2105  cf->cfParDeg = ntParDeg;
2106
2107  cf->iNumberOfParameters = rVar(R);
2108  cf->pParameterNames = R->names;
2109  cf->cfParameter = ntParameter;
2110
2111  if( nCoeff_is_Q(R->cf) )
2112    cf->cfClearContent = ntClearContent;
2113
2114  cf->cfClearDenominators = ntClearDenominators;
2115
2116  return FALSE;
2117}
Note: See TracBrowser for help on using the repository browser.