source: git/libpolys/polys/ext_fields/transext.cc @ 9a78eb

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