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

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