source: git/libpolys/polys/ext_fields/transext.cc @ 5a4e17

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