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

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