source: git/libpolys/polys/ext_fields/transext.cc @ 146c603

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