source: git/libpolys/polys/ext_fields/transext.cc @ df43d9

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