source: git/libpolys/polys/ext_fields/transext.cc @ 161e20

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