source: git/libpolys/polys/ext_fields/transext.cc

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