source: git/libpolys/polys/ext_fields/transext.cc @ 0c516e

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