source: git/libpolys/polys/ext_fields/transext.cc @ 28216f

spielwiese
Last change on this file since 28216f was 28216f, checked in by Martin Lee <martinlee84@…>, 12 years ago
fix: ntGetDenom & ntGetNumerator
  • Property mode set to 100644
File size: 42.4 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#define TRANSEXT_PRIVATES
34
35#include "config.h"
36#include <misc/auxiliary.h>
37
38#include <omalloc/omalloc.h>
39
40#include <reporter/reporter.h>
41
42#include <coeffs/coeffs.h>
43#include <coeffs/numbers.h>
44#include <coeffs/longrat.h>
45
46#include <polys/monomials/ring.h>
47#include <polys/monomials/p_polys.h>
48#include <polys/simpleideals.h>
49
50#ifdef HAVE_FACTORY
51#include <polys/clapsing.h>
52#include <polys/clapconv.h>
53#include <factory/factory.h>
54#endif
55
56#include "ext_fields/transext.h"
57#include "prCopy.h"
58
59/* constants for controlling the complexity of numbers */
60#define ADD_COMPLEXITY 1   /**< complexity increase due to + and - */
61#define MULT_COMPLEXITY 2   /**< complexity increase due to * and / */
62#define BOUND_COMPLEXITY 10   /**< maximum complexity of a number */
63
64#define NUMIS1(f) (p_IsConstant(NUM(f), cf->extRing) && \
65                   n_IsOne(p_GetCoeff(NUM(f), cf->extRing), \
66                           cf->extRing->cf))
67                   /**< TRUE iff num. represents 1 */
68#define COM(f) f->complexity
69
70
71#ifdef LDEBUG
72#define ntTest(a) ntDBTest(a,__FILE__,__LINE__,cf)
73BOOLEAN  ntDBTest(number a, const char *f, const int l, const coeffs r);
74#else
75#define ntTest(a) (TRUE)
76#endif
77
78/// Our own type!
79static const n_coeffType ID = n_transExt;
80
81/* polynomial ring in which the numerators and denominators of our
82   numbers live */
83#define ntRing cf->extRing
84
85/* coeffs object in which the coefficients of our numbers live;
86 * methods attached to ntCoeffs may be used to compute with the
87 * coefficients of our numbers, e.g., use ntCoeffs->nAdd to add
88 * coefficients of our numbers */
89#define ntCoeffs cf->extRing->cf
90
91
92
93omBin fractionObjectBin = omGetSpecBin(sizeof(fractionObject));
94
95/// forward declarations
96BOOLEAN  ntGreaterZero(number a, const coeffs cf);
97BOOLEAN  ntGreater(number a, number b, const coeffs cf);
98BOOLEAN  ntEqual(number a, number b, const coeffs cf);
99BOOLEAN  ntIsOne(number a, const coeffs cf);
100BOOLEAN  ntIsMOne(number a, const coeffs cf);
101BOOLEAN  ntIsZero(number a, const coeffs cf);
102number   ntInit(long i, const coeffs cf);
103int      ntInt(number &a, const coeffs cf);
104number   ntNeg(number a, const coeffs cf);
105number   ntInvers(number a, const coeffs cf);
106number   ntAdd(number a, number b, const coeffs cf);
107number   ntSub(number a, number b, const coeffs cf);
108number   ntMult(number a, number b, const coeffs cf);
109number   ntDiv(number a, number b, const coeffs cf);
110void     ntPower(number a, int exp, number *b, const coeffs cf);
111number   ntCopy(number a, const coeffs cf);
112void     ntWrite(number &a, const coeffs cf);
113number   ntRePart(number a, const coeffs cf);
114number   ntImPart(number a, const coeffs cf);
115number   ntGetDenom(number &a, const coeffs cf);
116number   ntGetNumerator(number &a, const coeffs cf);
117number   ntGcd(number a, number b, const coeffs cf);
118number   ntLcm(number a, number b, const coeffs cf);
119int      ntSize(number a, const coeffs cf);
120void     ntDelete(number * a, const coeffs cf);
121void     ntCoeffWrite(const coeffs cf, BOOLEAN details);
122number   ntIntDiv(number a, number b, const coeffs cf);
123const char * ntRead(const char *s, number *a, const coeffs cf);
124static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param);
125
126void heuristicGcdCancellation(number a, const coeffs cf);
127void definiteGcdCancellation(number a, const coeffs cf,
128                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed);
129void handleNestedFractionsOverQ(fraction f, const coeffs cf);
130
131#ifdef LDEBUG
132BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs cf)
133{
134  assume(getCoeffType(cf) == ID);
135  fraction t = (fraction)a;
136  if (IS0(t)) return TRUE;
137  assume(NUM(t) != NULL);   /**< t != 0 ==> numerator(t) != 0 */
138  p_Test(NUM(t), ntRing);
139  if (!DENIS1(t)) 
140  {
141    p_Test(DEN(t), ntRing);
142    if(p_IsConstant(DEN(t),ntRing) && (n_IsOne(pGetCoeff(DEN(t)),ntRing->cf)))
143    {
144      Print("?/1 in %s:%d\n",f,l);
145       return FALSE;
146    }
147  }
148  return TRUE;
149}
150#endif
151
152/* returns the bottom field in this field extension tower; if the tower
153   is flat, i.e., if there is no extension, then r itself is returned;
154   as a side-effect, the counter 'height' is filled with the height of
155   the extension tower (in case the tower is flat, 'height' is zero) */
156static coeffs nCoeff_bottom(const coeffs r, int &height)
157{
158  assume(r != NULL);
159  coeffs cf = r;
160  height = 0;
161  while (nCoeff_is_Extension(cf))
162  {
163    assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL);
164    cf = cf->extRing->cf;
165    height++;
166  }
167  return cf;
168}
169
170BOOLEAN ntIsZero(number a, const coeffs cf)
171{
172  ntTest(a);
173  return (IS0(a));
174}
175
176void ntDelete(number * a, const coeffs cf)
177{
178  fraction f = (fraction)(*a);
179  if (IS0(f)) return;
180  p_Delete(&NUM(f), ntRing);
181  if (!DENIS1(f)) p_Delete(&DEN(f), ntRing);
182  omFreeBin((ADDRESS)f, fractionObjectBin);
183  *a = NULL;
184}
185
186BOOLEAN ntEqual(number a, number b, const coeffs cf)
187{
188  ntTest(a); ntTest(b);
189
190  /// simple tests
191  if (a == b) return TRUE;
192  if ((IS0(a)) && (!IS0(b))) return FALSE;
193  if ((IS0(b)) && (!IS0(a))) return FALSE;
194
195  /// cheap test if gcd's have been cancelled in both numbers
196  fraction fa = (fraction)a;
197  fraction fb = (fraction)b;
198  if ((COM(fa) == 1) && (COM(fb) == 1))
199  {
200    poly f = p_Add_q(p_Copy(NUM(fa), ntRing),
201                     p_Neg(p_Copy(NUM(fb), ntRing), ntRing),
202                     ntRing);
203    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
204    if (DENIS1(fa) && DENIS1(fb))  return TRUE;
205    if (DENIS1(fa) && !DENIS1(fb)) return FALSE;
206    if (!DENIS1(fa) && DENIS1(fb)) return FALSE;
207    f = p_Add_q(p_Copy(DEN(fa), ntRing),
208                p_Neg(p_Copy(DEN(fb), ntRing), ntRing),
209                ntRing);
210    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
211    return TRUE;
212  }
213
214  /* default: the more expensive multiplication test
215              a/b = c/d  <==>  a*d = b*c */
216  poly f = p_Copy(NUM(fa), ntRing);
217  if (!DENIS1(fb)) f = p_Mult_q(f, p_Copy(DEN(fb), ntRing), ntRing);
218  poly g = p_Copy(NUM(fb), ntRing);
219  if (!DENIS1(fa)) g = p_Mult_q(g, p_Copy(DEN(fa), ntRing), ntRing);
220  poly h = p_Add_q(f, p_Neg(g, ntRing), ntRing);
221  if (h == NULL) return TRUE;
222  else
223  {
224    p_Delete(&h, ntRing);
225    return FALSE;
226  }
227}
228
229number ntCopy(number a, const coeffs cf)
230{
231  ntTest(a);
232  if (IS0(a)) return NULL;
233  fraction f = (fraction)a;
234  poly g = p_Copy(NUM(f), ntRing);
235  poly h = NULL; if (!DENIS1(f)) h = p_Copy(DEN(f), ntRing);
236  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
237  NUM(result) = g;
238  DEN(result) = h;
239  COM(result) = COM(f);
240  return (number)result;
241}
242
243number ntGetNumerator(number &a, const coeffs cf)
244{
245  ntTest(a);
246  definiteGcdCancellation(a, cf, FALSE);
247  if (IS0(a)) return NULL;
248  fraction f = (fraction)a;
249  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
250  BOOLEAN denis1= DENIS1 (f);
251  if (getCoeffType (ntCoeffs) == n_Q && !denis1)
252    handleNestedFractionsOverQ (f, cf);
253  NUM (result)= p_Copy (NUM (f), ntRing);
254  DEN (result) = NULL;
255  COM (result) = 0;
256  if (getCoeffType (ntCoeffs) == n_Q && denis1)
257  {
258    if (!p_IsConstant (NUM (result), ntRing) && pNext (NUM(result)) != NULL)
259      p_Cleardenom (NUM(result), ntRing);
260    else
261    {
262      number g= p_GetAllDenom (NUM (result), ntRing);
263      NUM (result)= p_Mult_nn (NUM (result), g, ntRing);
264    }
265  }
266  return (number)result;
267}
268
269number ntGetDenom(number &a, const coeffs cf)
270{
271  ntTest(a);
272  definiteGcdCancellation(a, cf, FALSE);
273  fraction f = (fraction)a;
274  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
275  number g;
276  if (IS0(f) || (DENIS1 (f) && getCoeffType (ntCoeffs) != n_Q))
277  {
278    NUM (result)= p_One(ntRing);
279    DEN (result)= NULL;
280    COM (result)= 0;
281  }
282  else if (DENIS1 (f))
283  {
284    poly num= p_Copy (NUM (f), ntRing);
285    if (!p_IsConstant (num, ntRing) && pNext(num) != NULL)
286      p_Cleardenom_n (num, ntRing, g);
287    else
288      g= p_GetAllDenom (num, ntRing);
289    result= (fraction) ntSetMap (ntRing->cf, cf) (g, ntRing->cf, cf);
290  }
291  else
292  {
293    if (getCoeffType (ntCoeffs) == n_Q)
294      handleNestedFractionsOverQ (f, cf);
295    NUM (result)= p_Copy (DEN (f), ntRing);
296    DEN (result) = NULL;
297    COM (result) = 0;
298  }
299  return (number)result;
300}
301
302BOOLEAN ntIsOne(number a, const coeffs cf)
303{
304  ntTest(a);
305  definiteGcdCancellation(a, cf, FALSE);
306  fraction f = (fraction)a;
307  return (f!=NULL) && DENIS1(f) && NUMIS1(f);
308}
309
310BOOLEAN ntIsMOne(number a, const coeffs cf)
311{
312  ntTest(a);
313  definiteGcdCancellation(a, cf, FALSE);
314  fraction f = (fraction)a;
315  if ((f==NULL) || (!DENIS1(f))) return FALSE;
316  poly g = NUM(f);
317  if (!p_IsConstant(g, ntRing)) return FALSE;
318  return n_IsMOne(p_GetCoeff(g, ntRing), ntCoeffs);
319}
320
321/// this is in-place, modifies a
322number ntNeg(number a, const coeffs cf)
323{
324  ntTest(a);
325  if (!IS0(a))
326  {
327    fraction f = (fraction)a;
328    NUM(f) = p_Neg(NUM(f), ntRing);
329  }
330  return a;
331}
332
333number ntImPart(number a, const coeffs cf)
334{
335  ntTest(a);
336  return NULL;
337}
338
339number ntInit_bigint(number longratBigIntNumber, const coeffs src, const coeffs cf)
340{
341  assume( cf != NULL );
342
343  const ring A = cf->extRing;
344
345  assume( A != NULL );
346
347  const coeffs C = A->cf;
348
349  assume( C != NULL );
350
351  number n = n_Init_bigint(longratBigIntNumber, src, C);
352
353  if ( n_IsZero(n, C) )
354  {
355    n_Delete(&n, C);
356    return NULL;
357  }   
358
359  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
360
361  NUM(result) = p_NSet(n, A);
362  DEN(result) = NULL;
363  COM(result) = 0;
364  return (number)result;
365}
366
367   
368number ntInit(long i, const coeffs cf)
369{
370  if (i == 0) return NULL;
371  else
372  {
373    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
374    NUM(result) = p_ISet(i, ntRing);
375    //DEN(result) = NULL; // done by omAlloc0Bin
376    //COM(result) = 0; // done by omAlloc0Bin
377    return (number)result;
378  }
379}
380
381number ntInit(poly p, const coeffs cf)
382{
383  if (p == 0) return NULL;
384  else
385  {
386    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
387    NUM(result) = p;
388    DEN(result) = NULL;
389    COM(result) = 0;
390    return (number)result;
391  }
392}
393
394int ntInt(number &a, const coeffs cf)
395{
396  ntTest(a);
397  if (IS0(a)) return 0;
398  definiteGcdCancellation(a, cf, FALSE);
399  fraction f = (fraction)a;
400  if (!DENIS1(f)) return 0;
401  if (!p_IsConstant(NUM(f), ntRing)) return 0;
402  return n_Int(p_GetCoeff(NUM(f), ntRing), ntCoeffs);
403}
404
405/* This method will only consider the numerators of a and b, without
406   cancelling gcd's before.
407   Moreover it may return TRUE only if one or both numerators
408   are zero or if their degrees are equal. Then TRUE is returned iff
409   coeff(numerator(a)) > coeff(numerator(b));
410   In all other cases, FALSE will be returned. */
411BOOLEAN ntGreater(number a, number b, const coeffs cf)
412{
413  ntTest(a); ntTest(b);
414  number aNumCoeff = NULL; int aNumDeg = 0;
415  number bNumCoeff = NULL; int bNumDeg = 0;
416  if (!IS0(a))
417  {
418    fraction fa = (fraction)a;
419    aNumDeg = p_Totaldegree(NUM(fa), ntRing);
420    aNumCoeff = p_GetCoeff(NUM(fa), ntRing);
421  }
422  if (!IS0(b))
423  {
424    fraction fb = (fraction)b;
425    bNumDeg = p_Totaldegree(NUM(fb), ntRing);
426    bNumCoeff = p_GetCoeff(NUM(fb), ntRing);
427  }
428  if (aNumDeg != bNumDeg) return FALSE;
429  else return n_Greater(aNumCoeff, bNumCoeff, ntCoeffs);
430}
431
432/* this method will only consider the numerator of a, without cancelling
433   the gcd before;
434   returns TRUE iff the leading coefficient of the numerator of a is > 0
435                    or the leading term of the numerator of a is not a
436                    constant */
437BOOLEAN ntGreaterZero(number a, const coeffs cf)
438{
439  ntTest(a);
440  if (IS0(a)) return FALSE;
441  fraction f = (fraction)a;
442  poly g = NUM(f);
443  return (n_GreaterZero(p_GetCoeff(g, ntRing), ntCoeffs) ||
444          (!p_LmIsConstant(g, ntRing)));
445}
446
447void ntCoeffWrite(const coeffs cf, BOOLEAN details)
448{
449  assume( cf != NULL );
450
451  const ring A = cf->extRing;
452
453  assume( A != NULL );
454  assume( A->cf != NULL );
455
456  n_CoeffWrite(A->cf, details);
457
458//  rWrite(A);
459
460  const int P = rVar(A);
461  assume( P > 0 );
462
463  Print("//   %d parameter    : ", P);
464
465  for (int nop=0; nop < P; nop ++)
466    Print("%s ", rRingVar(nop, A));
467
468  assume( A->minideal == NULL );
469
470  PrintS("\n//   minpoly        : 0\n");
471
472/*
473  PrintS("//   Coefficients live in the rational function field\n");
474  Print("//   K(");
475  for (int i = 0; i < rVar(ntRing); i++)
476  {
477    if (i > 0) PrintS(" ");
478    Print("%s", rRingVar(i, ntRing));
479  }
480  PrintS(") with\n");
481  PrintS("//   K: "); n_CoeffWrite(cf->extRing->cf);
482*/
483}
484
485number ntAdd(number a, number b, const coeffs cf)
486{
487  ntTest(a); ntTest(b);
488  if (IS0(a)) return ntCopy(b, cf);
489  if (IS0(b)) return ntCopy(a, cf);
490
491  fraction fa = (fraction)a;
492  fraction fb = (fraction)b;
493
494  poly g = p_Copy(NUM(fa), ntRing);
495  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
496  poly h = p_Copy(NUM(fb), ntRing);
497  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
498  g = p_Add_q(g, h, ntRing);
499
500  if (g == NULL) return NULL;
501
502  poly f;
503  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
504  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
505  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
506  else /* both denom's are != 1 */    f = p_Mult_q(p_Copy(DEN(fa), ntRing),
507                                                   p_Copy(DEN(fb), ntRing),
508                                                   ntRing);
509
510  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
511  NUM(result) = g;
512  DEN(result) = f;
513  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
514  heuristicGcdCancellation((number)result, cf);
515  return (number)result;
516}
517
518number ntSub(number a, number b, const coeffs cf)
519{
520  ntTest(a); ntTest(b);
521  if (IS0(a)) return ntNeg(ntCopy(b, cf), cf);
522  if (IS0(b)) return ntCopy(a, cf);
523
524  fraction fa = (fraction)a;
525  fraction fb = (fraction)b;
526
527  poly g = p_Copy(NUM(fa), ntRing);
528  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
529  poly h = p_Copy(NUM(fb), ntRing);
530  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
531  g = p_Add_q(g, p_Neg(h, ntRing), ntRing);
532
533  if (g == NULL) return NULL;
534
535  poly f;
536  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
537  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
538  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
539  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
540                                                   p_Copy(DEN(fb), ntRing),
541                                                   ntRing);
542
543  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
544  NUM(result) = g;
545  DEN(result) = f;
546  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
547  heuristicGcdCancellation((number)result, cf);
548  return (number)result;
549}
550
551number ntMult(number a, number b, const coeffs cf)
552{
553  ntTest(a); ntTest(b);
554  if (IS0(a) || IS0(b)) return NULL;
555
556  fraction fa = (fraction)a;
557  fraction fb = (fraction)b;
558
559  poly g = p_Copy(NUM(fa), ntRing);
560  poly h = p_Copy(NUM(fb), ntRing);
561  g = p_Mult_q(g, h, ntRing);
562
563  if (g == NULL) return NULL;   /* may happen due to zero divisors */
564
565  poly f;
566  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
567  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
568  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
569  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
570                                                   p_Copy(DEN(fb), ntRing),
571                                                   ntRing);
572
573  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
574  NUM(result) = g;
575  DEN(result) = f;
576  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
577  heuristicGcdCancellation((number)result, cf);
578  return (number)result;
579}
580
581number ntDiv(number a, number b, const coeffs cf)
582{
583  ntTest(a); ntTest(b);
584  if (IS0(a)) return NULL;
585  if (IS0(b)) WerrorS(nDivBy0);
586
587  fraction fa = (fraction)a;
588  fraction fb = (fraction)b;
589
590  poly g = p_Copy(NUM(fa), ntRing);
591  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
592
593  if (g == NULL) return NULL;   /* may happen due to zero divisors */
594
595  poly f = p_Copy(NUM(fb), ntRing);
596  if (!DENIS1(fa)) f = p_Mult_q(f, p_Copy(DEN(fa), ntRing), ntRing);
597
598  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
599  NUM(result) = g;
600  if (!p_IsConstant(f,ntRing) || !n_IsOne(pGetCoeff(f),ntRing->cf))
601    DEN(result) = f;
602  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
603  heuristicGcdCancellation((number)result, cf);
604  return (number)result;
605}
606
607/* 0^0 = 0;
608   for |exp| <= 7 compute power by a simple multiplication loop;
609   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
610      p^13 = p^1 * p^4 * p^8, where we utilise that
611      p^(2^(k+1)) = p^(2^k) * p^(2^k);
612   intermediate cancellation is controlled by the in-place method
613   heuristicGcdCancellation; see there.
614*/
615void ntPower(number a, int exp, number *b, const coeffs cf)
616{
617  ntTest(a);
618
619  /* special cases first */
620  if (IS0(a))
621  {
622    if (exp >= 0) *b = NULL;
623    else          WerrorS(nDivBy0);
624  }
625  else if (exp ==  0) { *b = ntInit(1, cf); return;}
626  else if (exp ==  1) { *b = ntCopy(a, cf); return;}
627  else if (exp == -1) { *b = ntInvers(a, cf); return;}
628
629  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
630
631  /* now compute a^expAbs */
632  number pow; number t;
633  if (expAbs <= 7)
634  {
635    pow = ntCopy(a, cf);
636    for (int i = 2; i <= expAbs; i++)
637    {
638      t = ntMult(pow, a, cf);
639      ntDelete(&pow, cf);
640      pow = t;
641      heuristicGcdCancellation(pow, cf);
642    }
643  }
644  else
645  {
646    pow = ntInit(1, cf);
647    number factor = ntCopy(a, cf);
648    while (expAbs != 0)
649    {
650      if (expAbs & 1)
651      {
652        t = ntMult(pow, factor, cf);
653        ntDelete(&pow, cf);
654        pow = t;
655        heuristicGcdCancellation(pow, cf);
656      }
657      expAbs = expAbs / 2;
658      if (expAbs != 0)
659      {
660        t = ntMult(factor, factor, cf);
661        ntDelete(&factor, cf);
662        factor = t;
663        heuristicGcdCancellation(factor, cf);
664      }
665    }
666    ntDelete(&factor, cf);
667  }
668
669  /* invert if original exponent was negative */
670  if (exp < 0)
671  {
672    t = ntInvers(pow, cf);
673    ntDelete(&pow, cf);
674    pow = t;
675  }
676  *b = pow;
677}
678
679/* assumes that cf represents the rationals, i.e. Q, and will only
680   be called in that case;
681   assumes furthermore that f != NULL and that the denominator of f != 1;
682   generally speaking, this method removes denominators in the rational
683   coefficients of the numerator and denominator of 'a';
684   more concretely, the following normalizations will be performed,
685   where t^alpha denotes a monomial in the transcendental variables t_k
686   (1) if 'a' is of the form
687          (sum_alpha a_alpha/b_alpha * t^alpha)
688          -------------------------------------
689            (sum_beta c_beta/d_beta * t^beta)
690       with integers a_alpha, b_alpha, c_beta, d_beta, then both the
691       numerator and the denominator will be multiplied by the LCM of
692       the b_alpha's and the d_beta's (if this LCM is != 1),
693   (2) if 'a' is - e.g. after having performed step (1) - of the form
694          (sum_alpha a_alpha * t^alpha)
695          -----------------------------
696            (sum_beta c_beta * t^beta)
697       with integers a_alpha, c_beta, and with a non-constant denominator,
698       then both the numerator and the denominator will be divided by the
699       GCD of the a_alpha's and the c_beta's (if this GCD is != 1),
700   (3) if 'a' is - e.g. after having performed steps (1) and (2) - of the
701       form
702          (sum_alpha a_alpha * t^alpha)
703          -----------------------------
704                        c
705       with integers a_alpha, and c != 1, then 'a' will be replaced by
706       (sum_alpha a_alpha/c * t^alpha);
707   this procedure does not alter COM(f) (this has to be done by the
708   calling procedure);
709   modifies f */
710void handleNestedFractionsOverQ(fraction f, const coeffs cf)
711{
712  assume(nCoeff_is_Q(ntCoeffs));
713  assume(!IS0(f));
714  assume(!DENIS1(f));
715
716  if (!p_IsConstant(DEN(f), ntRing))
717  { /* step (1); see documentation of this procedure above */
718    p_Normalize(NUM(f), ntRing);
719    p_Normalize(DEN(f), ntRing);
720    number lcmOfDenominators = n_Init(1, ntCoeffs);
721    number c; number tmp;
722    poly p = NUM(f);
723    /* careful when using n_Lcm!!! It computes the lcm of the numerator
724       of the 1st argument and the denominator of the 2nd!!! */
725    while (p != NULL)
726    {
727      c = p_GetCoeff(p, ntRing);
728      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
729      n_Delete(&lcmOfDenominators, ntCoeffs);
730      lcmOfDenominators = tmp;
731      pIter(p);
732    }
733    p = DEN(f);
734    while (p != NULL)
735    {
736      c = p_GetCoeff(p, ntRing);
737      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
738      n_Delete(&lcmOfDenominators, ntCoeffs);
739      lcmOfDenominators = tmp;
740      pIter(p);
741    }
742    if (!n_IsOne(lcmOfDenominators, ntCoeffs))
743    { /* multiply NUM(f) and DEN(f) with lcmOfDenominators */
744      NUM(f) = p_Mult_nn(NUM(f), lcmOfDenominators, ntRing);
745      p_Normalize(NUM(f), ntRing);
746      DEN(f) = p_Mult_nn(DEN(f), lcmOfDenominators, ntRing);
747      p_Normalize(DEN(f), ntRing);
748    }
749    n_Delete(&lcmOfDenominators, ntCoeffs);
750    if (!p_IsConstant(DEN(f), ntRing))
751    { /* step (2); see documentation of this procedure above */
752      p = NUM(f);
753      number gcdOfCoefficients = n_Copy(p_GetCoeff(p, ntRing), ntCoeffs);
754      pIter(p);
755      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
756      {
757        c = p_GetCoeff(p, ntRing);
758        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
759        n_Delete(&gcdOfCoefficients, ntCoeffs);
760        gcdOfCoefficients = tmp;
761        pIter(p);
762      }
763      p = DEN(f);
764      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
765      {
766        c = p_GetCoeff(p, ntRing);
767        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
768        n_Delete(&gcdOfCoefficients, ntCoeffs);
769        gcdOfCoefficients = tmp;
770        pIter(p);
771      }
772      if (!n_IsOne(gcdOfCoefficients, ntCoeffs))
773      { /* divide NUM(f) and DEN(f) by gcdOfCoefficients */
774        number inverseOfGcdOfCoefficients = n_Invers(gcdOfCoefficients,
775                                                     ntCoeffs);
776        NUM(f) = p_Mult_nn(NUM(f), inverseOfGcdOfCoefficients, ntRing);
777        p_Normalize(NUM(f), ntRing);
778        DEN(f) = p_Mult_nn(DEN(f), inverseOfGcdOfCoefficients, ntRing);
779        p_Normalize(DEN(f), ntRing);
780        n_Delete(&inverseOfGcdOfCoefficients, ntCoeffs);
781      }
782      n_Delete(&gcdOfCoefficients, ntCoeffs);
783    }
784  }
785  if (p_IsConstant(DEN(f), ntRing) &&
786      (!n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs)))
787  { /* step (3); see documentation of this procedure above */
788    number inverseOfDen = n_Invers(p_GetCoeff(DEN(f), ntRing), ntCoeffs);
789    NUM(f) = p_Mult_nn(NUM(f), inverseOfDen, ntRing);
790    n_Delete(&inverseOfDen, ntCoeffs);
791    p_Delete(&DEN(f), ntRing);
792    DEN(f) = NULL;
793  }
794
795  /* Now, due to the above computations, DEN(f) may have become the
796     1-polynomial which needs to be represented by NULL: */
797  if ((DEN(f) != NULL) &&
798      p_IsConstant(DEN(f), ntRing) &&
799      n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
800  {
801    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
802  }
803}
804
805/* modifies a */
806void heuristicGcdCancellation(number a, const coeffs cf)
807{
808  ntTest(a);
809  if (IS0(a)) return;
810
811  fraction f = (fraction)a;
812  if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
813
814  /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
815  if (p_EqualPolys(NUM(f), DEN(f), ntRing))
816  { /* numerator and denominator are both != 1 */
817    p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
818    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
819    COM(f) = 0;
820    return;
821  }
822
823  if (COM(f) <= BOUND_COMPLEXITY) return;
824  else definiteGcdCancellation(a, cf, TRUE);
825}
826
827/* modifies a */
828void definiteGcdCancellation(number a, const coeffs cf,
829                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed)
830{
831  ntTest(a);
832
833  fraction f = (fraction)a;
834
835  if (!simpleTestsHaveAlreadyBeenPerformed)
836  {
837    if (IS0(a)) return;
838    if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
839
840    /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
841    if (p_EqualPolys(NUM(f), DEN(f), ntRing))
842    { /* numerator and denominator are both != 1 */
843      p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
844      p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
845      COM(f) = 0;
846      return;
847    }
848  }
849
850#ifdef HAVE_FACTORY
851  /* singclap_gcd destroys its arguments; we hence need copies: */
852  poly pNum = p_Copy(NUM(f), ntRing);
853  poly pDen = p_Copy(DEN(f), ntRing);
854
855  /* Note that, over Q, singclap_gcd will remove the denominators in all
856     rational coefficients of pNum and pDen, before starting to compute
857     the gcd. Thus, we do not need to ensure that the coefficients of
858     pNum and pDen live in Z; they may well be elements of Q\Z. */
859  poly pGcd = singclap_gcd(pNum, pDen, cf->extRing);
860  if (p_IsConstant(pGcd, ntRing) &&
861      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
862  { /* gcd = 1; nothing to cancel;
863       Suppose the given rational function field is over Q. Although the
864       gcd is 1, we may have produced fractional coefficients in NUM(f),
865       DEN(f), or both, due to previous arithmetics. The next call will
866       remove those nested fractions, in case there are any. */
867    if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
868  }
869  else
870  { /* We divide both NUM(f) and DEN(f) by the gcd which is known
871       to be != 1. */
872    poly newNum = singclap_pdivide(NUM(f), pGcd, ntRing);
873    p_Delete(&NUM(f), ntRing);
874    NUM(f) = newNum;
875    poly newDen = singclap_pdivide(DEN(f), pGcd, ntRing);
876    p_Delete(&DEN(f), ntRing);
877    DEN(f) = newDen;
878    if (p_IsConstant(DEN(f), ntRing) &&
879        n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
880    {
881      /* DEN(f) = 1 needs to be represented by NULL! */
882      p_Delete(&DEN(f), ntRing);
883      newDen = NULL;
884    }
885    else
886    { /* Note that over Q, by cancelling the gcd, we may have produced
887         fractional coefficients in NUM(f), DEN(f), or both. The next
888         call will remove those nested fractions, in case there are
889         any. */
890      if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
891    }
892  }
893  COM(f) = 0;
894  p_Delete(&pGcd, ntRing);
895#endif /* HAVE_FACTORY */
896}
897
898/* modifies a */
899void ntWrite(number &a, const coeffs cf)
900{
901  ntTest(a);
902  definiteGcdCancellation(a, cf, FALSE);
903  if (IS0(a))
904    StringAppendS("0");
905  else
906  {
907    fraction f = (fraction)a;
908    // stole logic from napWrite from kernel/longtrans.cc of legacy singular
909    BOOLEAN omitBrackets = p_IsConstant(NUM(f), ntRing);
910    if (!omitBrackets) StringAppendS("(");
911    p_String0(NUM(f), ntRing, ntRing);
912    if (!omitBrackets) StringAppendS(")");
913    if (!DENIS1(f))
914    {
915      StringAppendS("/");
916      omitBrackets = p_IsConstant(DEN(f), ntRing);
917      if (!omitBrackets) StringAppendS("(");
918      p_String0(DEN(f), ntRing, ntRing);
919      if (!omitBrackets) StringAppendS(")");
920    }
921  }
922}
923
924const char * ntRead(const char *s, number *a, const coeffs cf)
925{
926  poly p;
927  const char * result = p_Read(s, p, ntRing);
928  if (p == NULL) { *a = NULL; return result; }
929  else
930  {
931    fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
932    NUM(f) = p;
933    DEN(f) = NULL;
934    COM(f) = 0;
935    *a = (number)f;
936    return result;
937  }
938}
939
940void ntNormalize (number &a, const coeffs cf)
941{
942  definiteGcdCancellation(a, cf, FALSE);
943}
944
945/* expects *param to be castable to TransExtInfo */
946static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
947{
948  if (ID != n) return FALSE;
949  TransExtInfo *e = (TransExtInfo *)param;
950  /* for rational function fields we expect the underlying
951     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
952     this expectation is based on the assumption that we have properly
953     registered cf and perform reference counting rather than creating
954     multiple copies of the same coefficient field/domain/ring */
955  return (ntRing == e->r);
956}
957
958number ntLcm(number a, number b, const coeffs cf)
959{
960  ntTest(a); ntTest(b);
961  fraction fb = (fraction)b;
962  if ((b==NULL)||(DEN(fb)==NULL)) return ntCopy(a,cf);
963#ifdef HAVE_FACTORY
964  fraction fa = (fraction)a;
965  /* singclap_gcd destroys its arguments; we hence need copies: */
966  poly pa = p_Copy(NUM(fa), ntRing);
967  poly pb = p_Copy(DEN(fb), ntRing);
968
969  /* Note that, over Q, singclap_gcd will remove the denominators in all
970     rational coefficients of pa and pb, before starting to compute
971     the gcd. Thus, we do not need to ensure that the coefficients of
972     pa and pb live in Z; they may well be elements of Q\Z. */
973  poly pGcd = singclap_gcd(pa, pb, cf->extRing);
974  if (p_IsConstant(pGcd, ntRing) &&
975      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
976  { /* gcd = 1; return pa*pb*/
977    p_Delete(&pGcd,ntRing);
978    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
979    NUM(result) = pp_Mult_qq(NUM(fa),DEN(fb),ntRing);
980    return (number)result;
981  }
982  else
983  { /* return pa*pb/gcd */
984    poly newNum = singclap_pdivide(NUM(fa), pGcd, ntRing);
985    p_Delete(&pGcd,ntRing);
986    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
987    NUM(result) = p_Mult_q(p_Copy(DEN(fb),ntRing),newNum,ntRing);
988    return (number)result;
989  }
990#else
991  Print("// factory needed: transext.cc:ntLcm\n");
992  return NULL;
993#endif /* HAVE_FACTORY */
994  return NULL;
995}
996
997number ntGcd(number a, number b, const coeffs cf)
998{
999  ntTest(a); ntTest(b);
1000  if (a==NULL) return ntCopy(b,cf);
1001  if (b==NULL) return ntCopy(a,cf);
1002#ifdef HAVE_FACTORY
1003  fraction fa = (fraction)a;
1004  fraction fb = (fraction)b;
1005  /* singclap_gcd destroys its arguments; we hence need copies: */
1006  poly pa = p_Copy(NUM(fa), ntRing);
1007  poly pb = p_Copy(NUM(fb), ntRing);
1008
1009  /* Note that, over Q, singclap_gcd will remove the denominators in all
1010     rational coefficients of pa and pb, before starting to compute
1011     the gcd. Thus, we do not need to ensure that the coefficients of
1012     pa and pb live in Z; they may well be elements of Q\Z. */
1013  poly pGcd = singclap_gcd(pa, pb, cf->extRing);
1014  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1015  NUM(result) = pGcd;
1016  return (number)result;
1017#else
1018  Print("// factory needed: transext.cc:ntGcd\n");
1019  return NULL;
1020#endif /* HAVE_FACTORY */
1021}
1022
1023int ntSize(number a, const coeffs cf)
1024{
1025  ntTest(a);
1026  if (IS0(a)) return -1;
1027  /* this has been taken from the old implementation of field extensions,
1028     where we computed the sum of the degrees and the numbers of terms in
1029     the numerator and denominator of a; so we leave it at that, for the
1030     time being */
1031  fraction f = (fraction)a;
1032  poly p = NUM(f);
1033  int noOfTerms = 0;
1034  int numDegree = 0;
1035  while (p != NULL)
1036  {
1037    noOfTerms++;
1038    int d = 0;
1039    for (int i = 1; i <= rVar(ntRing); i++)
1040      d += p_GetExp(p, i, ntRing);
1041    if (d > numDegree) numDegree = d;
1042    pIter(p);
1043  }
1044  int denDegree = 0;
1045  if (!DENIS1(f))
1046  {
1047    p = DEN(f);
1048    while (p != NULL)
1049    {
1050      noOfTerms++;
1051      int d = 0;
1052      for (int i = 1; i <= rVar(ntRing); i++)
1053        d += p_GetExp(p, i, ntRing);
1054      if (d > denDegree) denDegree = d;
1055      pIter(p);
1056    }
1057  }
1058  return numDegree + denDegree + noOfTerms;
1059}
1060
1061number ntInvers(number a, const coeffs cf)
1062{
1063  ntTest(a);
1064  if (IS0(a)) WerrorS(nDivBy0);
1065  fraction f = (fraction)a;
1066  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1067  poly g;
1068  if (DENIS1(f)) g = p_One(ntRing);
1069  else           g = p_Copy(DEN(f), ntRing);
1070  NUM(result) = g;
1071  DEN(result) = p_Copy(NUM(f), ntRing);
1072  COM(result) = COM(f);
1073  return (number)result;
1074}
1075
1076/* assumes that src = Q, dst = Q(t_1, ..., t_s) */
1077number ntMap00(number a, const coeffs src, const coeffs dst)
1078{
1079  if (n_IsZero(a, src)) return NULL;
1080  assume(src == dst->extRing->cf);
1081  poly p = p_Init(dst->extRing);
1082  number na=n_Copy(a, src);
1083  n_Normalize(na, src);
1084  p_SetCoeff0(p, na, dst->extRing);
1085  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1086  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1087  return (number)f;
1088}
1089
1090/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
1091number ntMapP0(number a, const coeffs src, const coeffs dst)
1092{
1093  if (n_IsZero(a, src)) return NULL;
1094  /* mapping via intermediate int: */
1095  int n = n_Int(a, src);
1096  number q = n_Init(n, dst->extRing->cf);
1097  poly p;
1098  if (n_IsZero(q, dst->extRing->cf))
1099  {
1100    n_Delete(&q, dst->extRing->cf);
1101    return NULL;
1102  }
1103  p = p_One(dst->extRing);
1104  p_SetCoeff(p, q, dst->extRing);
1105  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1106  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1107  return (number)f;
1108}
1109
1110/* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
1111                       src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
1112number ntCopyMap(number a, const coeffs cf, const coeffs dst)
1113{
1114//  if (n_IsZero(a, cf)) return NULL;
1115   
1116  ntTest(a);
1117
1118  if (IS0(a)) return NULL;
1119   
1120  const ring rSrc = cf->extRing;
1121  const ring rDst = dst->extRing;
1122 
1123  if( rSrc == rDst )
1124    return ntCopy(a, dst); // USUALLY WRONG!
1125 
1126  fraction f = (fraction)a;
1127  poly g = prCopyR(NUM(f), rSrc, rDst);
1128   
1129  poly h = NULL;
1130 
1131  if (!DENIS1(f))
1132     h = prCopyR(DEN(f), rSrc, rDst);
1133   
1134  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1135   
1136  NUM(result) = g;
1137  DEN(result) = h;
1138  COM(result) = COM(f);
1139  return (number)result; 
1140}
1141
1142number ntCopyAlg(number a, const coeffs cf, const coeffs dst)
1143{
1144  if (n_IsZero(a, cf)) return NULL;
1145   
1146  fraction f = (fraction)omAlloc0Bin(fractionObjectBin); 
1147  // DEN(f) = NULL; COM(f) = 0;
1148  NUM(f) = prCopyR((poly)a, cf->extRing, dst->extRing);
1149  return (number)f;
1150}
1151
1152/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
1153number ntMap0P(number a, const coeffs src, const coeffs dst)
1154{
1155  if (n_IsZero(a, src)) return NULL;
1156  int p = rChar(dst->extRing);
1157  number q = nlModP(a, src, dst->extRing->cf);
1158
1159  if (n_IsZero(q, dst->extRing->cf))
1160  {
1161    n_Delete(&q, dst->extRing->cf);
1162    return NULL;
1163  }
1164 
1165  poly g = p_NSet(q, dst->extRing);
1166
1167  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1168  NUM(f) = g; // DEN(f) = NULL; COM(f) = 0;
1169  return (number)f;
1170}
1171
1172/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
1173number ntMapPP(number a, const coeffs src, const coeffs dst)
1174{
1175  if (n_IsZero(a, src)) return NULL;
1176  assume(src == dst->extRing->cf);
1177  poly p = p_One(dst->extRing);
1178  p_SetCoeff(p, n_Copy(a, src), dst->extRing);
1179  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1180  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1181  return (number)f;
1182}
1183
1184/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
1185number ntMapUP(number a, const coeffs src, const coeffs dst)
1186{
1187  if (n_IsZero(a, src)) return NULL;
1188  /* mapping via intermediate int: */
1189  int n = n_Int(a, src);
1190  number q = n_Init(n, dst->extRing->cf);
1191  poly p;
1192  if (n_IsZero(q, dst->extRing->cf))
1193  {
1194    n_Delete(&q, dst->extRing->cf);
1195    return NULL;
1196  }
1197  p = p_One(dst->extRing);
1198  p_SetCoeff(p, q, dst->extRing);
1199  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1200  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1201  return (number)f;
1202}
1203
1204nMapFunc ntSetMap(const coeffs src, const coeffs dst)
1205{
1206  /* dst is expected to be a rational function field */
1207  assume(getCoeffType(dst) == ID);
1208
1209  int h = 0; /* the height of the extension tower given by dst */
1210  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
1211  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
1212
1213  /* for the time being, we only provide maps if h = 1 and if b is Q or
1214     some field Z/pZ: */
1215  if (h==0)
1216  {
1217    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
1218      return ntMap00;                                 /// Q       -->  Q(T)
1219    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
1220      return ntMapP0;                                 /// Z/p     -->  Q(T)
1221    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
1222      return ntMap0P;                                 /// Q       --> Z/p(T)
1223    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
1224    {
1225      if (src->ch == dst->ch) return ntMapPP;         /// Z/p     --> Z/p(T)
1226      else return ntMapUP;                            /// Z/u     --> Z/p(T)
1227    }
1228  }
1229  if (h != 1) return NULL;
1230  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
1231
1232  /* Let T denote the sequence of transcendental extension variables, i.e.,
1233     K[t_1, ..., t_s] =: K[T];
1234     Let moreover, for any such sequence T, T' denote any subsequence of T
1235     of the form t_1, ..., t_w with w <= s. */
1236
1237  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
1238
1239  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
1240  {
1241    if (rVar(src->extRing) > rVar(dst->extRing)) 
1242       return NULL;
1243     
1244    for (int i = 0; i < rVar(src->extRing); i++)
1245      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0) 
1246         return NULL;
1247     
1248    if (src->type==n_transExt)
1249       return ntCopyMap;          /// Q(T')   --> Q(T)
1250    else
1251       return ntCopyAlg;
1252  }
1253
1254  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
1255  {
1256    if (rVar(src->extRing) > rVar(dst->extRing)) 
1257       return NULL;
1258     
1259    for (int i = 0; i < rVar(src->extRing); i++)
1260      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0) 
1261         return NULL;
1262     
1263    if (src->type==n_transExt)
1264       return ntCopyMap;         /// Z/p(T') --> Z/p(T)
1265    else
1266       return ntCopyAlg;
1267  }
1268
1269  return NULL;                                 /// default
1270}
1271#if 0
1272nMapFunc ntSetMap_T(const coeffs src, const coeffs dst)
1273{
1274  nMapFunc n=ntSetMap(src,dst);
1275  if (n==ntCopyAlg) printf("n=ntCopyAlg\n");
1276  else if (n==ntCopyMap) printf("n=ntCopyMap\n");
1277  else if (n==ntMapUP) printf("n=ntMapUP\n");
1278  else if (n==ntMap0P) printf("n=ntMap0P\n");
1279  else if (n==ntMapP0) printf("n=ntMapP0\n");
1280  else if (n==ntMap00) printf("n=ntMap00\n");
1281  else if (n==NULL) printf("n=NULL\n");
1282  else printf("n=?\n");
1283  return n;
1284}
1285#endif
1286
1287void ntKillChar(coeffs cf)
1288{
1289  if ((--cf->extRing->ref) == 0)
1290    rDelete(cf->extRing);
1291}
1292#ifdef HAVE_FACTORY
1293number ntConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
1294{
1295  if (n.isZero()) return NULL;
1296  poly p=convFactoryPSingP(n,ntRing);
1297  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1298  NUM(result) = p;
1299  //DEN(result) = NULL; // done by omAlloc0Bin
1300  //COM(result) = 0; // done by omAlloc0Bin
1301  return (number)result;
1302}
1303CanonicalForm ntConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs cf )
1304{
1305  ntTest(n);
1306  if (IS0(n)) return CanonicalForm(0);
1307
1308  fraction f = (fraction)n;
1309  return convSingPFactoryP(NUM(f),ntRing);
1310}
1311#endif
1312
1313int ntParDeg(number a, const coeffs cf)
1314{
1315  if (IS0(a)) return -1;
1316  fraction fa = (fraction)a;
1317  return cf->extRing->pFDeg(NUM(fa),cf->extRing);
1318}
1319
1320BOOLEAN ntInitChar(coeffs cf, void * infoStruct)
1321{
1322
1323  assume( infoStruct != NULL );
1324
1325  TransExtInfo *e = (TransExtInfo *)infoStruct;
1326
1327  assume( e->r                != NULL);      // extRing;
1328  assume( e->r->cf            != NULL);      // extRing->cf;
1329  assume( e->r->minideal == NULL );
1330
1331  assume( cf != NULL );
1332  assume(getCoeffType(cf) == ID);                // coeff type;
1333
1334  cf->extRing           = e->r;
1335  cf->extRing->ref ++; // increase the ref.counter for the ground poly. ring!
1336  cf->factoryVarOffset = cf->extRing->cf->factoryVarOffset+rVar(cf->extRing);
1337
1338  /* propagate characteristic up so that it becomes
1339     directly accessible in cf: */
1340  cf->ch = cf->extRing->cf->ch;
1341
1342  cf->cfGreaterZero  = ntGreaterZero;
1343  cf->cfGreater      = ntGreater;
1344  cf->cfEqual        = ntEqual;
1345  cf->cfIsZero       = ntIsZero;
1346  cf->cfIsOne        = ntIsOne;
1347  cf->cfIsMOne       = ntIsMOne;
1348  cf->cfInit         = ntInit;
1349  cf->cfInit_bigint  = ntInit_bigint;
1350  cf->cfInt          = ntInt;
1351  cf->cfNeg          = ntNeg;
1352  cf->cfAdd          = ntAdd;
1353  cf->cfSub          = ntSub;
1354  cf->cfMult         = ntMult;
1355  cf->cfDiv          = ntDiv;
1356  cf->cfExactDiv     = ntDiv;
1357  cf->cfPower        = ntPower;
1358  cf->cfCopy         = ntCopy;
1359  cf->cfWrite        = ntWrite;
1360  cf->cfRead         = ntRead;
1361  cf->cfNormalize    = ntNormalize;
1362  cf->cfDelete       = ntDelete;
1363  cf->cfSetMap       = ntSetMap;
1364  cf->cfGetDenom     = ntGetDenom;
1365  cf->cfGetNumerator = ntGetNumerator;
1366  cf->cfRePart       = ntCopy;
1367  cf->cfImPart       = ntImPart;
1368  cf->cfCoeffWrite   = ntCoeffWrite;
1369#ifdef LDEBUG
1370  cf->cfDBTest       = ntDBTest;
1371#endif
1372  cf->cfGcd          = ntGcd;
1373  cf->cfLcm          = ntLcm;
1374  cf->cfSize         = ntSize;
1375  cf->nCoeffIsEqual  = ntCoeffIsEqual;
1376  cf->cfInvers       = ntInvers;
1377  cf->cfIntDiv       = ntDiv;
1378  cf->cfKillChar     = ntKillChar;
1379
1380#ifndef HAVE_FACTORY
1381  PrintS("// Warning: The 'factory' module is not available.\n");
1382  PrintS("//          Hence gcd's cannot be cancelled in any\n");
1383  PrintS("//          computed fraction!\n");
1384#else
1385  cf->convFactoryNSingN =ntConvFactoryNSingN;
1386  cf->convSingNFactoryN =ntConvSingNFactoryN;
1387#endif
1388  cf->cfParDeg = ntParDeg;
1389
1390  return FALSE;
1391}
1392
1393
1394number ntParam(const short iParameter, const coeffs cf)
1395{
1396  assume(getCoeffType(cf) == ID);
1397
1398  const ring R = cf->extRing;
1399  assume( R != NULL );
1400  assume( 0 < iParameter && iParameter <= rVar(R) );
1401
1402  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
1403
1404//  return (number) p;
1405
1406  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1407  NUM(f) = p;
1408  DEN(f) = NULL;
1409  COM(f) = 0;
1410
1411  return (number)f;
1412}
1413
1414
1415/// if m == var(i)/1 => return i,
1416int ntIsParam(number m, const coeffs cf)
1417{
1418  assume(getCoeffType(cf) == ID);
1419
1420  const ring R = cf->extRing;
1421  assume( R != NULL );
1422
1423  fraction f = (fraction)m;
1424
1425  if( DEN(f) != NULL )
1426    return 0;
1427
1428  return p_Var( NUM(f), R );
1429}
Note: See TracBrowser for help on using the repository browser.