source: git/libpolys/polys/ext_fields/transext.cc @ 644f81

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