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

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