source: git/libpolys/polys/ext_fields/transext.cc @ 2e4ec14

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