source: git/libpolys/polys/ext_fields/transext.cc @ 6a1aa7

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