source: git/libpolys/polys/ext_fields/transext.cc @ 1090a98

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