source: git/libpolys/polys/ext_fields/transext.cc @ 0635d51

spielwiese
Last change on this file since 0635d51 was 0635d51, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: missing content over Z in ntLcm, ntGcd
  • Property mode set to 100644
File size: 60.3 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  poly pGcd;
1403  if (nCoeff_is_Q(ntCoeffs))
1404  {
1405    if (p_IsConstant(pa,ntRing) && p_IsConstant(pb,ntRing))
1406    {
1407      pGcd = pa;
1408      p_SetCoeff (pGcd, n_Gcd (pGetCoeff(pGcd), pGetCoeff(pb), ntCoeffs), ntRing);
1409    }
1410    else
1411    {
1412      number contentpa, contentpb, tmp;
1413
1414      contentpb= p_GetCoeff(pb, ntRing);
1415      pIter(pb);
1416      while (pb != NULL)
1417      {
1418        tmp = n_Gcd(contentpb, p_GetCoeff(pb, ntRing) , ntCoeffs);
1419        n_Delete(&contentpb, ntCoeffs);
1420        contentpb = tmp;
1421        pIter(pb);
1422      }
1423
1424      contentpa= p_GetCoeff(pa, ntRing);
1425      pIter(pa);
1426      while (pa != NULL)
1427      {
1428        tmp = n_Gcd(contentpa, p_GetCoeff(pa, ntRing), ntCoeffs);
1429        n_Delete(&contentpa, ntCoeffs);
1430        contentpa = tmp;
1431        pIter(pa);
1432      }
1433
1434      tmp= n_Gcd (contentpb, contentpa, ntCoeffs);
1435      n_Delete(&contentpa, ntCoeffs);
1436      n_Delete(&contentpb, ntCoeffs);
1437      contentpa= tmp;
1438      p_Delete(&pb, ntRing);
1439      p_Delete(&pa, ntRing);
1440
1441      /* singclap_gcd destroys its arguments; we hence need copies: */
1442      pGcd = singclap_gcd(p_Copy(NUM(fa),ntRing), p_Copy(DEN(fb),ntRing), ntRing);
1443      pGcd= p_Mult_nn (pGcd, contentpa, ntRing);
1444      n_Delete(&contentpa, ntCoeffs);
1445    }
1446  }
1447  else
1448    pGcd = singclap_gcd(pa, pb, cf->extRing);
1449
1450  /* Note that, over Q, singclap_gcd will remove the denominators in all
1451     rational coefficients of pa and pb, before starting to compute
1452     the gcd. Thus, we do not need to ensure that the coefficients of
1453     pa and pb live in Z; they may well be elements of Q\Z. */
1454
1455  if (p_IsConstant(pGcd, ntRing) &&
1456      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
1457  { /* gcd = 1; return pa*pb*/
1458    p_Delete(&pGcd,ntRing);
1459    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1460    NUM(result) = pp_Mult_qq(NUM(fa),DEN(fb),ntRing);
1461
1462    ntTest((number)result); // !!!!
1463
1464    return (number)result;
1465  }
1466
1467
1468  /* return pa*pb/gcd */
1469    poly newNum = singclap_pdivide(NUM(fa), pGcd, ntRing);
1470    p_Delete(&pGcd,ntRing);
1471    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1472    NUM(result) = p_Mult_q(p_Copy(DEN(fb),ntRing),newNum,ntRing);
1473    ntTest((number)result); // !!!!
1474    return (number)result;
1475
1476#else
1477  Print("// factory needed: transext.cc:ntLcm\n");
1478  return NULL;
1479#endif /* HAVE_FACTORY */
1480  return NULL;
1481}
1482
1483number ntGcd(number a, number b, const coeffs cf)
1484{
1485  ntTest(a);
1486  ntTest(b);
1487  if (a==NULL) return ntCopy(b,cf);
1488  if (b==NULL) return ntCopy(a,cf);
1489#ifdef HAVE_FACTORY
1490  fraction fa = (fraction)a;
1491  fraction fb = (fraction)b;
1492
1493  poly pa = p_Copy(NUM(fa), ntRing);
1494  poly pb = p_Copy(NUM(fb), ntRing);
1495
1496  poly pGcd;
1497  if (nCoeff_is_Q(ntCoeffs))
1498  {
1499    if (p_IsConstant(pa,ntRing) && p_IsConstant(pb,ntRing))
1500    {
1501      pGcd = pa;
1502      p_SetCoeff (pGcd, n_Gcd (pGetCoeff(pGcd), pGetCoeff(pb), ntCoeffs), ntRing);
1503    }
1504    else
1505    {
1506      number contentpa, contentpb, tmp;
1507
1508      contentpb= p_GetCoeff(pb, ntRing);
1509      pIter(pb);
1510      while (pb != NULL)
1511      {
1512        tmp = n_Gcd(contentpb, p_GetCoeff(pb, ntRing) , ntCoeffs);
1513        n_Delete(&contentpb, ntCoeffs);
1514        contentpb = tmp;
1515        pIter(pb);
1516      }
1517
1518      contentpa= p_GetCoeff(pa, ntRing);
1519      pIter(pa);
1520      while (pa != NULL)
1521      {
1522        tmp = n_Gcd(contentpa, p_GetCoeff(pa, ntRing), ntCoeffs);
1523        n_Delete(&contentpa, ntCoeffs);
1524        contentpa = tmp;
1525        pIter(pa);
1526      }
1527
1528      tmp= n_Gcd (contentpb, contentpa, ntCoeffs);
1529      n_Delete(&contentpa, ntCoeffs);
1530      n_Delete(&contentpb, ntCoeffs);
1531      contentpa= tmp;
1532      p_Delete(&pb, ntRing);
1533      p_Delete(&pa, ntRing);
1534
1535      /* singclap_gcd destroys its arguments; we hence need copies: */
1536      pGcd = singclap_gcd(p_Copy(NUM(fa),ntRing), p_Copy(NUM(fb),ntRing), ntRing);
1537      pGcd= p_Mult_nn (pGcd, contentpa, ntRing);
1538      n_Delete(&contentpa, ntCoeffs);
1539    }
1540  }
1541  else
1542    pGcd = singclap_gcd(pa, pb, cf->extRing);
1543  /* Note that, over Q, singclap_gcd will remove the denominators in all
1544     rational coefficients of pa and pb, before starting to compute
1545     the gcd. Thus, we do not need to ensure that the coefficients of
1546     pa and pb live in Z; they may well be elements of Q\Z. */
1547
1548  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1549  NUM(result) = pGcd;
1550  ntTest((number)result); // !!!!
1551  return (number)result;
1552#else
1553  Print("// factory needed: transext.cc:ntGcd\n");
1554  return NULL;
1555#endif /* HAVE_FACTORY */
1556}
1557
1558int ntSize(number a, const coeffs cf)
1559{
1560  ntTest(a);
1561  if (IS0(a)) return -1;
1562  /* this has been taken from the old implementation of field extensions,
1563     where we computed the sum of the degrees and the numbers of terms in
1564     the numerator and denominator of a; so we leave it at that, for the
1565     time being */
1566  fraction f = (fraction)a;
1567  poly p = NUM(f);
1568  int noOfTerms = 0;
1569  int numDegree = 0;
1570  while (p != NULL)
1571  {
1572    noOfTerms++;
1573    int d = 0;
1574    for (int i = 1; i <= rVar(ntRing); i++)
1575      d += p_GetExp(p, i, ntRing);
1576    if (d > numDegree) numDegree = d;
1577    pIter(p);
1578  }
1579  int denDegree = 0;
1580  if (!DENIS1(f))
1581  {
1582    p = DEN(f);
1583    while (p != NULL)
1584    {
1585      noOfTerms++;
1586      int d = 0;
1587      for (int i = 1; i <= rVar(ntRing); i++)
1588        d += p_GetExp(p, i, ntRing);
1589      if (d > denDegree) denDegree = d;
1590      pIter(p);
1591    }
1592  }
1593  ntTest(a); // !!!!
1594  return numDegree + denDegree + noOfTerms;
1595}
1596
1597number ntInvers(number a, const coeffs cf)
1598{
1599  ntTest(a);
1600  if (IS0(a))
1601  {
1602    WerrorS(nDivBy0);
1603    return NULL;
1604  }
1605  fraction f = (fraction)a;
1606  assume( f != NULL );
1607
1608  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1609
1610  assume( NUM(f) != NULL );
1611  const poly den = DEN(f);
1612
1613  if (den == NULL)
1614    NUM(result) = p_One(ntRing);
1615  else
1616    NUM(result) = p_Copy(den, ntRing);
1617
1618  if( !NUMIS1(f) )
1619  {
1620    poly num_f=NUM(f);
1621    BOOLEAN neg= !n_GreaterZero(pGetCoeff(num_f),ntCoeffs);
1622    if (neg)
1623    {
1624      num_f=p_Neg(p_Copy(num_f, ntRing), ntRing);
1625      NUM(result)=p_Neg(NUM(result), ntRing);
1626    }
1627    else
1628    {
1629      num_f=p_Copy(num_f, ntRing);
1630    }
1631    DEN(result) = num_f;
1632    COM(result) = COM(f);
1633    if (neg)
1634    {
1635      if (p_IsOne(num_f, ntRing))
1636      {
1637        DEN(result)=NULL;
1638        COM(result) = 0;
1639        p_Delete(&num_f,ntRing);
1640      }
1641    }
1642  }
1643  else
1644  {
1645    DEN(result) = NULL;
1646    COM(result) = 0;
1647  }
1648  ntTest((number)result); // !!!!
1649  return (number)result;
1650}
1651
1652/* assumes that src = Q, dst = Q(t_1, ..., t_s) */
1653number ntMap00(number a, const coeffs src, const coeffs dst)
1654{
1655  if (n_IsZero(a, src)) return NULL;
1656  assume(n_Test(a, src));
1657  assume(src == dst->extRing->cf);
1658  if ((SR_HDL(a) & SR_INT) || (a->s==3))
1659  {
1660    number res=ntInit(p_NSet(n_Copy(a, src), dst->extRing), dst);
1661    n_Test(res,dst);
1662    return res;
1663  }
1664  number nn=nlGetDenom(a,src);
1665  number zz=nlGetNumerator(a,src);
1666  number res=ntInit(p_NSet(zz,dst->extRing), dst);
1667  fraction ff=(fraction)res;
1668  if (n_IsOne(nn,src)) DEN(ff)=NULL;
1669  else                 DEN(ff)=p_NSet(nn,dst->extRing);
1670  n_Test((number)ff,dst);
1671  return (number)ff;
1672}
1673
1674/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
1675number ntMapP0(number a, const coeffs src, const coeffs dst)
1676{
1677  if (n_IsZero(a, src)) return NULL;
1678  assume(n_Test(a, src));
1679  /* mapping via intermediate int: */
1680  int n = n_Int(a, src);
1681  number q = n_Init(n, dst->extRing->cf);
1682  if (n_IsZero(q, dst->extRing->cf))
1683  {
1684    n_Delete(&q, dst->extRing->cf);
1685    return NULL;
1686  }
1687  return ntInit(p_NSet(q, dst->extRing), dst);
1688}
1689
1690/* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
1691                       src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
1692number ntCopyMap(number a, const coeffs cf, const coeffs dst)
1693{
1694//  if (n_IsZero(a, cf)) return NULL;
1695
1696  ntTest(a);
1697
1698  if (IS0(a)) return NULL;
1699
1700  const ring rSrc = cf->extRing;
1701  const ring rDst = dst->extRing;
1702
1703  if( rSrc == rDst )
1704    return ntCopy(a, dst); // USUALLY WRONG!
1705
1706  fraction f = (fraction)a;
1707  poly g = prCopyR(NUM(f), rSrc, rDst);
1708
1709  poly h = NULL;
1710
1711  if (!DENIS1(f))
1712     h = prCopyR(DEN(f), rSrc, rDst);
1713
1714  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1715
1716  NUM(result) = g;
1717  DEN(result) = h;
1718  COM(result) = COM(f);
1719  assume(n_Test((number)result, dst));
1720  return (number)result;
1721}
1722
1723number ntCopyAlg(number a, const coeffs cf, const coeffs dst)
1724{
1725  assume( n_Test(a, cf) );
1726  if (n_IsZero(a, cf)) return NULL;
1727
1728  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1729  // DEN(f) = NULL; COM(f) = 0;
1730  NUM(f) = prCopyR((poly)a, cf->extRing, dst->extRing);
1731  assume(n_Test((number)f, dst));
1732  return (number)f;
1733}
1734
1735/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
1736number ntMap0P(number a, const coeffs src, const coeffs dst)
1737{
1738  assume( n_Test(a, src) );
1739  if (n_IsZero(a, src)) return NULL;
1740  // int p = rChar(dst->extRing);
1741  number q = nlModP(a, src, dst->extRing->cf);
1742
1743  if (n_IsZero(q, dst->extRing->cf))
1744  {
1745    n_Delete(&q, dst->extRing->cf);
1746    return NULL;
1747  }
1748
1749  poly g = p_NSet(q, dst->extRing);
1750
1751  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1752  NUM(f) = g; // DEN(f) = NULL; COM(f) = 0;
1753  assume(n_Test((number)f, dst));
1754  return (number)f;
1755}
1756
1757/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
1758number ntMapPP(number a, const coeffs src, const coeffs dst)
1759{
1760  assume( n_Test(a, src) );
1761  if (n_IsZero(a, src)) return NULL;
1762  assume(src == dst->extRing->cf);
1763  poly p = p_One(dst->extRing);
1764  p_SetCoeff(p, n_Copy(a, src), dst->extRing);
1765  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1766  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1767  assume(n_Test((number)f, dst));
1768  return (number)f;
1769}
1770
1771/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
1772number ntMapUP(number a, const coeffs src, const coeffs dst)
1773{
1774  assume( n_Test(a, src) );
1775  if (n_IsZero(a, src)) return NULL;
1776  /* mapping via intermediate int: */
1777  int n = n_Int(a, src);
1778  number q = n_Init(n, dst->extRing->cf);
1779  poly p;
1780  if (n_IsZero(q, dst->extRing->cf))
1781  {
1782    n_Delete(&q, dst->extRing->cf);
1783    return NULL;
1784  }
1785  p = p_One(dst->extRing);
1786  p_SetCoeff(p, q, dst->extRing);
1787  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1788  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1789  assume(n_Test((number)f, dst));
1790  return (number)f;
1791}
1792
1793nMapFunc ntSetMap(const coeffs src, const coeffs dst)
1794{
1795  /* dst is expected to be a rational function field */
1796  assume(getCoeffType(dst) == ID);
1797
1798  if( src == dst ) return ndCopyMap;
1799
1800  int h = 0; /* the height of the extension tower given by dst */
1801  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
1802  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
1803
1804  /* for the time being, we only provide maps if h = 1 and if b is Q or
1805     some field Z/pZ: */
1806  if (h==0)
1807  {
1808    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
1809      return ntMap00;                                 /// Q       -->  Q(T)
1810    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
1811      return ntMapP0;                                 /// Z/p     -->  Q(T)
1812    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
1813      return ntMap0P;                                 /// Q       --> Z/p(T)
1814    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
1815    {
1816      if (src->ch == dst->ch) return ntMapPP;         /// Z/p     --> Z/p(T)
1817      else return ntMapUP;                            /// Z/u     --> Z/p(T)
1818    }
1819  }
1820  if (h != 1) return NULL;
1821  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
1822
1823  /* Let T denote the sequence of transcendental extension variables, i.e.,
1824     K[t_1, ..., t_s] =: K[T];
1825     Let moreover, for any such sequence T, T' denote any subsequence of T
1826     of the form t_1, ..., t_w with w <= s. */
1827
1828  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
1829
1830  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
1831  {
1832    if (rVar(src->extRing) > rVar(dst->extRing))
1833       return NULL;
1834
1835    for (int i = 0; i < rVar(src->extRing); i++)
1836      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0)
1837         return NULL;
1838
1839    if (src->type==n_transExt)
1840       return ntCopyMap;          /// Q(T')   --> Q(T)
1841    else
1842       return ntCopyAlg;
1843  }
1844
1845  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
1846  {
1847    if (rVar(src->extRing) > rVar(dst->extRing))
1848       return NULL;
1849
1850    for (int i = 0; i < rVar(src->extRing); i++)
1851      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0)
1852         return NULL;
1853
1854    if (src->type==n_transExt)
1855       return ntCopyMap;         /// Z/p(T') --> Z/p(T)
1856    else
1857       return ntCopyAlg;
1858  }
1859
1860  return NULL;                                 /// default
1861}
1862#if 0
1863nMapFunc ntSetMap_T(const coeffs src, const coeffs dst)
1864{
1865  nMapFunc n=ntSetMap(src,dst);
1866  if (n==ntCopyAlg) printf("n=ntCopyAlg\n");
1867  else if (n==ntCopyMap) printf("n=ntCopyMap\n");
1868  else if (n==ntMapUP) printf("n=ntMapUP\n");
1869  else if (n==ntMap0P) printf("n=ntMap0P\n");
1870  else if (n==ntMapP0) printf("n=ntMapP0\n");
1871  else if (n==ntMap00) printf("n=ntMap00\n");
1872  else if (n==NULL) printf("n=NULL\n");
1873  else printf("n=?\n");
1874  return n;
1875}
1876#endif
1877
1878void ntKillChar(coeffs cf)
1879{
1880  if ((--cf->extRing->ref) == 0)
1881    rDelete(cf->extRing);
1882}
1883#ifdef HAVE_FACTORY
1884number ntConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
1885{
1886  if (n.isZero()) return NULL;
1887  poly p=convFactoryPSingP(n,ntRing);
1888  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1889  NUM(result) = p;
1890  //DEN(result) = NULL; // done by omAlloc0Bin
1891  //COM(result) = 0; // done by omAlloc0Bin
1892  ntTest((number)result);
1893  return (number)result;
1894}
1895CanonicalForm ntConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
1896{
1897  ntTest(n);
1898  if (IS0(n)) return CanonicalForm(0);
1899
1900  fraction f = (fraction)n;
1901  return convSingPFactoryP(NUM(f),ntRing);
1902}
1903#endif
1904
1905static int ntParDeg(number a, const coeffs cf)
1906{
1907  ntTest(a);
1908  if (IS0(a)) return -1;
1909  fraction fa = (fraction)a;
1910  return cf->extRing->pFDeg(NUM(fa),cf->extRing);
1911}
1912
1913/// return the specified parameter as a number in the given trans.ext.
1914static number ntParameter(const int iParameter, const coeffs cf)
1915{
1916  assume(getCoeffType(cf) == ID);
1917
1918  const ring R = cf->extRing;
1919  assume( R != NULL );
1920  assume( 0 < iParameter && iParameter <= rVar(R) );
1921
1922  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
1923  p_Test(p,R);
1924
1925  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1926  NUM(f) = p;
1927  DEN(f) = NULL;
1928  COM(f) = 0;
1929
1930  ntTest((number)f);
1931
1932  return (number)f;
1933}
1934
1935/// if m == var(i)/1 => return i,
1936int ntIsParam(number m, const coeffs cf)
1937{
1938  ntTest(m);
1939  assume(getCoeffType(cf) == ID);
1940
1941  const ring R = cf->extRing;
1942  assume( R != NULL );
1943
1944  fraction f = (fraction)m;
1945
1946  if( DEN(f) != NULL )
1947    return 0;
1948
1949  return p_Var( NUM(f), R );
1950}
1951
1952struct NTNumConverter
1953{
1954  static inline poly convert(const number& n)
1955  {
1956    // suitable for trans. ext. numbers that are fractions of polys
1957    return NUM((fraction)n); // return the numerator
1958  }
1959};
1960
1961
1962static void ntClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1963{
1964  assume(cf != NULL);
1965  assume(getCoeffType(cf) == ID);
1966  // all coeffs are given by fractions of polynomails over integers!!!
1967  // without denominators!!!
1968
1969  const ring   R = cf->extRing;
1970  assume(R != NULL);
1971  const coeffs Q = R->cf;
1972  assume(Q != NULL);
1973  assume(nCoeff_is_Q(Q));
1974
1975
1976  numberCollectionEnumerator.Reset();
1977
1978  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
1979  {
1980    c = ntInit(1, cf);
1981    return;
1982  }
1983
1984  // all coeffs are given by integers after returning from this routine
1985
1986  // part 1, collect product of all denominators /gcds
1987  poly cand = NULL;
1988
1989  do
1990  {
1991    number &n = numberCollectionEnumerator.Current();
1992
1993    ntNormalize(n, cf);
1994
1995    fraction f = (fraction)n;
1996
1997    assume( f != NULL );
1998
1999    const poly den = DEN(f);
2000
2001    assume( den == NULL ); // ?? / 1 ?
2002
2003    const poly num = NUM(f);
2004
2005    if( cand == NULL )
2006      cand = p_Copy(num, R);
2007    else
2008      cand = singclap_gcd(cand, p_Copy(num, R), R); // gcd(cand, num)
2009
2010    if( p_IsConstant(cand, R) )
2011      break;
2012  }
2013  while( numberCollectionEnumerator.MoveNext() ) ;
2014
2015
2016  // part2: all coeffs = all coeffs * cand
2017  if( cand != NULL )
2018  {
2019  if( !p_IsConstant(cand, R) )
2020  {
2021    c = ntInit(cand, cf);
2022    numberCollectionEnumerator.Reset();
2023    while (numberCollectionEnumerator.MoveNext() )
2024    {
2025      number &n = numberCollectionEnumerator.Current();
2026      const number t = ntDiv(n, c, cf); // TODO: rewrite!?
2027      ntDelete(&n, cf);
2028      n = t;
2029    }
2030  } // else NUM (result) = p_One(R);
2031  else { p_Delete(&cand, R); cand = NULL; }
2032  }
2033
2034  // Quick and dirty fix for constant content clearing: consider numerators???
2035  CRecursivePolyCoeffsEnumerator<NTNumConverter> itr(numberCollectionEnumerator); // recursively treat the NUM(numbers) as polys!
2036  number cc;
2037
2038  n_ClearContent(itr, cc, Q);
2039  number g = ntInit(p_NSet(cc, R), cf);
2040
2041  if( cand != NULL )
2042  {
2043    number gg = ntMult(g, c, cf);
2044    ntDelete(&g, cf);
2045    ntDelete(&c, cf); c = gg;
2046  } else
2047    c = g;
2048  ntTest(c);
2049}
2050
2051static void ntClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
2052{
2053  assume(cf != NULL);
2054  assume(getCoeffType(cf) == ID); // both over Q(a) and Zp(a)!
2055  // all coeffs are given by fractions of polynomails over integers!!!
2056
2057  numberCollectionEnumerator.Reset();
2058
2059  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2060  {
2061    c = ntInit(1, cf);
2062    return;
2063  }
2064
2065  // all coeffs are given by integers after returning from this routine
2066
2067  // part 1, collect product of all denominators /gcds
2068  poly cand = NULL;
2069
2070  const ring R = cf->extRing;
2071  assume(R != NULL);
2072
2073  const coeffs Q = R->cf;
2074  assume(Q != NULL);
2075//  assume(nCoeff_is_Q(Q));
2076
2077  do
2078  {
2079    number &n = numberCollectionEnumerator.Current();
2080
2081    ntNormalize(n, cf);
2082
2083    fraction f = (fraction)ntGetDenom (n, cf);
2084
2085    assume( f != NULL );
2086
2087    const poly den = NUM(f);
2088
2089    if( den == NULL ) // ?? / 1 ?
2090      continue;
2091
2092    if( cand == NULL )
2093      cand = p_Copy(den, R);
2094    else
2095    {
2096      // cand === LCM( cand, den )!!!!
2097      // NOTE: maybe it's better to make the product and clearcontent afterwards!?
2098      // TODO: move the following to factory?
2099      poly gcd = singclap_gcd(p_Copy(cand, R), p_Copy(den, R), R); // gcd(cand, den) is monic no mater leading coeffs! :((((
2100      if (nCoeff_is_Q (Q))
2101      {
2102        number LcGcd= n_Gcd (p_GetCoeff (cand, R), p_GetCoeff(den, R), Q);
2103        gcd = p_Mult_nn(gcd, LcGcd, R);
2104        n_Delete(&LcGcd,Q);
2105      }
2106//      assume( n_IsOne(pGetCoeff(gcd), Q) ); // TODO: this may be wrong...
2107      cand = p_Mult_q(cand, p_Copy(den, R), R); // cand *= den
2108      const poly t = singclap_pdivide( cand, gcd, R ); // cand' * den / gcd(cand', den)
2109      p_Delete(&cand, R);
2110      p_Delete(&gcd, R);
2111      cand = t;
2112    }
2113  }
2114  while( numberCollectionEnumerator.MoveNext() );
2115
2116  if( cand == NULL )
2117  {
2118    c = ntInit(1, cf);
2119    return;
2120  }
2121
2122  c = ntInit(cand, cf);
2123
2124  numberCollectionEnumerator.Reset();
2125
2126  number d = NULL;
2127
2128  while (numberCollectionEnumerator.MoveNext() )
2129  {
2130    number &n = numberCollectionEnumerator.Current();
2131    number t = ntMult(n, c, cf); // TODO: rewrite!?
2132    ntDelete(&n, cf);
2133
2134    ntNormalize(t, cf); // TODO: needed?
2135    n = t;
2136
2137    fraction f = (fraction)t;
2138    assume( f != NULL );
2139
2140    const poly den = DEN(f);
2141
2142    if( den != NULL ) // ?? / ?? ?
2143    {
2144      assume( p_IsConstant(den, R) );
2145      assume( pNext(den) == NULL );
2146
2147      if( d == NULL )
2148        d = n_Copy(pGetCoeff(den), Q);
2149      else
2150      {
2151        number g = n_Lcm(d, pGetCoeff(den), Q);
2152        n_Delete(&d, Q); d = g;
2153      }
2154    }
2155  }
2156
2157  if( d != NULL )
2158  {
2159    numberCollectionEnumerator.Reset();
2160    while (numberCollectionEnumerator.MoveNext() )
2161    {
2162      number &n = numberCollectionEnumerator.Current();
2163      fraction f = (fraction)n;
2164
2165      assume( f != NULL );
2166
2167      const poly den = DEN(f);
2168
2169      if( den == NULL ) // ?? / 1 ?
2170        NUM(f) = p_Mult_nn(NUM(f), d, R);
2171      else
2172      {
2173        assume( p_IsConstant(den, R) );
2174        assume( pNext(den) == NULL );
2175
2176        number ddd = n_Div(d, pGetCoeff(den), Q); // but be an integer now!!!
2177        NUM(f) = p_Mult_nn(NUM(f), ddd, R);
2178        n_Delete(&ddd, Q);
2179
2180        p_Delete(&DEN(f), R);
2181        DEN(f) = NULL; // TODO: check if this is needed!?
2182      }
2183
2184      assume( DEN(f) == NULL );
2185    }
2186
2187    NUM(c) = p_Mult_nn(NUM(c), d, R);
2188    n_Delete(&d, Q);
2189  }
2190
2191
2192  ntTest(c);
2193}
2194
2195BOOLEAN ntInitChar(coeffs cf, void * infoStruct)
2196{
2197
2198  assume( infoStruct != NULL );
2199
2200  TransExtInfo *e = (TransExtInfo *)infoStruct;
2201
2202  assume( e->r                != NULL);      // extRing;
2203  assume( e->r->cf            != NULL);      // extRing->cf;
2204  assume( e->r->qideal == NULL );
2205
2206  assume( cf != NULL );
2207  assume(getCoeffType(cf) == ID);                // coeff type;
2208
2209  ring R = e->r;
2210  assume(R != NULL);
2211
2212  R->ref ++; // increase the ref.counter for the ground poly. ring!
2213
2214  cf->extRing           = R;
2215  /* propagate characteristic up so that it becomes
2216     directly accessible in cf: */
2217  cf->ch = R->cf->ch;
2218  cf->factoryVarOffset = R->cf->factoryVarOffset + rVar(R);
2219
2220  cf->cfGreaterZero  = ntGreaterZero;
2221  cf->cfGreater      = ntGreater;
2222  cf->cfEqual        = ntEqual;
2223  cf->cfIsZero       = ntIsZero;
2224  cf->cfIsOne        = ntIsOne;
2225  cf->cfIsMOne       = ntIsMOne;
2226  cf->cfInit         = ntInit;
2227  cf->cfInit_bigint  = ntInit_bigint;
2228  cf->cfInt          = ntInt;
2229  cf->cfNeg          = ntNeg;
2230  cf->cfAdd          = ntAdd;
2231  cf->cfSub          = ntSub;
2232  cf->cfMult         = ntMult;
2233  cf->cfDiv          = ntDiv;
2234  cf->cfExactDiv     = ntDiv;
2235  cf->cfPower        = ntPower;
2236  cf->cfCopy         = ntCopy;
2237  cf->cfWriteLong    = ntWriteLong;
2238  cf->cfRead         = ntRead;
2239  cf->cfNormalize    = ntNormalize;
2240  cf->cfDelete       = ntDelete;
2241  cf->cfSetMap       = ntSetMap;
2242  cf->cfGetDenom     = ntGetDenom;
2243  cf->cfGetNumerator = ntGetNumerator;
2244  cf->cfRePart       = ntCopy;
2245  cf->cfImPart       = ntImPart;
2246  cf->cfCoeffWrite   = ntCoeffWrite;
2247#ifdef LDEBUG
2248  cf->cfDBTest       = ntDBTest;
2249#endif
2250  cf->cfGcd          = ntGcd;
2251  cf->cfLcm          = ntLcm;
2252  cf->cfSize         = ntSize;
2253  cf->nCoeffIsEqual  = ntCoeffIsEqual;
2254  cf->cfInvers       = ntInvers;
2255  cf->cfIntDiv       = ntDiv;
2256  cf->cfKillChar     = ntKillChar;
2257
2258  if( rCanShortOut(ntRing) )
2259    cf->cfWriteShort = ntWriteShort;
2260  else
2261    cf->cfWriteShort = ntWriteLong;
2262
2263#ifndef HAVE_FACTORY
2264  PrintS("// Warning: The 'factory' module is not available.\n");
2265  PrintS("//          Hence gcd's cannot be cancelled in any\n");
2266  PrintS("//          computed fraction!\n");
2267#else
2268  cf->convFactoryNSingN =ntConvFactoryNSingN;
2269  cf->convSingNFactoryN =ntConvSingNFactoryN;
2270#endif
2271  cf->cfParDeg = ntParDeg;
2272
2273  cf->iNumberOfParameters = rVar(R);
2274  cf->pParameterNames = R->names;
2275  cf->cfParameter = ntParameter;
2276
2277  if( nCoeff_is_Q(R->cf) )
2278    cf->cfClearContent = ntClearContent;
2279
2280  cf->cfClearDenominators = ntClearDenominators;
2281
2282  return FALSE;
2283}
2284
2285template class CRecursivePolyCoeffsEnumerator<NTNumConverter>;
2286template class IEnumerator<snumber*>;
Note: See TracBrowser for help on using the repository browser.