source: git/libpolys/polys/ext_fields/transext.cc @ 44d5ad

spielwiese
Last change on this file since 44d5ad was 0afa07, checked in by Frank Seelisch <seelisch@…>, 13 years ago
removed cfPar, cfParDeg, n_Par, n_ParDeg, ndPar, ndParDeg from coeffs
  • Property mode set to 100644
File size: 34.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: numbers in a rational function field K(t_1, .., t_s) with
7*           transcendental variables t_1, ..., t_s, where s >= 1.
8*           Denoting the implemented coeffs object by cf, then these numbers
9*           are represented as quotients of polynomials living in the
10*           polynomial ring K[t_1, .., t_s] represented by cf->extring.
11*
12*           An element of K(t_1, .., t_s) may have numerous representations,
13*           due to the possibility of common polynomial factors in the
14*           numerator and denominator. This problem is handled by a
15*           cancellation heuristic: Each number "knows" its complexity
16*           which is 0 if and only if common factors have definitely been
17*           cancelled, and some positive integer otherwise.
18*           Each arithmetic operation of two numbers with complexities c1
19*           and c2 will result in a number of complexity c1 + c2 + some
20*           penalty (specific for each arithmetic operation; see constants
21*           in the *.h file). Whenever the resulting complexity exceeds a
22*           certain threshold (see constant in the *.h file), then the
23*           cancellation heuristic will call 'factory' to compute the gcd
24*           and cancel it out in the given number. (This definite cancel-
25*           lation will also be performed at the beginning of ntWrite,
26*           ensuring that any output is free of common factors.
27*           For the special case of K = Q (i.e., when computing over the
28*           rationals), this definite cancellation procedure will also take
29*           care of nested fractions: If there are fractional coefficients
30*           in the numerator or denominator of a number, then this number
31*           is being replaced by a quotient of two polynomials over Z, or
32*           - if the denominator is a constant - by a polynomial over Q.
33*/
34
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#endif
53
54#include <polys/ext_fields/transext.h>
55
56/// our type has been defined as a macro in transext.h
57/// and is accessible by 'ntID'
58
59extern omBin fractionObjectBin = omGetSpecBin(sizeof(fractionObject));
60
61/// forward declarations
62BOOLEAN  ntGreaterZero(number a, const coeffs cf); 
63BOOLEAN  ntGreater(number a, number b, const coeffs cf);
64BOOLEAN  ntEqual(number a, number b, const coeffs cf);
65BOOLEAN  ntIsOne(number a, const coeffs cf);
66BOOLEAN  ntIsMOne(number a, const coeffs cf);
67BOOLEAN  ntIsZero(number a, const coeffs cf);
68number   ntInit(int i, const coeffs cf);
69int      ntInt(number &a, const coeffs cf);
70number   ntNeg(number a, const coeffs cf);
71number   ntInvers(number a, const coeffs cf);
72number   ntAdd(number a, number b, const coeffs cf);
73number   ntSub(number a, number b, const coeffs cf);
74number   ntMult(number a, number b, const coeffs cf);
75number   ntDiv(number a, number b, const coeffs cf);
76void     ntPower(number a, int exp, number *b, const coeffs cf);
77number   ntCopy(number a, const coeffs cf);
78void     ntWrite(number &a, const coeffs cf);
79number   ntRePart(number a, const coeffs cf);
80number   ntImPart(number a, const coeffs cf);
81number   ntGetDenom(number &a, const coeffs cf);
82number   ntGetNumerator(number &a, const coeffs cf);
83number   ntGcd(number a, number b, const coeffs cf);
84number   ntLcm(number a, number b, const coeffs cf);
85int      ntSize(number a, const coeffs cf);
86void     ntDelete(number * a, const coeffs cf);
87void     ntCoeffWrite(const coeffs cf);
88number   ntIntDiv(number a, number b, const coeffs cf);
89const char * ntRead(const char *s, number *a, const coeffs cf);
90static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param);
91
92void heuristicGcdCancellation(number a, const coeffs cf);
93void definiteGcdCancellation(number a, const coeffs cf,
94                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed);
95void handleNestedFractionsOverQ(fraction f, const coeffs cf);
96
97#ifdef LDEBUG
98BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs cf)
99{
100  assume(getCoeffType(cf) == ntID);
101  fraction t = (fraction)a;
102  if (IS0(t)) return TRUE;
103  assume(NUM(t) != NULL);   /**< t != 0 ==> numerator(t) != 0 */
104  p_Test(NUM(t), ntRing);
105  if (!DENIS1(t)) p_Test(DEN(t), ntRing);
106  return TRUE;
107}
108#endif
109
110/* returns the bottom field in this field extension tower; if the tower
111   is flat, i.e., if there is no extension, then r itself is returned;
112   as a side-effect, the counter 'height' is filled with the height of
113   the extension tower (in case the tower is flat, 'height' is zero) */
114static coeffs nCoeff_bottom(const coeffs r, int &height)
115{
116  assume(r != NULL);
117  coeffs cf = r;
118  height = 0;
119  while (nCoeff_is_Extension(cf))
120  {
121    assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL);
122    cf = cf->extRing->cf;
123    height++;
124  }
125  return cf;
126}
127
128BOOLEAN ntIsZero(number a, const coeffs cf)
129{
130  ntTest(a);
131  return (IS0(a));
132}
133
134void ntDelete(number * a, const coeffs cf)
135{
136  fraction f = (fraction)(*a);
137  if (IS0(f)) return;
138  p_Delete(&NUM(f), ntRing);
139  if (!DENIS1(f)) p_Delete(&DEN(f), ntRing);
140  omFreeBin((ADDRESS)f, fractionObjectBin);
141  *a = NULL;
142}
143
144BOOLEAN ntEqual(number a, number b, const coeffs cf)
145{
146  ntTest(a); ntTest(b);
147 
148  /// simple tests
149  if (a == b) return TRUE;
150  if ((IS0(a)) && (!IS0(b))) return FALSE;
151  if ((IS0(b)) && (!IS0(a))) return FALSE;
152 
153  /// cheap test if gcd's have been cancelled in both numbers
154  fraction fa = (fraction)a;
155  fraction fb = (fraction)b;
156  if ((COM(fa) == 1) && (COM(fb) == 1))
157  {
158    poly f = p_Add_q(p_Copy(NUM(fa), ntRing),
159                     p_Neg(p_Copy(NUM(fb), ntRing), ntRing),
160                     ntRing);
161    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
162    if (DENIS1(fa) && DENIS1(fb))  return TRUE;
163    if (DENIS1(fa) && !DENIS1(fb)) return FALSE;
164    if (!DENIS1(fa) && DENIS1(fb)) return FALSE;
165    f = p_Add_q(p_Copy(DEN(fa), ntRing),
166                p_Neg(p_Copy(DEN(fb), ntRing), ntRing),
167                ntRing);
168    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
169    return TRUE;
170  }
171 
172  /* default: the more expensive multiplication test
173              a/b = c/d  <==>  a*d = b*c */
174  poly f = p_Copy(NUM(fa), ntRing);
175  if (!DENIS1(fb)) f = p_Mult_q(f, p_Copy(DEN(fb), ntRing), ntRing);
176  poly g = p_Copy(NUM(fb), ntRing);
177  if (!DENIS1(fa)) g = p_Mult_q(g, p_Copy(DEN(fa), ntRing), ntRing);
178  poly h = p_Add_q(f, p_Neg(g, ntRing), ntRing);
179  if (h == NULL) return TRUE;
180  else
181  {
182    p_Delete(&h, ntRing);
183    return FALSE;
184  }
185}
186
187number ntCopy(number a, const coeffs cf)
188{
189  ntTest(a);
190  if (IS0(a)) return NULL;
191  fraction f = (fraction)a;
192  poly g = p_Copy(NUM(f), ntRing);
193  poly h = NULL; if (!DENIS1(f)) h = p_Copy(DEN(f), ntRing);
194  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
195  NUM(result) = g;
196  DEN(result) = h;
197  COM(result) = COM(f);
198  return (number)result;
199}
200
201number ntGetNumerator(number &a, const coeffs cf)
202{
203  ntTest(a);
204  definiteGcdCancellation(a, cf, FALSE);
205  if (IS0(a)) return NULL;
206  fraction f = (fraction)a;
207  poly g = p_Copy(NUM(f), ntRing);
208  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
209  NUM(result) = g;
210  DEN(result) = NULL;
211  COM(result) = 0;
212  return (number)result;
213}
214
215number ntGetDenom(number &a, const coeffs cf)
216{
217  ntTest(a);
218  definiteGcdCancellation(a, cf, FALSE);
219  fraction f = (fraction)a;
220  poly g;
221  if (IS0(f) || DENIS1(f)) g = p_One(ntRing);
222  else g = p_Copy(DEN(f), ntRing);
223  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
224  NUM(result) = g;
225  DEN(result) = NULL;
226  COM(result) = 0;
227  return (number)result;
228}
229
230BOOLEAN ntIsOne(number a, const coeffs cf)
231{
232  ntTest(a);
233  definiteGcdCancellation(a, cf, FALSE);
234  fraction f = (fraction)a;
235  return DENIS1(f) && NUMIS1(f);
236}
237
238BOOLEAN ntIsMOne(number a, const coeffs cf)
239{
240  ntTest(a);
241  definiteGcdCancellation(a, cf, FALSE);
242  fraction f = (fraction)a;
243  if (!DENIS1(f)) return FALSE;
244  poly g = NUM(f);
245  if (!p_IsConstant(g, ntRing)) return FALSE;
246  return n_IsMOne(p_GetCoeff(g, ntRing), ntCoeffs);
247}
248
249/// this is in-place, modifies a
250number ntNeg(number a, const coeffs cf)
251{
252  ntTest(a);
253  if (!IS0(a))
254  {
255    fraction f = (fraction)a;
256    NUM(f) = p_Neg(NUM(f), ntRing);
257  }
258  return a;
259}
260
261number ntImPart(number a, const coeffs cf)
262{
263  ntTest(a);
264  return NULL;
265}
266
267number ntInit(int i, const coeffs cf)
268{
269  if (i == 0) return NULL;
270  else
271  {
272    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
273    NUM(result) = p_ISet(i, ntRing);
274    DEN(result) = NULL;
275    COM(result) = 0;
276    return (number)result;
277  }
278}
279
280int ntInt(number &a, const coeffs cf)
281{
282  ntTest(a);
283  if (IS0(a)) return 0;
284  definiteGcdCancellation(a, cf, FALSE);
285  fraction f = (fraction)a;
286  if (!DENIS1(f)) return 0;
287  if (!p_IsConstant(NUM(f), ntRing)) return 0;
288  return n_Int(p_GetCoeff(NUM(f), ntRing), ntCoeffs);
289}
290
291/* This method will only consider the numerators of a and b, without
292   cancelling gcd's before.
293   Moreover it may return TRUE only if one or both numerators
294   are zero or if their degrees are equal. Then TRUE is returned iff
295   coeff(numerator(a)) > coeff(numerator(b));
296   In all other cases, FALSE will be returned. */
297BOOLEAN ntGreater(number a, number b, const coeffs cf)
298{
299  ntTest(a); ntTest(b);
300  number aNumCoeff = NULL; int aNumDeg = 0;
301  number bNumCoeff = NULL; int bNumDeg = 0;
302  if (!IS0(a))
303  {
304    fraction fa = (fraction)a;
305    aNumDeg = p_Totaldegree(NUM(fa), ntRing);
306    aNumCoeff = p_GetCoeff(NUM(fa), ntRing);
307  }
308  if (!IS0(b))
309  {
310    fraction fb = (fraction)b;
311    bNumDeg = p_Totaldegree(NUM(fb), ntRing);
312    bNumCoeff = p_GetCoeff(NUM(fb), ntRing);
313  }
314  if (aNumDeg != bNumDeg) return FALSE;
315  else return n_Greater(aNumCoeff, bNumCoeff, ntCoeffs);
316}
317
318/* this method will only consider the numerator of a, without cancelling
319   the gcd before;
320   returns TRUE iff the leading coefficient of the numerator of a is > 0
321                    or the leading term of the numerator of a is not a
322                    constant */
323BOOLEAN ntGreaterZero(number a, const coeffs cf)
324{
325  ntTest(a);
326  if (IS0(a)) return FALSE;
327  fraction f = (fraction)a;
328  poly g = NUM(f);
329  return (n_GreaterZero(p_GetCoeff(g, ntRing), ntCoeffs) ||
330          (!p_LmIsConstant(g, ntRing)));
331}
332
333void ntCoeffWrite(const coeffs cf)
334{
335  PrintS("//   Coefficients live in the rational function field\n");
336  Print("//   K(");
337  for (int i = 0; i < rVar(ntRing); i++)
338  {
339    if (i > 0) PrintS(", ");
340    Print("%s", rRingVar(i, ntRing));
341  }
342  PrintS(") with\n");
343  PrintS("//   K: "); n_CoeffWrite(cf->extRing->cf);
344}
345
346number ntAdd(number a, number b, const coeffs cf)
347{
348  ntTest(a); ntTest(b);
349  if (IS0(a)) return ntCopy(b, cf);
350  if (IS0(b)) return ntCopy(a, cf);
351 
352  fraction fa = (fraction)a;
353  fraction fb = (fraction)b;
354 
355  poly g = p_Copy(NUM(fa), ntRing);
356  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
357  poly h = p_Copy(NUM(fb), ntRing);
358  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
359  g = p_Add_q(g, h, ntRing);
360 
361  if (g == NULL) return NULL;
362 
363  poly f;
364  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
365  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
366  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
367  else /* both denom's are != 1 */    f = p_Mult_q(p_Copy(DEN(fa), ntRing),
368                                                   p_Copy(DEN(fb), ntRing),
369                                                   ntRing);
370 
371  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
372  NUM(result) = g;
373  DEN(result) = f;
374  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
375  heuristicGcdCancellation((number)result, cf);
376  return (number)result;
377}
378
379number ntSub(number a, number b, const coeffs cf)
380{
381  ntTest(a); ntTest(b);
382  if (IS0(a)) return ntNeg(ntCopy(b, cf), cf);
383  if (IS0(b)) return ntCopy(a, cf);
384 
385  fraction fa = (fraction)a;
386  fraction fb = (fraction)b;
387 
388  poly g = p_Copy(NUM(fa), ntRing);
389  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
390  poly h = p_Copy(NUM(fb), ntRing);
391  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
392  g = p_Add_q(g, p_Neg(h, ntRing), ntRing);
393 
394  if (g == NULL) return NULL;
395 
396  poly f;
397  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
398  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
399  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
400  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
401                                                   p_Copy(DEN(fb), ntRing),
402                                                   ntRing);
403 
404  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
405  NUM(result) = g;
406  DEN(result) = f;
407  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
408  heuristicGcdCancellation((number)result, cf);
409  return (number)result;
410}
411
412number ntMult(number a, number b, const coeffs cf)
413{
414  ntTest(a); ntTest(b);
415  if (IS0(a) || IS0(b)) return NULL;
416 
417  fraction fa = (fraction)a;
418  fraction fb = (fraction)b;
419 
420  poly g = p_Copy(NUM(fa), ntRing);
421  poly h = p_Copy(NUM(fb), ntRing);
422  g = p_Mult_q(g, h, ntRing);
423 
424  if (g == NULL) return NULL;   /* may happen due to zero divisors */
425 
426  poly f;
427  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
428  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
429  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
430  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
431                                                   p_Copy(DEN(fb), ntRing),
432                                                   ntRing);
433 
434  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
435  NUM(result) = g;
436  DEN(result) = f;
437  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
438  heuristicGcdCancellation((number)result, cf);
439  return (number)result;
440}
441
442number ntDiv(number a, number b, const coeffs cf)
443{
444  ntTest(a); ntTest(b);
445  if (IS0(a)) return NULL;
446  if (IS0(b)) WerrorS(nDivBy0);
447 
448  fraction fa = (fraction)a;
449  fraction fb = (fraction)b;
450 
451  poly g = p_Copy(NUM(fa), ntRing);
452  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
453 
454  if (g == NULL) return NULL;   /* may happen due to zero divisors */
455 
456  poly f = p_Copy(NUM(fb), ntRing);
457  if (!DENIS1(fa)) f = p_Mult_q(f, p_Copy(DEN(fa), ntRing), ntRing);
458 
459  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
460  NUM(result) = g;
461  DEN(result) = f;
462  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
463  heuristicGcdCancellation((number)result, cf);
464  return (number)result;
465}
466
467/* 0^0 = 0;
468   for |exp| <= 7 compute power by a simple multiplication loop;
469   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
470      p^13 = p^1 * p^4 * p^8, where we utilise that
471      p^(2^(k+1)) = p^(2^k) * p^(2^k);
472   intermediate cancellation is controlled by the in-place method
473   heuristicGcdCancellation; see there.
474*/
475void ntPower(number a, int exp, number *b, const coeffs cf)
476{
477  ntTest(a);
478 
479  /* special cases first */
480  if (IS0(a))
481  {
482    if (exp >= 0) *b = NULL;
483    else          WerrorS(nDivBy0);
484  }
485  else if (exp ==  0) *b = ntInit(1, cf);
486  else if (exp ==  1) *b = ntCopy(a, cf);
487  else if (exp == -1) *b = ntInvers(a, cf);
488 
489  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
490 
491  /* now compute a^expAbs */
492  number pow; number t;
493  if (expAbs <= 7)
494  {
495    pow = ntCopy(a, cf);
496    for (int i = 2; i <= expAbs; i++)
497    {
498      t = ntMult(pow, a, cf);
499      ntDelete(&pow, cf);
500      pow = t;
501      heuristicGcdCancellation(pow, cf);
502    }
503  }
504  else
505  {
506    pow = ntInit(1, cf);
507    number factor = ntCopy(a, cf);
508    while (expAbs != 0)
509    {
510      if (expAbs & 1)
511      {
512        t = ntMult(pow, factor, cf);
513        ntDelete(&pow, cf);
514        pow = t;
515        heuristicGcdCancellation(pow, cf);
516      }
517      expAbs = expAbs / 2;
518      if (expAbs != 0)
519      {
520        t = ntMult(factor, factor, cf);
521        ntDelete(&factor, cf);
522        factor = t;
523        heuristicGcdCancellation(factor, cf);
524      }
525    }
526    ntDelete(&factor, cf);
527  }
528 
529  /* invert if original exponent was negative */
530  if (exp < 0)
531  {
532    t = ntInvers(pow, cf);
533    ntDelete(&pow, cf);
534    pow = t;
535  }
536  *b = pow;
537}
538
539/* assumes that cf represents the rationals, i.e. Q, and will only
540   be called in that case;
541   assumes furthermore that f != NULL and that the denominator of f != 1;
542   generally speaking, this method removes denominators in the rational
543   coefficients of the numerator and denominator of 'a';
544   more concretely, the following normalizations will be performed,
545   where t^alpha denotes a monomial in the transcendental variables t_k
546   (1) if 'a' is of the form
547          (sum_alpha a_alpha/b_alpha * t^alpha)
548          -------------------------------------
549            (sum_beta c_beta/d_beta * t^beta)
550       with integers a_alpha, b_alpha, c_beta, d_beta, then both the
551       numerator and the denominator will be multiplied by the LCM of
552       the b_alpha's and the d_beta's (if this LCM is != 1),
553   (2) if 'a' is - e.g. after having performed step (1) - of the form
554          (sum_alpha a_alpha * t^alpha)
555          -----------------------------
556            (sum_beta c_beta * t^beta)
557       with integers a_alpha, c_beta, and with a non-constant denominator,
558       then both the numerator and the denominator will be divided by the
559       GCD of the a_alpha's and the c_beta's (if this GCD is != 1),
560   (3) if 'a' is - e.g. after having performed steps (1) and (2) - of the
561       form
562          (sum_alpha a_alpha * t^alpha)
563          -----------------------------
564                        c
565       with integers a_alpha, and c != 1, then 'a' will be replaced by
566       (sum_alpha a_alpha/c * t^alpha);
567   this procedure does not alter COM(f) (this has to be done by the
568   calling procedure);
569   modifies f */
570void handleNestedFractionsOverQ(fraction f, const coeffs cf)
571{
572  assume(nCoeff_is_Q(ntCoeffs));
573  assume(!IS0(f));
574  assume(!DENIS1(f));
575 
576  if (!p_IsConstant(DEN(f), ntRing))
577  { /* step (1); see documentation of this procedure above */
578    p_Normalize(NUM(f), ntRing);
579    p_Normalize(DEN(f), ntRing);
580    number lcmOfDenominators = n_Init(1, ntCoeffs);
581    number c; number tmp;
582    poly p = NUM(f);
583    /* careful when using n_Lcm!!! It computes the lcm of the numerator
584       of the 1st argument and the denominator of the 2nd!!! */
585    while (p != NULL)
586    {
587      c = p_GetCoeff(p, ntRing);
588      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
589      n_Delete(&lcmOfDenominators, ntCoeffs);
590      lcmOfDenominators = tmp;
591      pIter(p);
592    }
593    p = DEN(f);
594    while (p != NULL)
595    {
596      c = p_GetCoeff(p, ntRing);
597      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
598      n_Delete(&lcmOfDenominators, ntCoeffs);
599      lcmOfDenominators = tmp;
600      pIter(p);
601    }
602    if (!n_IsOne(lcmOfDenominators, ntCoeffs))
603    { /* multiply NUM(f) and DEN(f) with lcmOfDenominators */
604      NUM(f) = p_Mult_nn(NUM(f), lcmOfDenominators, ntRing);
605      p_Normalize(NUM(f), ntRing);
606      DEN(f) = p_Mult_nn(DEN(f), lcmOfDenominators, ntRing);
607      p_Normalize(DEN(f), ntRing);
608    }
609    n_Delete(&lcmOfDenominators, ntCoeffs);
610    if (!p_IsConstant(DEN(f), ntRing))
611    { /* step (2); see documentation of this procedure above */
612      p = NUM(f);
613      number gcdOfCoefficients = n_Copy(p_GetCoeff(p, ntRing), ntCoeffs);
614      pIter(p);
615      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
616      {
617        c = p_GetCoeff(p, ntRing);
618        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
619        n_Delete(&gcdOfCoefficients, ntCoeffs);
620        gcdOfCoefficients = tmp;
621        pIter(p);
622      }
623      p = DEN(f);
624      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
625      {
626        c = p_GetCoeff(p, ntRing);
627        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
628        n_Delete(&gcdOfCoefficients, ntCoeffs);
629        gcdOfCoefficients = tmp;
630        pIter(p);
631      }
632      if (!n_IsOne(gcdOfCoefficients, ntCoeffs))
633      { /* divide NUM(f) and DEN(f) by gcdOfCoefficients */
634        number inverseOfGcdOfCoefficients = n_Invers(gcdOfCoefficients,
635                                                     ntCoeffs);
636        NUM(f) = p_Mult_nn(NUM(f), inverseOfGcdOfCoefficients, ntRing);
637        p_Normalize(NUM(f), ntRing);
638        DEN(f) = p_Mult_nn(DEN(f), inverseOfGcdOfCoefficients, ntRing);
639        p_Normalize(DEN(f), ntRing);
640        n_Delete(&inverseOfGcdOfCoefficients, ntCoeffs);
641      }
642      n_Delete(&gcdOfCoefficients, ntCoeffs);
643    }
644  }
645  if (p_IsConstant(DEN(f), ntRing) &&
646      (!n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs)))
647  { /* step (3); see documentation of this procedure above */
648    number inverseOfDen = n_Invers(p_GetCoeff(DEN(f), ntRing), ntCoeffs);
649    NUM(f) = p_Mult_nn(NUM(f), inverseOfDen, ntRing);
650    n_Delete(&inverseOfDen, ntCoeffs);
651    p_Delete(&DEN(f), ntRing);
652    DEN(f) = NULL;
653  }
654 
655  /* Now, due to the above computations, DEN(f) may have become the
656     1-polynomial which needs to be represented by NULL: */
657  if ((DEN(f) != NULL) &&
658      p_IsConstant(DEN(f), ntRing) &&
659      n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
660  {
661    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
662  }
663}
664
665/* modifies a */
666void heuristicGcdCancellation(number a, const coeffs cf)
667{
668  ntTest(a);
669  if (IS0(a)) return;
670 
671  fraction f = (fraction)a;
672  if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
673 
674  /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
675  if (p_EqualPolys(NUM(f), DEN(f), ntRing))
676  { /* numerator and denominator are both != 1 */
677    p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
678    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
679    COM(f) = 0;
680    return;
681  }
682 
683  if (COM(f) <= BOUND_COMPLEXITY) return;
684  else definiteGcdCancellation(a, cf, TRUE);
685}
686
687/* modifies a */
688void definiteGcdCancellation(number a, const coeffs cf,
689                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed)
690{
691  ntTest(a);
692 
693  fraction f = (fraction)a;
694 
695  if (!simpleTestsHaveAlreadyBeenPerformed)
696  {
697    if (IS0(a)) return;
698    if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
699 
700    /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
701    if (p_EqualPolys(NUM(f), DEN(f), ntRing))
702    { /* numerator and denominator are both != 1 */
703      p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
704      p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
705      COM(f) = 0;
706      return;
707    }
708  }
709 
710#ifdef HAVE_FACTORY 
711  /* singclap_gcd destroys its arguments; we hence need copies: */
712  poly pNum = p_Copy(NUM(f), ntRing);
713  poly pDen = p_Copy(DEN(f), ntRing);
714 
715  /* Note that, over Q, singclap_gcd will remove the denominators in all
716     rational coefficients of pNum and pDen, before starting to compute
717     the gcd. Thus, we do not need to ensure that the coefficients of
718     pNum and pDen live in Z; they may well be elements of Q\Z. */
719  poly pGcd = singclap_gcd(pNum, pDen, cf->extRing);
720  if (p_IsConstant(pGcd, ntRing) &&
721      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
722  { /* gcd = 1; nothing to cancel;
723       Suppose the given rational function field is over Q. Although the
724       gcd is 1, we may have produced fractional coefficients in NUM(f),
725       DEN(f), or both, due to previous arithmetics. The next call will
726       remove those nested fractions, in case there are any. */
727    if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
728  }
729  else
730  { /* We divide both NUM(f) and DEN(f) by the gcd which is known
731       to be != 1. */
732    poly newNum = singclap_pdivide(NUM(f), pGcd, ntRing);
733    p_Delete(&NUM(f), ntRing);
734    NUM(f) = newNum;
735    poly newDen = singclap_pdivide(DEN(f), pGcd, ntRing);
736    p_Delete(&DEN(f), ntRing);
737    DEN(f) = newDen;
738    if (p_IsConstant(DEN(f), ntRing) &&
739        n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
740    {
741      /* DEN(f) = 1 needs to be represented by NULL! */
742      p_Delete(&DEN(f), ntRing);
743      newDen = NULL;
744    }
745    else
746    { /* Note that over Q, by cancelling the gcd, we may have produced
747         fractional coefficients in NUM(f), DEN(f), or both. The next
748         call will remove those nested fractions, in case there are
749         any. */
750      if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
751    }
752  }
753  COM(f) = 0;
754  p_Delete(&pGcd, ntRing);
755#endif /* HAVE_FACTORY */
756}
757
758/* modifies a */
759void ntWrite(number &a, const coeffs cf)
760{
761  ntTest(a);
762  definiteGcdCancellation(a, cf, FALSE);
763  if (IS0(a))
764    StringAppendS("0");
765  else
766  {
767    fraction f = (fraction)a;
768    BOOLEAN useBrackets = (!p_IsConstant(NUM(f), ntRing)) ||
769                          (!n_GreaterZero(p_GetCoeff(NUM(f), ntRing),
770                                          ntCoeffs));
771    if (useBrackets) StringAppendS("(");
772    p_String0(NUM(f), ntRing, ntRing);
773    if (useBrackets) StringAppendS(")");
774    if (!DENIS1(f))
775    {
776      StringAppendS("/");
777      useBrackets = (!p_IsConstant(DEN(f), ntRing)) ||
778                    (!n_GreaterZero(p_GetCoeff(DEN(f), ntRing), ntCoeffs));
779      if (useBrackets) StringAppendS("(");
780      p_String0(DEN(f), ntRing, ntRing);
781      if (useBrackets) StringAppendS(")");
782    }
783  }
784}
785
786const char * ntRead(const char *s, number *a, const coeffs cf)
787{
788  poly p;
789  const char * result = p_Read(s, p, ntRing);
790  if (p == NULL) { *a = NULL; return result; }
791  else
792  {
793    fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
794    NUM(f) = p;
795    DEN(f) = NULL;
796    COM(f) = 0;
797    *a = (number)f;
798    return result;
799  }
800}
801
802/* expects *param to be castable to TransExtInfo */
803static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
804{
805  if (ntID != n) return FALSE;
806  TransExtInfo *e = (TransExtInfo *)param;
807  /* for rational function fields we expect the underlying
808     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
809     this expectation is based on the assumption that we have properly
810     registered cf and perform reference counting rather than creating
811     multiple copies of the same coefficient field/domain/ring */
812  return (ntRing == e->r);
813}
814
815number ntLcm(number a, number b, const coeffs cf)
816{
817  ntTest(a); ntTest(b);
818  /* TO BE IMPLEMENTED!
819     for the time, we simply return NULL, representing the number zero */
820  Print("// TO BE IMPLEMENTED: transext.cc:ntLcm\n");
821  return NULL;
822}
823
824number ntGcd(number a, number b, const coeffs cf)
825{
826  ntTest(a); ntTest(b);
827  /* TO BE IMPLEMENTED!
828     for the time, we simply return NULL, representing the number zero */
829  Print("// TO BE IMPLEMENTED: transext.cc:ntGcd\n");
830  return NULL;
831}
832
833int ntSize(number a, const coeffs cf)
834{
835  ntTest(a);
836  if (IS0(a)) return -1;
837  /* this has been taken from the old implementation of field extensions,
838     where we computed the sum of the degrees and the numbers of terms in
839     the numerator and denominator of a; so we leave it at that, for the
840     time being */
841  fraction f = (fraction)a;
842  poly p = NUM(f);
843  int noOfTerms = 0;
844  int numDegree = 0;
845  while (p != NULL)
846  {
847    noOfTerms++;
848    int d = 0;
849    for (int i = 1; i <= rVar(ntRing); i++)
850      d += p_GetExp(p, i, ntRing);
851    if (d > numDegree) numDegree = d;
852    pIter(p);
853  }
854  int denDegree = 0;
855  if (!DENIS1(f))
856  {
857    p = DEN(f);
858    while (p != NULL)
859    {
860      noOfTerms++;
861      int d = 0;
862      for (int i = 1; i <= rVar(ntRing); i++)
863        d += p_GetExp(p, i, ntRing);
864      if (d > denDegree) denDegree = d;
865      pIter(p);
866    }
867  }
868  return numDegree + denDegree + noOfTerms;
869}
870
871number ntInvers(number a, const coeffs cf)
872{
873  ntTest(a);
874  if (IS0(a)) WerrorS(nDivBy0);
875  fraction f = (fraction)a;
876  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
877  poly g;
878  if (DENIS1(f)) g = p_One(ntRing);
879  else           g = p_Copy(DEN(f), ntRing);
880  NUM(result) = g;
881  DEN(result) = p_Copy(NUM(f), ntRing);
882  COM(result) = COM(f);
883  return (number)result;
884}
885
886/* assumes that src = Q, dst = Q(t_1, ..., t_s) */
887number ntMap00(number a, const coeffs src, const coeffs dst)
888{
889  if (n_IsZero(a, src)) return NULL;
890  assume(src == dst->extRing->cf);
891  poly p = p_One(dst->extRing);
892  p_SetCoeff(p, ntCopy(a, src), dst->extRing);
893  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
894  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
895  return (number)f;
896}
897
898/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
899number ntMapP0(number a, const coeffs src, const coeffs dst)
900{
901  if (n_IsZero(a, src)) return NULL;
902  /* mapping via intermediate int: */
903  int n = n_Int(a, src);
904  number q = n_Init(n, dst->extRing->cf);
905  poly p;
906  if (n_IsZero(q, dst->extRing->cf))
907  {
908    n_Delete(&q, dst->extRing->cf);
909    return NULL;
910  }
911  p = p_One(dst->extRing);
912  p_SetCoeff(p, q, dst->extRing);
913  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
914  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
915  return (number)f;
916}
917
918/* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
919                       src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
920number ntCopyMap(number a, const coeffs src, const coeffs dst)
921{
922  return ntCopy(a, dst);
923}
924
925/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
926number ntMap0P(number a, const coeffs src, const coeffs dst)
927{
928  if (n_IsZero(a, src)) return NULL;
929  int p = rChar(dst->extRing);
930  int n = nlModP(a, p, src);
931  number q = n_Init(n, dst->extRing->cf);
932  poly g;
933  if (n_IsZero(q, dst->extRing->cf))
934  {
935    n_Delete(&q, dst->extRing->cf);
936    return NULL;
937  }
938  g = p_One(dst->extRing);
939  p_SetCoeff(g, q, dst->extRing);
940  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
941  NUM(f) = g; DEN(f) = NULL; COM(f) = 0;
942  return (number)f;
943}
944
945/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
946number ntMapPP(number a, const coeffs src, const coeffs dst)
947{
948  if (n_IsZero(a, src)) return NULL;
949  assume(src == dst->extRing->cf);
950  poly p = p_One(dst->extRing);
951  p_SetCoeff(p, ntCopy(a, src), dst->extRing);
952  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
953  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
954  return (number)f;
955}
956
957/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
958number ntMapUP(number a, const coeffs src, const coeffs dst)
959{
960  if (n_IsZero(a, src)) return NULL;
961  /* mapping via intermediate int: */
962  int n = n_Int(a, src);
963  number q = n_Init(n, dst->extRing->cf);
964  poly p;
965  if (n_IsZero(q, dst->extRing->cf))
966  {
967    n_Delete(&q, dst->extRing->cf);
968    return NULL;
969  }
970  p = p_One(dst->extRing);
971  p_SetCoeff(p, q, dst->extRing);
972  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
973  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
974  return (number)f;
975}
976
977nMapFunc ntSetMap(const coeffs src, const coeffs dst)
978{
979  /* dst is expected to be a rational function field */
980  assume(getCoeffType(dst) == ntID);
981 
982  int h = 0; /* the height of the extension tower given by dst */
983  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
984 
985  /* for the time being, we only provide maps if h = 1 and if b is Q or
986     some field Z/pZ: */
987  if (h != 1) return NULL;
988  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
989 
990  /* Let T denote the sequence of transcendental extension variables, i.e.,
991     K[t_1, ..., t_s] =: K[T];
992     Let moreover, for any such sequence T, T' denote any subsequence of T
993     of the form t_1, ..., t_w with w <= s. */
994 
995  if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
996    return ntMap00;                                      /// Q       -->  Q(T)
997 
998  if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
999    return ntMapP0;                                      /// Z/p     -->  Q(T)
1000 
1001  if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
1002    return ntMap0P;                                      /// Q       --> Z/p(T)
1003 
1004  if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
1005  {
1006    if (src->ch == dst->ch) return ntMapPP;              /// Z/p     --> Z/p(T)
1007    else return ntMapUP;                                 /// Z/u     --> Z/p(T)
1008  }
1009 
1010  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
1011  if (h != 1) return NULL;
1012  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
1013 
1014  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
1015  {
1016    if (rVar(src->extRing) > rVar(dst->extRing)) return NULL;
1017    for (int i = 0; i < rVar(src->extRing); i++)
1018      if (strcmp(rRingVar(i, src->extRing),
1019                 rRingVar(i, dst->extRing)) != 0) return NULL;
1020      return ntCopyMap;                                  /// Q(T')   --> Q(T)
1021  }
1022 
1023  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
1024  {
1025    if (rVar(src->extRing) > rVar(dst->extRing)) return NULL;
1026    for (int i = 0; i < rVar(src->extRing); i++)
1027      if (strcmp(rRingVar(i, src->extRing),
1028                 rRingVar(i, dst->extRing)) != 0) return NULL;
1029      return ntCopyMap;                                  /// Z/p(T') --> Z/p(T)
1030  }
1031 
1032  return NULL;                                           /// default
1033}
1034
1035BOOLEAN ntInitChar(coeffs cf, void * infoStruct)
1036{ 
1037  TransExtInfo *e = (TransExtInfo *)infoStruct;
1038  /// first check whether cf->extRing != NULL and delete old ring???
1039  cf->extRing           = e->r;
1040  cf->extRing->minideal = NULL;
1041
1042  assume(cf->extRing                != NULL);      // extRing;
1043  assume(cf->extRing->cf            != NULL);      // extRing->cf;
1044  assume(getCoeffType(cf) == ntID);                // coeff type;
1045 
1046  /* propagate characteristic up so that it becomes
1047     directly accessible in cf: */
1048  cf->ch = cf->extRing->cf->ch;
1049 
1050  cf->cfGreaterZero  = ntGreaterZero;
1051  cf->cfGreater      = ntGreater;
1052  cf->cfEqual        = ntEqual;
1053  cf->cfIsZero       = ntIsZero;
1054  cf->cfIsOne        = ntIsOne;
1055  cf->cfIsMOne       = ntIsMOne;
1056  cf->cfInit         = ntInit;
1057  cf->cfInt          = ntInt;
1058  cf->cfNeg          = ntNeg;
1059  cf->cfAdd          = ntAdd;
1060  cf->cfSub          = ntSub;
1061  cf->cfMult         = ntMult;
1062  cf->cfDiv          = ntDiv;
1063  cf->cfExactDiv     = ntDiv;
1064  cf->cfPower        = ntPower;
1065  cf->cfCopy         = ntCopy;
1066  cf->cfWrite        = ntWrite;
1067  cf->cfRead         = ntRead;
1068  cf->cfDelete       = ntDelete;
1069  cf->cfSetMap       = ntSetMap;
1070  cf->cfGetDenom     = ntGetDenom;
1071  cf->cfGetNumerator = ntGetNumerator;
1072  cf->cfRePart       = ntCopy;
1073  cf->cfImPart       = ntImPart;
1074  cf->cfCoeffWrite   = ntCoeffWrite;
1075  cf->cfDBTest       = ntDBTest;
1076  cf->cfGcd          = ntGcd;
1077  cf->cfLcm          = ntLcm;
1078  cf->cfSize         = ntSize;
1079  cf->nCoeffIsEqual  = ntCoeffIsEqual;
1080  cf->cfInvers       = ntInvers;
1081  cf->cfIntDiv       = ntDiv;
1082 
1083#ifndef HAVE_FACTORY
1084  PrintS("// Warning: The 'factory' module is not available.\n");
1085  PrintS("//          Hence gcd's cannot be cancelled in any\n");
1086  PrintS("//          computed fraction!\n");
1087#endif
1088 
1089  return FALSE;
1090}
Note: See TracBrowser for help on using the repository browser.