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

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