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

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