source: git/libpolys/polys/ext_fields/transext.cc @ 2d3091c

spielwiese
Last change on this file since 2d3091c was 2d3091c, checked in by Frank Seelisch <seelisch@…>, 13 years ago
test for trac ticket #308 with new implementation (result is instantaneous)
  • Property mode set to 100644
File size: 27.1 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 in the polynomial
10*           ring K[t_1, .., t_s] represented by cf->extring.
11*/
12
13#include "config.h"
14#include <misc/auxiliary.h>
15
16#include <omalloc/omalloc.h>
17
18#include <reporter/reporter.h>
19
20#include <coeffs/coeffs.h>
21#include <coeffs/numbers.h>
22#include <coeffs/longrat.h>
23
24#include <polys/monomials/ring.h>
25#include <polys/monomials/p_polys.h>
26#include <polys/simpleideals.h>
27
28#ifdef HAVE_FACTORY
29#include <polys/clapsing.h>
30#endif
31
32#include <polys/ext_fields/transext.h>
33
34/// our type has been defined as a macro in transext.h
35/// and is accessible by 'ntID'
36
37extern omBin fractionObjectBin = omGetSpecBin(sizeof(fractionObject));
38
39/// forward declarations
40BOOLEAN  ntGreaterZero(number a, const coeffs cf); 
41BOOLEAN  ntGreater(number a, number b, const coeffs cf);
42BOOLEAN  ntEqual(number a, number b, const coeffs cf);
43BOOLEAN  ntIsOne(number a, const coeffs cf);
44BOOLEAN  ntIsMOne(number a, const coeffs cf);
45BOOLEAN  ntIsZero(number a, const coeffs cf);
46number   ntInit(int i, const coeffs cf);
47int      ntInt(number &a, const coeffs cf);
48number   ntNeg(number a, const coeffs cf);
49number   ntInvers(number a, const coeffs cf);
50number   ntPar(int i, const coeffs cf);
51number   ntAdd(number a, number b, const coeffs cf);
52number   ntSub(number a, number b, const coeffs cf);
53number   ntMult(number a, number b, const coeffs cf);
54number   ntDiv(number a, number b, const coeffs cf);
55void     ntPower(number a, int exp, number *b, const coeffs cf);
56number   ntCopy(number a, const coeffs cf);
57void     ntWrite(number &a, const coeffs cf);
58number   ntRePart(number a, const coeffs cf);
59number   ntImPart(number a, const coeffs cf);
60number   ntGetDenom(number &a, const coeffs cf);
61number   ntGetNumerator(number &a, const coeffs cf);
62number   ntGcd(number a, number b, const coeffs cf);
63number   ntLcm(number a, number b, const coeffs cf);
64int      ntSize(number a, const coeffs cf);
65void     ntDelete(number * a, const coeffs cf);
66void     ntCoeffWrite(const coeffs cf);
67number   ntIntDiv(number a, number b, const coeffs cf);
68const char * ntRead(const char *s, number *a, const coeffs cf);
69static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param);
70
71void heuristicGcdCancellation(number a, const coeffs cf);
72void definiteGcdCancellation(number a, const coeffs cf,
73                             BOOLEAN skipSimpleTests);
74
75#ifdef LDEBUG
76BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs cf)
77{
78  assume(getCoeffType(cf) == ntID);
79  fraction t = (fraction)a;
80  if (IS0(t)) return TRUE;
81  assume(NUM(t) != NULL);   /**< t != 0 ==> numerator(t) != 0 */
82  p_Test(NUM(t), ntRing);
83  if (!DENIS1(t)) p_Test(DEN(t), ntRing);
84  return TRUE;
85}
86#endif
87
88/* returns the bottom field in this field extension tower; if the tower
89   is flat, i.e., if there is no extension, then r itself is returned;
90   as a side-effect, the counter 'height' is filled with the height of
91   the extension tower (in case the tower is flat, 'height' is zero) */
92static coeffs nCoeff_bottom(const coeffs r, int &height)
93{
94  assume(r != NULL);
95  coeffs cf = r;
96  height = 0;
97  while (nCoeff_is_Extension(cf))
98  {
99    assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL);
100    cf = cf->extRing->cf;
101    height++;
102  }
103  return cf;
104}
105
106BOOLEAN ntIsZero(number a, const coeffs cf)
107{
108  ntTest(a);
109  return (IS0(a));
110}
111
112void ntDelete(number * a, const coeffs cf)
113{
114  fraction f = (fraction)(*a);
115  if (IS0(f)) return;
116  p_Delete(&NUM(f), ntRing);
117  if (!DENIS1(f)) p_Delete(&DEN(f), ntRing);
118  omFreeBin((ADDRESS)f, fractionObjectBin);
119  *a = NULL;
120}
121
122BOOLEAN ntEqual(number a, number b, const coeffs cf)
123{
124  ntTest(a); ntTest(b);
125 
126  /// simple tests
127  if (a == b) return TRUE;
128  if ((IS0(a)) && (!IS0(b))) return FALSE;
129  if ((IS0(b)) && (!IS0(a))) return FALSE;
130 
131  /// cheap test if gcd's have been cancelled in both numbers
132  fraction fa = (fraction)a;
133  fraction fb = (fraction)b;
134  if ((COM(fa) == 1) && (COM(fb) == 1))
135  {
136    poly f = p_Add_q(p_Copy(NUM(fa), ntRing),
137                     p_Neg(p_Copy(NUM(fb), ntRing), ntRing),
138                     ntRing);
139    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
140    if (DENIS1(fa) && DENIS1(fb))  return TRUE;
141    if (DENIS1(fa) && !DENIS1(fb)) return FALSE;
142    if (!DENIS1(fa) && DENIS1(fb)) return FALSE;
143    f = p_Add_q(p_Copy(DEN(fa), ntRing),
144                p_Neg(p_Copy(DEN(fb), ntRing), ntRing),
145                ntRing);
146    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
147    return TRUE;
148  }
149 
150  /* default: the more expensive multiplication test
151              a/b = c/d  <==>  a*d = b*c */
152  poly f = p_Copy(NUM(fa), ntRing);
153  if (!DENIS1(fb)) f = p_Mult_q(f, p_Copy(DEN(fb), ntRing), ntRing);
154  poly g = p_Copy(NUM(fb), ntRing);
155  if (!DENIS1(fa)) g = p_Mult_q(g, p_Copy(DEN(fa), ntRing), ntRing);
156  poly h = p_Add_q(f, p_Neg(g, ntRing), ntRing);
157  if (h == NULL) return TRUE;
158  else
159  {
160    p_Delete(&h, ntRing);
161    return FALSE;
162  }
163}
164
165number ntCopy(number a, const coeffs cf)
166{
167  ntTest(a);
168  if (IS0(a)) return NULL;
169  fraction f = (fraction)a;
170  poly g = p_Copy(NUM(f), ntRing);
171  poly h = NULL; if (!DENIS1(f)) h = p_Copy(DEN(f), ntRing);
172  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
173  NUM(result) = g;
174  DEN(result) = h;
175  COM(result) = COM(f);
176  return (number)result;
177}
178
179number ntGetNumerator(number &a, const coeffs cf)
180{
181  ntTest(a);
182  definiteGcdCancellation(a, cf, FALSE);
183  if (IS0(a)) return NULL;
184  fraction f = (fraction)a;
185  poly g = p_Copy(NUM(f), ntRing);
186  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
187  NUM(result) = g;
188  DEN(result) = NULL;
189  COM(result) = 0;
190  return (number)result;
191}
192
193number ntGetDenom(number &a, const coeffs cf)
194{
195  ntTest(a);
196  definiteGcdCancellation(a, cf, FALSE);
197  fraction f = (fraction)a;
198  poly g;
199  if (IS0(f) || DENIS1(f)) g = p_One(ntRing);
200  else g = p_Copy(DEN(f), ntRing);
201  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
202  NUM(result) = g;
203  DEN(result) = NULL;
204  COM(result) = 0;
205  return (number)result;
206}
207
208BOOLEAN ntIsOne(number a, const coeffs cf)
209{
210  ntTest(a);
211  definiteGcdCancellation(a, cf, FALSE);
212  fraction f = (fraction)a;
213  return DENIS1(f) && NUMIS1(f);
214}
215
216BOOLEAN ntIsMOne(number a, const coeffs cf)
217{
218  ntTest(a);
219  definiteGcdCancellation(a, cf, FALSE);
220  fraction f = (fraction)a;
221  if (!DENIS1(f)) return FALSE;
222  poly g = NUM(f);
223  if (!p_IsConstant(g, ntRing)) return FALSE;
224  return n_IsMOne(p_GetCoeff(g, ntRing), ntCoeffs);
225}
226
227/// this is in-place, modifies a
228number ntNeg(number a, const coeffs cf)
229{
230  ntTest(a);
231  if (!IS0(a))
232  {
233    fraction f = (fraction)a;
234    NUM(f) = p_Neg(NUM(f), ntRing);
235  }
236  return a;
237}
238
239number ntImPart(number a, const coeffs cf)
240{
241  ntTest(a);
242  return NULL;
243}
244
245number ntInit(int i, const coeffs cf)
246{
247  if (i == 0) return NULL;
248  else
249  {
250    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
251    NUM(result) = p_ISet(i, ntRing);
252    DEN(result) = NULL;
253    COM(result) = 0;
254    return (number)result;
255  }
256}
257
258int ntInt(number &a, const coeffs cf)
259{
260  ntTest(a);
261  if (IS0(a)) return 0;
262  definiteGcdCancellation(a, cf, FALSE);
263  fraction f = (fraction)a;
264  if (!DENIS1(f)) return 0;
265  if (!p_IsConstant(NUM(f), ntRing)) return 0;
266  return n_Int(p_GetCoeff(NUM(f), ntRing), ntCoeffs);
267}
268
269/* This method will only consider the numerators of a and b, without
270   cancelling gcd's before.
271   Moreover it may return TRUE only if one or both numerators
272   are zero or if their degrees are equal. Then TRUE is returned iff
273   coeff(numerator(a)) > coeff(numerator(b));
274   In all other cases, FALSE will be returned. */
275BOOLEAN ntGreater(number a, number b, const coeffs cf)
276{
277  ntTest(a); ntTest(b);
278  number aNumCoeff = NULL; int aNumDeg = 0;
279  number bNumCoeff = NULL; int bNumDeg = 0;
280  if (!IS0(a))
281  {
282    fraction fa = (fraction)a;
283    aNumDeg = p_Totaldegree(NUM(fa), ntRing);
284    aNumCoeff = p_GetCoeff(NUM(fa), ntRing);
285  }
286  if (!IS0(b))
287  {
288    fraction fb = (fraction)b;
289    bNumDeg = p_Totaldegree(NUM(fb), ntRing);
290    bNumCoeff = p_GetCoeff(NUM(fb), ntRing);
291  }
292  if (aNumDeg != bNumDeg) return FALSE;
293  else return n_Greater(aNumCoeff, bNumCoeff, ntCoeffs);
294}
295
296/* this method will only consider the numerator of a, without cancelling
297   the gcd before;
298   returns TRUE iff the leading coefficient of the numerator of a is > 0
299                    or the leading term of the numerator of a is not a
300                    constant */
301BOOLEAN ntGreaterZero(number a, const coeffs cf)
302{
303  ntTest(a);
304  if (IS0(a)) return FALSE;
305  fraction f = (fraction)a;
306  poly g = NUM(f);
307  return (n_GreaterZero(p_GetCoeff(g, ntRing), ntCoeffs) ||
308          (!p_LmIsConstant(g, ntRing)));
309}
310
311void ntCoeffWrite(const coeffs cf)
312{
313  PrintS("//   Coefficients live in the rational function field\n");
314  Print("//   K(");
315  for (int i = 0; i < rVar(ntRing); i++)
316  {
317    if (i > 0) PrintS(", ");
318    Print("%s", rRingVar(i, ntRing));
319  }
320  PrintS(") with\n");
321  PrintS("//   K: "); n_CoeffWrite(cf->extRing->cf);
322}
323
324/* the i-th parameter */
325number ntPar(int i, const coeffs cf)
326{
327  assume((1 <= i) && (i <= rVar(ntRing)));
328  poly p = p_ISet(1, ntRing);
329  p_SetExp(p, i, 1, ntRing);
330  p_Setm(p, ntRing);
331  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
332  NUM(result) = p;
333  DEN(result) = NULL;
334  COM(result) = 0;
335  return (number)result;
336}
337
338number ntAdd(number a, number b, const coeffs cf)
339{
340  ntTest(a); ntTest(b);
341  if (IS0(a)) return ntCopy(b, cf);
342  if (IS0(b)) return ntCopy(a, cf);
343 
344  fraction fa = (fraction)a;
345  fraction fb = (fraction)b;
346 
347  poly g = p_Copy(NUM(fa), ntRing);
348  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
349  poly h = p_Copy(NUM(fb), ntRing);
350  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
351  g = p_Add_q(g, h, ntRing);
352 
353  if (g == NULL) return NULL;
354 
355  poly f;
356  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
357  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
358  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
359  else /* both denom's are != 1 */    f = p_Mult_q(p_Copy(DEN(fa), ntRing),
360                                                   p_Copy(DEN(fb), ntRing),
361                                                   ntRing);
362 
363  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
364  NUM(result) = g;
365  DEN(result) = f;
366  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
367  heuristicGcdCancellation((number)result, cf);
368  return (number)result;
369}
370
371number ntSub(number a, number b, const coeffs cf)
372{
373  ntTest(a); ntTest(b);
374  if (IS0(a)) return ntNeg(ntCopy(b, cf), cf);
375  if (IS0(b)) return ntCopy(a, cf);
376 
377  fraction fa = (fraction)a;
378  fraction fb = (fraction)b;
379 
380  poly g = p_Copy(NUM(fa), ntRing);
381  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
382  poly h = p_Copy(NUM(fb), ntRing);
383  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
384  g = p_Add_q(g, p_Neg(h, ntRing), ntRing);
385 
386  if (g == NULL) return NULL;
387 
388  poly f;
389  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
390  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
391  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
392  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
393                                                   p_Copy(DEN(fb), ntRing),
394                                                   ntRing);
395 
396  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
397  NUM(result) = g;
398  DEN(result) = f;
399  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
400  heuristicGcdCancellation((number)result, cf);
401  return (number)result;
402}
403
404number ntMult(number a, number b, const coeffs cf)
405{
406  ntTest(a); ntTest(b);
407  if (IS0(a) || IS0(b)) return NULL;
408 
409  fraction fa = (fraction)a;
410  fraction fb = (fraction)b;
411 
412  poly g = p_Copy(NUM(fa), ntRing);
413  poly h = p_Copy(NUM(fb), ntRing);
414  g = p_Mult_q(g, h, ntRing);
415 
416  if (g == NULL) return NULL;   /* may happen due to zero divisors */
417 
418  poly f;
419  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
420  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
421  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
422  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
423                                                   p_Copy(DEN(fb), ntRing),
424                                                   ntRing);
425 
426  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
427  NUM(result) = g;
428  DEN(result) = f;
429  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
430  heuristicGcdCancellation((number)result, cf);
431  return (number)result;
432}
433
434number ntDiv(number a, number b, const coeffs cf)
435{
436  ntTest(a); ntTest(b);
437  if (IS0(a)) return NULL;
438  if (IS0(b)) WerrorS(nDivBy0);
439 
440  fraction fa = (fraction)a;
441  fraction fb = (fraction)b;
442 
443  poly g = p_Copy(NUM(fa), ntRing);
444  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
445 
446  if (g == NULL) return NULL;   /* may happen due to zero divisors */
447 
448  poly f = p_Copy(NUM(fb), ntRing);
449  if (!DENIS1(fa)) f = p_Mult_q(f, p_Copy(DEN(fa), ntRing), ntRing);
450 
451  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
452  NUM(result) = g;
453  DEN(result) = f;
454  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
455  heuristicGcdCancellation((number)result, cf);
456  return (number)result;
457}
458
459/* 0^0 = 0;
460   for |exp| <= 7 compute power by a simple multiplication loop;
461   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
462      p^13 = p^1 * p^4 * p^8, where we utilise that
463      p^(2^(k+1)) = p^(2^k) * p^(2^k);
464   intermediate cancellation is controlled by the in-place method
465   heuristicGcdCancellation; see there.
466*/
467void ntPower(number a, int exp, number *b, const coeffs cf)
468{
469  ntTest(a);
470 
471  /* special cases first */
472  if (IS0(a))
473  {
474    if (exp >= 0) *b = NULL;
475    else          WerrorS(nDivBy0);
476  }
477  else if (exp ==  0) *b = ntInit(1, cf);
478  else if (exp ==  1) *b = ntCopy(a, cf);
479  else if (exp == -1) *b = ntInvers(a, cf);
480 
481  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
482 
483  /* now compute a^expAbs */
484  number pow; number t;
485  if (expAbs <= 7)
486  {
487    pow = ntCopy(a, cf);
488    for (int i = 2; i <= expAbs; i++)
489    {
490      t = ntMult(pow, a, cf);
491      ntDelete(&pow, cf);
492      pow = t;
493      heuristicGcdCancellation(pow, cf);
494    }
495  }
496  else
497  {
498    pow = ntInit(1, cf);
499    number factor = ntCopy(a, cf);
500    while (expAbs != 0)
501    {
502      if (expAbs & 1)
503      {
504        t = ntMult(pow, factor, cf);
505        ntDelete(&pow, cf);
506        pow = t;
507        heuristicGcdCancellation(pow, cf);
508      }
509      expAbs = expAbs / 2;
510      if (expAbs != 0)
511      {
512        t = ntMult(factor, factor, cf);
513        ntDelete(&factor, cf);
514        factor = t;
515        heuristicGcdCancellation(factor, cf);
516      }
517    }
518    ntDelete(&factor, cf);
519  }
520 
521  /* invert if original exponent was negative */
522  if (exp < 0)
523  {
524    t = ntInvers(pow, cf);
525    ntDelete(&pow, cf);
526    pow = t;
527  }
528  *b = pow;
529}
530
531/* modifies a */
532void heuristicGcdCancellation(number a, const coeffs cf)
533{
534  ntTest(a);
535  if (IS0(a)) return;
536 
537  fraction f = (fraction)a;
538  if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
539 
540  /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
541  if (p_EqualPolys(NUM(f), DEN(f), ntRing))
542  { /* numerator and denominator are both != 1 */
543    p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
544    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
545    COM(f) = 0;
546    return;
547  }
548 
549  if (COM(f) <= BOUND_COMPLEXITY) return;
550  else definiteGcdCancellation(a, cf, TRUE);
551}
552
553/* modifies a */
554void definiteGcdCancellation(number a, const coeffs cf,
555                             BOOLEAN skipSimpleTests)
556{
557  ntTest(a);
558 
559  fraction f = (fraction)a;
560 
561  if (!skipSimpleTests)
562  {
563    if (IS0(a)) return;
564    if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
565 
566    /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
567    if (p_EqualPolys(NUM(f), DEN(f), ntRing))
568    { /* numerator and denominator are both != 1 */
569      p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
570      p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
571      COM(f) = 0;
572      return;
573    }
574  }
575 
576#ifdef HAVE_FACTORY
577  /* singclap_gcd destroys its arguments. We hence need copies: */
578  poly pNum = p_Copy(NUM(f), ntRing);
579  poly pDen = p_Copy(DEN(f), ntRing);
580  poly pGcd = singclap_gcd(pNum, pDen, cf->extRing);
581  if (p_IsConstant(pGcd, ntRing) &&
582      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
583  {
584      /* gcd = 1; nothing to cancel */
585  }
586  else
587  {
588      poly newNum = singclap_pdivide(NUM(f), pGcd, ntRing);
589      p_Delete(&NUM(f), ntRing);
590      NUM(f) = newNum;
591      poly newDen = singclap_pdivide(DEN(f), pGcd, ntRing);
592      p_Delete(&DEN(f), ntRing);
593      if (p_IsConstant(newDen, ntRing) &&
594          n_IsOne(p_GetCoeff(newDen, ntRing), ntCoeffs))
595      {
596        /* newDen = 1 needs to be represented by NULL */
597        p_Delete(&newDen, ntRing);
598        newDen = NULL;
599      }
600      DEN(f) = newDen;
601  }
602  COM(f) = 0;
603  p_Delete(&pGcd, ntRing);
604#endif /* HAVE_FACTORY */
605}
606
607/* modifies a */
608void ntWrite(number &a, const coeffs cf)
609{
610  ntTest(a);
611  definiteGcdCancellation(a, cf, FALSE);
612  if (IS0(a))
613    StringAppendS("0");
614  else
615  {
616    fraction f = (fraction)a;
617    BOOLEAN useBrackets = (!p_IsConstant(NUM(f), ntRing)) ||
618                          (!n_GreaterZero(p_GetCoeff(NUM(f), ntRing),
619                                          ntCoeffs));
620    if (useBrackets) StringAppendS("(");
621    p_String0(NUM(f), ntRing, ntRing);
622    if (useBrackets) StringAppendS(")");
623    if (!DENIS1(f))
624    {
625      StringAppendS("/");
626      useBrackets = (!p_IsConstant(DEN(f), ntRing)) ||
627                    (!n_GreaterZero(p_GetCoeff(DEN(f), ntRing), ntCoeffs));
628      if (useBrackets) StringAppendS("(");
629      p_String0(DEN(f), ntRing, ntRing);
630      if (useBrackets) StringAppendS(")");
631    }
632  }
633}
634
635const char * ntRead(const char *s, number *a, const coeffs cf)
636{
637  poly p;
638  const char * result = p_Read(s, p, ntRing);
639  if (p == NULL) { *a = NULL; return result; }
640  else
641  {
642    fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
643    NUM(f) = p;
644    DEN(f) = NULL;
645    COM(f) = 0;
646    *a = (number)f;
647    return result;
648  }
649}
650
651/* expects *param to be castable to TransExtInfo */
652static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
653{
654  if (ntID != n) return FALSE;
655  TransExtInfo *e = (TransExtInfo *)param;
656  /* for rational function fields we expect the underlying
657     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
658     this expectation is based on the assumption that we have properly
659     registered cf and perform reference counting rather than creating
660     multiple copies of the same coefficient field/domain/ring */
661  return (ntRing == e->r);
662}
663
664number ntLcm(number a, number b, const coeffs cf)
665{
666  ntTest(a); ntTest(b);
667  /* TO BE IMPLEMENTED!
668     for the time, we simply return NULL, representing the number zero */
669  Print("// TO BE IMPLEMENTED: transext.cc:ntLcm\n");
670  return NULL;
671}
672
673number ntGcd(number a, number b, const coeffs cf)
674{
675  ntTest(a); ntTest(b);
676  /* TO BE IMPLEMENTED!
677     for the time, we simply return NULL, representing the number zero */
678  Print("// TO BE IMPLEMENTED: transext.cc:ntGcd\n");
679  return NULL;
680}
681
682int ntSize(number a, const coeffs cf)
683{
684  ntTest(a);
685  if (IS0(a)) return -1;
686  /* this has been taken from the old implementation of field extensions,
687     where we computed the sum of the degrees and the numbers of terms in
688     the numerator and denominator of a; so we leave it at that, for the
689     time being */
690  fraction f = (fraction)a;
691  poly p = NUM(f);
692  int noOfTerms = 0;
693  int numDegree = 0;
694  while (p != NULL)
695  {
696    noOfTerms++;
697    int d = 0;
698    for (int i = 1; i <= rVar(ntRing); i++)
699      d += p_GetExp(p, i, ntRing);
700    if (d > numDegree) numDegree = d;
701    pIter(p);
702  }
703  int denDegree = 0;
704  if (!DENIS1(f))
705  {
706    p = DEN(f);
707    while (p != NULL)
708    {
709      noOfTerms++;
710      int d = 0;
711      for (int i = 1; i <= rVar(ntRing); i++)
712        d += p_GetExp(p, i, ntRing);
713      if (d > denDegree) denDegree = d;
714      pIter(p);
715    }
716  }
717  return numDegree + denDegree + noOfTerms;
718}
719
720number ntInvers(number a, const coeffs cf)
721{
722  ntTest(a);
723  if (IS0(a)) WerrorS(nDivBy0);
724  fraction f = (fraction)a;
725  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
726  poly g;
727  if (DENIS1(f)) g = p_One(ntRing);
728  else           g = p_Copy(DEN(f), ntRing);
729  NUM(result) = g;
730  DEN(result) = p_Copy(NUM(f), ntRing);
731  COM(result) = COM(f);
732  return (number)result;
733}
734
735/* assumes that src = Q, dst = Q(t_1, ..., t_s) */
736number ntMap00(number a, const coeffs src, const coeffs dst)
737{
738  if (n_IsZero(a, src)) return NULL;
739  assume(src == dst->extRing->cf);
740  poly p = p_One(dst->extRing);
741  p_SetCoeff(p, ntCopy(a, src), dst->extRing);
742  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
743  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
744  return (number)f;
745}
746
747/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
748number ntMapP0(number a, const coeffs src, const coeffs dst)
749{
750  if (n_IsZero(a, src)) return NULL;
751  /* mapping via intermediate int: */
752  int n = n_Int(a, src);
753  number q = n_Init(n, dst->extRing->cf);
754  poly p;
755  if (n_IsZero(q, dst->extRing->cf))
756  {
757    n_Delete(&q, dst->extRing->cf);
758    return NULL;
759  }
760  p = p_One(dst->extRing);
761  p_SetCoeff(p, q, dst->extRing);
762  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
763  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
764  return (number)f;
765}
766
767/* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
768                       src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
769number ntCopyMap(number a, const coeffs src, const coeffs dst)
770{
771  return ntCopy(a, dst);
772}
773
774/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
775number ntMap0P(number a, const coeffs src, const coeffs dst)
776{
777  if (n_IsZero(a, src)) return NULL;
778  int p = rChar(dst->extRing);
779  int n = nlModP(a, p, src);
780  number q = n_Init(n, dst->extRing->cf);
781  poly g;
782  if (n_IsZero(q, dst->extRing->cf))
783  {
784    n_Delete(&q, dst->extRing->cf);
785    return NULL;
786  }
787  g = p_One(dst->extRing);
788  p_SetCoeff(g, q, dst->extRing);
789  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
790  NUM(f) = g; DEN(f) = NULL; COM(f) = 0;
791  return (number)f;
792}
793
794/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
795number ntMapPP(number a, const coeffs src, const coeffs dst)
796{
797  if (n_IsZero(a, src)) return NULL;
798  assume(src == dst->extRing->cf);
799  poly p = p_One(dst->extRing);
800  p_SetCoeff(p, ntCopy(a, src), dst->extRing);
801  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
802  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
803  return (number)f;
804}
805
806/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
807number ntMapUP(number a, const coeffs src, const coeffs dst)
808{
809  if (n_IsZero(a, src)) return NULL;
810  /* mapping via intermediate int: */
811  int n = n_Int(a, src);
812  number q = n_Init(n, dst->extRing->cf);
813  poly p;
814  if (n_IsZero(q, dst->extRing->cf))
815  {
816    n_Delete(&q, dst->extRing->cf);
817    return NULL;
818  }
819  p = p_One(dst->extRing);
820  p_SetCoeff(p, q, dst->extRing);
821  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
822  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
823  return (number)f;
824}
825
826nMapFunc ntSetMap(const coeffs src, const coeffs dst)
827{
828  /* dst is expected to be a rational function field */
829  assume(getCoeffType(dst) == ntID);
830 
831  int h = 0; /* the height of the extension tower given by dst */
832  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
833 
834  /* for the time being, we only provide maps if h = 1 and if b is Q or
835     some field Z/pZ: */
836  if (h != 1) return NULL;
837  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
838 
839  /* Let T denote the sequence of transcendental extension variables, i.e.,
840     K[t_1, ..., t_s] =: K[T];
841     Let moreover, for any such sequence T, T' denote any subsequence of T
842     of the form t_1, ..., t_w with w <= s. */
843 
844  if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
845    return ntMap00;                                      /// Q       -->  Q(T)
846 
847  if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
848    return ntMapP0;                                      /// Z/p     -->  Q(T)
849 
850  if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
851    return ntMap0P;                                      /// Q       --> Z/p(T)
852 
853  if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
854  {
855    if (src->ch == dst->ch) return ntMapPP;              /// Z/p     --> Z/p(T)
856    else return ntMapUP;                                 /// Z/u     --> Z/p(T)
857  }
858 
859  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
860  if (h != 1) return NULL;
861  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
862 
863  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
864  {
865    if (rVar(src->extRing) > rVar(dst->extRing)) return NULL;
866    for (int i = 0; i < rVar(src->extRing); i++)
867      if (strcmp(rRingVar(i, src->extRing),
868                 rRingVar(i, dst->extRing)) != 0) return NULL;
869      return ntCopyMap;                                  /// Q(T')   --> Q(T)
870  }
871 
872  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
873  {
874    if (rVar(src->extRing) > rVar(dst->extRing)) return NULL;
875    for (int i = 0; i < rVar(src->extRing); i++)
876      if (strcmp(rRingVar(i, src->extRing),
877                 rRingVar(i, dst->extRing)) != 0) return NULL;
878      return ntCopyMap;                                  /// Z/p(T') --> Z/p(T)
879  }
880 
881  return NULL;                                           /// default
882}
883
884BOOLEAN ntInitChar(coeffs cf, void * infoStruct)
885{ 
886  TransExtInfo *e = (TransExtInfo *)infoStruct;
887  /// first check whether cf->extRing != NULL and delete old ring???
888  cf->extRing           = e->r;
889  cf->extRing->minideal = NULL;
890
891  assume(cf->extRing                != NULL);      // extRing;
892  assume(cf->extRing->cf            != NULL);      // extRing->cf;
893  assume(getCoeffType(cf) == ntID);                // coeff type;
894 
895  /* propagate characteristic up so that it becomes
896     directly accessible in cf: */
897  cf->ch = cf->extRing->cf->ch;
898 
899  cf->cfGreaterZero  = ntGreaterZero;
900  cf->cfGreater      = ntGreater;
901  cf->cfEqual        = ntEqual;
902  cf->cfIsZero       = ntIsZero;
903  cf->cfIsOne        = ntIsOne;
904  cf->cfIsMOne       = ntIsMOne;
905  cf->cfInit         = ntInit;
906  cf->cfInt          = ntInt;
907  cf->cfNeg          = ntNeg;
908  cf->cfPar          = ntPar;
909  cf->cfAdd          = ntAdd;
910  cf->cfSub          = ntSub;
911  cf->cfMult         = ntMult;
912  cf->cfDiv          = ntDiv;
913  cf->cfExactDiv     = ntDiv;
914  cf->cfPower        = ntPower;
915  cf->cfCopy         = ntCopy;
916  cf->cfWrite        = ntWrite;
917  cf->cfRead         = ntRead;
918  cf->cfDelete       = ntDelete;
919  cf->cfSetMap       = ntSetMap;
920  cf->cfGetDenom     = ntGetDenom;
921  cf->cfGetNumerator = ntGetNumerator;
922  cf->cfRePart       = ntCopy;
923  cf->cfImPart       = ntImPart;
924  cf->cfCoeffWrite   = ntCoeffWrite;
925  cf->cfDBTest       = ntDBTest;
926  cf->cfGcd          = ntGcd;
927  cf->cfLcm          = ntLcm;
928  cf->cfSize         = ntSize;
929  cf->nCoeffIsEqual  = ntCoeffIsEqual;
930  cf->cfInvers       = ntInvers;
931  cf->cfIntDiv       = ntDiv;
932 
933#ifndef HAVE_FACTORY
934  PrintS("// Warning: The 'factory' module is not available.\n");
935  PrintS("//          Hence gcd's cannot be cancelled in any\n");
936  PrintS("//          computed fraction!\n");
937#endif
938 
939  return FALSE;
940}
Note: See TracBrowser for help on using the repository browser.