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

spielwiese
Last change on this file since d208a80 was 95c826, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: some simple bugs in comparison
  • Property mode set to 100644
File size: 55.6 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*           TODO: the description above needs a major update!!!
34*/
35#define TRANSEXT_PRIVATES
36
37#include "config.h"
38#include <misc/auxiliary.h>
39
40#include <omalloc/omalloc.h>
41
42#include <reporter/reporter.h>
43
44#include <coeffs/coeffs.h>
45#include <coeffs/numbers.h>
46#include <coeffs/longrat.h>
47
48#include <polys/monomials/ring.h>
49#include <polys/monomials/p_polys.h>
50#include <polys/simpleideals.h>
51
52#ifdef HAVE_FACTORY
53#include <polys/clapsing.h>
54#include <polys/clapconv.h>
55#include <factory/factory.h>
56#endif
57
58#include <polys/ext_fields/transext.h>
59#include <polys/prCopy.h>
60
61#include <polys/PolyEnumerator.h>
62
63
64/* constants for controlling the complexity of numbers */
65#define ADD_COMPLEXITY 1   /**< complexity increase due to + and - */
66#define MULT_COMPLEXITY 2   /**< complexity increase due to * and / */
67#define BOUND_COMPLEXITY 10   /**< maximum complexity of a number */
68
69
70static inline BOOLEAN p_IsOne(const poly p, const ring R)
71{
72  assume( p_Test(p, R) );
73  return (p_IsConstant(p, R) && n_IsOne(p_GetCoeff(p, R), R->cf));
74}
75
76/// TRUE iff num. represents 1
77#define NUMIS1(f) (p_IsOne(NUM(f), cf->extRing))
78
79#define COM(f) f->complexity
80
81
82#ifdef LDEBUG
83#define ntTest(a) assume(ntDBTest(a,__FILE__,__LINE__,cf))
84BOOLEAN  ntDBTest(number a, const char *f, const int l, const coeffs r);
85#else
86#define ntTest(a) do {} while (0)
87#endif
88
89/// Our own type!
90static const n_coeffType ID = n_transExt;
91
92/* polynomial ring in which the numerators and denominators of our
93   numbers live */
94#define ntRing cf->extRing
95
96/* coeffs object in which the coefficients of our numbers live;
97 * methods attached to ntCoeffs may be used to compute with the
98 * coefficients of our numbers, e.g., use ntCoeffs->nAdd to add
99 * coefficients of our numbers */
100#define ntCoeffs cf->extRing->cf
101
102
103omBin fractionObjectBin = omGetSpecBin(sizeof(fractionObject));
104
105/// forward declarations
106BOOLEAN  ntGreaterZero(number a, const coeffs cf);
107BOOLEAN  ntGreater(number a, number b, const coeffs cf);
108BOOLEAN  ntEqual(number a, number b, const coeffs cf);
109BOOLEAN  ntIsOne(number a, const coeffs cf);
110BOOLEAN  ntIsMOne(number a, const coeffs cf);
111BOOLEAN  ntIsZero(number a, const coeffs cf);
112number   ntInit(long i, const coeffs cf);
113int      ntInt(number &a, const coeffs cf);
114number   ntNeg(number a, const coeffs cf);
115number   ntInvers(number a, const coeffs cf);
116number   ntAdd(number a, number b, const coeffs cf);
117number   ntSub(number a, number b, const coeffs cf);
118number   ntMult(number a, number b, const coeffs cf);
119number   ntDiv(number a, number b, const coeffs cf);
120void     ntPower(number a, int exp, number *b, const coeffs cf);
121number   ntCopy(number a, const coeffs cf);
122void     ntWriteLong(number &a, const coeffs cf);
123void     ntWriteShort(number &a, const coeffs cf);
124number   ntRePart(number a, const coeffs cf);
125number   ntImPart(number a, const coeffs cf);
126number   ntGetDenom(number &a, const coeffs cf);
127number   ntGetNumerator(number &a, const coeffs cf);
128number   ntGcd(number a, number b, const coeffs cf);
129number   ntLcm(number a, number b, const coeffs cf);
130int      ntSize(number a, const coeffs cf);
131void     ntDelete(number * a, const coeffs cf);
132void     ntCoeffWrite(const coeffs cf, BOOLEAN details);
133number   ntIntDiv(number a, number b, const coeffs cf);
134const char * ntRead(const char *s, number *a, const coeffs cf);
135static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param);
136
137void heuristicGcdCancellation(number a, const coeffs cf);
138void definiteGcdCancellation(number a, const coeffs cf,
139                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed);
140void handleNestedFractionsOverQ(fraction f, const coeffs cf);
141
142#ifdef LDEBUG
143BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs cf)
144{
145  assume(getCoeffType(cf) == ID);
146
147  if (IS0(a)) return TRUE;
148
149  const fraction t = (fraction)a;
150
151  const poly num = NUM(t);
152  assume(num != NULL);   /**< t != 0 ==> numerator(t) != 0 */
153  assume( _p_Test(num, ntRing,1) );
154
155  const poly den = DEN(t);
156
157  if (den != NULL) // !DENIS1(f)
158  {
159    assume( _p_Test(den, ntRing,1) );
160
161    if(p_IsConstant(den, ntRing) && (n_IsOne(pGetCoeff(den), ntRing->cf)))
162    {
163      Print("?/1 in %s:%d\n",f,l);
164      return FALSE;
165    }
166
167    if( !n_GreaterZero(pGetCoeff(den), ntRing->cf) )
168    {
169      Print("negative sign of DEN. of a fraction in %s:%d\n",f,l);
170      return FALSE;
171    }
172
173    // test that den is over integers!?
174
175  } else
176  {  // num != NULL // den == NULL
177
178//    if( COM(t) != 0 )
179//    {
180//      Print("?//NULL with non-zero complexity: %d in %s:%d\n", COM(t), f, l);
181//      return FALSE;
182//    }
183    // test that nume is over integers!?
184  }
185  return TRUE;
186}
187#endif
188
189/* returns the bottom field in this field extension tower; if the tower
190   is flat, i.e., if there is no extension, then r itself is returned;
191   as a side-effect, the counter 'height' is filled with the height of
192   the extension tower (in case the tower is flat, 'height' is zero) */
193static coeffs nCoeff_bottom(const coeffs r, int &height)
194{
195  assume(r != NULL);
196  coeffs cf = r;
197  height = 0;
198  while (nCoeff_is_Extension(cf))
199  {
200    assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL);
201    cf = cf->extRing->cf;
202    height++;
203  }
204  return cf;
205}
206
207BOOLEAN ntIsZero(number a, const coeffs cf)
208{
209  ntTest(a); // !!!
210  return (IS0(a));
211}
212
213void ntDelete(number * a, const coeffs cf)
214{
215  ntTest(*a); // !!!
216  fraction f = (fraction)(*a);
217  if (IS0(f)) return;
218  p_Delete(&NUM(f), ntRing);
219  if (!DENIS1(f)) p_Delete(&DEN(f), ntRing);
220  omFreeBin((ADDRESS)f, fractionObjectBin);
221  *a = NULL;
222}
223
224BOOLEAN ntEqual(number a, number b, const coeffs cf)
225{
226  ntTest(a);
227  ntTest(b);
228
229  /// simple tests
230  if (a == b) return TRUE;
231  if ((IS0(a)) && (!IS0(b))) return FALSE;
232  if ((IS0(b)) && (!IS0(a))) return FALSE;
233
234  /// cheap test if gcd's have been cancelled in both numbers
235  fraction fa = (fraction)a;
236  fraction fb = (fraction)b;
237  if ((COM(fa) == 1) && (COM(fb) == 1))
238  {
239    poly f = p_Add_q(p_Copy(NUM(fa), ntRing),
240                     p_Neg(p_Copy(NUM(fb), ntRing), ntRing),
241                     ntRing);
242    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
243    if (DENIS1(fa) && DENIS1(fb))  return TRUE;
244    if (DENIS1(fa) && !DENIS1(fb)) return FALSE;
245    if (!DENIS1(fa) && DENIS1(fb)) return FALSE;
246    f = p_Add_q(p_Copy(DEN(fa), ntRing),
247                p_Neg(p_Copy(DEN(fb), ntRing), ntRing),
248                ntRing);
249    if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
250    return TRUE;
251  }
252
253  /* default: the more expensive multiplication test
254              a/b = c/d  <==>  a*d = b*c */
255  poly f = p_Copy(NUM(fa), ntRing);
256  if (!DENIS1(fb)) f = p_Mult_q(f, p_Copy(DEN(fb), ntRing), ntRing);
257  poly g = p_Copy(NUM(fb), ntRing);
258  if (!DENIS1(fa)) g = p_Mult_q(g, p_Copy(DEN(fa), ntRing), ntRing);
259  poly h = p_Add_q(f, p_Neg(g, ntRing), ntRing);
260  if (h == NULL) return TRUE;
261  else
262  {
263    p_Delete(&h, ntRing);
264    return FALSE;
265  }
266}
267
268number ntCopy(number a, const coeffs cf)
269{
270  ntTest(a); // !!!
271  if (IS0(a)) return NULL;
272  fraction f = (fraction)a;
273  poly g = p_Copy(NUM(f), ntRing);
274  poly h = NULL; if (!DENIS1(f)) h = p_Copy(DEN(f), ntRing);
275  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
276  NUM(result) = g;
277  DEN(result) = h;
278  COM(result) = COM(f);
279  ntTest((number)result);
280  return (number)result;
281}
282
283/// TODO: normalization of a!?
284number ntGetNumerator(number &a, const coeffs cf)
285{
286  ntTest(a);
287  definiteGcdCancellation(a, cf, FALSE);
288
289  if (IS0(a)) return NULL;
290
291  fraction f = (fraction)a;
292  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
293
294  const BOOLEAN denis1= DENIS1 (f);
295
296  if (getCoeffType (ntCoeffs) == n_Q && !denis1)
297    handleNestedFractionsOverQ (f, cf);
298
299  if (getCoeffType (ntCoeffs) == n_Q && denis1)
300  {
301    assume( DEN (f) == NULL );
302
303    number g;
304    // TODO/NOTE: the following should not be necessary (due to
305    // Hannes!) as NUM (f) should be over Z!!!
306    CPolyCoeffsEnumerator itr(NUM(f));
307
308
309    n_ClearDenominators(itr, g, ntRing->cf);
310
311    if( !n_GreaterZero(g, ntRing->cf) )
312    {
313      NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
314      g = n_Neg(g, ntRing->cf);
315    }
316
317    // g should be a positive integer now!
318    assume( n_GreaterZero(g, ntRing->cf) );
319
320    if( !n_IsOne(g, ntRing->cf) )
321    {
322      DEN (f) = p_NSet(g, ntRing); // update COM(f)???
323      COM (f) ++;
324      assume( DEN (f) != NULL );
325    }
326    else
327      n_Delete(&g, ntRing->cf);
328
329    ntTest(a);
330  }
331
332  // Call ntNormalize instead of above?!?
333
334  NUM (result) = p_Copy (NUM (f), ntRing); // ???
335  DEN (result) = NULL;
336  COM (result) = 0;
337
338  ntTest((number)result);
339  return (number)result;
340}
341
342/// TODO: normalization of a!?
343number ntGetDenom(number &a, const coeffs cf)
344{
345  ntTest(a);
346  definiteGcdCancellation(a, cf, FALSE);
347  fraction f = (fraction)a;
348
349  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
350  DEN (result)= NULL;
351  COM (result)= 0;
352
353
354  const BOOLEAN denis1 = DENIS1 (f);
355
356  if( IS0(f) || (denis1 && getCoeffType (ntCoeffs) != n_Q) ) // */1 or 0
357  {
358    NUM (result)= p_One(ntRing);
359    ntTest((number)result);
360    return (number)result;
361  }
362
363  if (!denis1) // */* / Q
364  {
365    assume( DEN (f) != NULL );
366
367    if (getCoeffType (ntCoeffs) == n_Q)
368      handleNestedFractionsOverQ (f, cf);
369
370    ntTest(a);
371
372    if( DEN (f) != NULL ) // is it ?? // 1 now???
373    {
374      assume( !p_IsOne(DEN (f), ntRing) );
375
376      NUM (result) = p_Copy (DEN (f), ntRing);
377      ntTest((number)result);
378      return (number)result;
379    }
380//    NUM (result) = p_One(ntRing); // NOTE: just in order to be sure...
381  }
382
383  // */1 / Q
384  assume( getCoeffType (ntCoeffs) == n_Q );
385  assume( DEN (f) == NULL );
386
387  number g;
388//    poly num= p_Copy (NUM (f), ntRing); // ???
389
390
391  // TODO/NOTE: the following should not be necessary (due to
392  // Hannes!) as NUM (f) should be over Z!!!
393  CPolyCoeffsEnumerator itr(NUM(f));
394
395  n_ClearDenominators(itr, g, ntRing->cf); // may return -1 :(((
396
397  if( !n_GreaterZero(g, ntRing->cf) )
398  {
399//     NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
400//     g = n_Neg(g, ntRing->cf);
401    NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
402    g = n_Neg(g, ntRing->cf);
403  }
404
405  // g should be a positive integer now!
406  assume( n_GreaterZero(g, ntRing->cf) );
407
408  if( !n_IsOne(g, ntRing->cf) )
409  {
410    assume( n_GreaterZero(g, ntRing->cf) );
411    assume( !n_IsOne(g, ntRing->cf) );
412
413    DEN (f) = p_NSet(g, ntRing); // update COM(f)???
414    assume( DEN (f) != NULL );
415    COM (f) ++;
416
417    NUM (result)= p_Copy (DEN (f), ntRing);
418  }
419  else
420  { // common denom == 1?
421    NUM (result)= p_NSet(g, ntRing); // p_Copy (DEN (f), ntRing);
422//  n_Delete(&g, ntRing->cf);
423  }
424
425//    if (!p_IsConstant (num, ntRing) && pNext(num) != NULL)
426//    else
427//      g= p_GetAllDenom (num, ntRing);
428//    result= (fraction) ntSetMap (ntRing->cf, cf) (g, ntRing->cf, cf);
429
430  ntTest((number)result);
431  return (number)result;
432}
433
434BOOLEAN ntIsOne(number a, const coeffs cf)
435{
436  ntTest(a); // !!!
437  definiteGcdCancellation(a, cf, FALSE);
438  fraction f = (fraction)a;
439  return (f!=NULL) && DENIS1(f) && NUMIS1(f);
440}
441
442BOOLEAN ntIsMOne(number a, const coeffs cf)
443{
444  ntTest(a);
445  definiteGcdCancellation(a, cf, FALSE);
446  fraction f = (fraction)a;
447  if ((f==NULL) || (!DENIS1(f))) return FALSE;
448  poly g = NUM(f);
449  if (!p_IsConstant(g, ntRing)) return FALSE;
450  return n_IsMOne(p_GetCoeff(g, ntRing), ntCoeffs);
451}
452
453/// this is in-place, modifies a
454number ntNeg(number a, const coeffs cf)
455{
456  ntTest(a);
457  if (!IS0(a))
458  {
459    fraction f = (fraction)a;
460    NUM(f) = p_Neg(NUM(f), ntRing);
461  }
462  ntTest(a);
463  return a;
464}
465
466number ntImPart(number a, const coeffs cf)
467{
468  ntTest(a);
469  return NULL;
470}
471
472number ntInit_bigint(number longratBigIntNumber, const coeffs src, const coeffs cf)
473{
474  assume( cf != NULL );
475
476  const ring A = cf->extRing;
477
478  assume( A != NULL );
479
480  const coeffs C = A->cf;
481
482  assume( C != NULL );
483
484  number n = n_Init_bigint(longratBigIntNumber, src, C);
485
486  if ( n_IsZero(n, C) )
487  {
488    n_Delete(&n, C);
489    return NULL;
490  }
491
492  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
493
494  number den = n_GetDenom(n, C);
495
496  assume( n_GreaterZero(den, C) );
497
498  if( n_IsOne(den, C) )
499  {
500     NUM(result) = p_NSet(n, A);
501     DEN(result) = NULL;
502     n_Delete(&den, C);
503  } else
504  {
505     DEN(result) = p_NSet(den, A);
506     NUM(result) = p_NSet(n_GetNumerator(n, C), A);
507     n_Delete(&n, C);
508  }
509
510  COM(result) = 0;
511
512  ntTest((number)result);
513
514  return (number)result;
515}
516
517
518number ntInit(long i, const coeffs cf)
519{
520  if (i != 0)
521  {
522    poly p=p_ISet(i, ntRing);
523    if (p!=NULL)
524    {
525      fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
526      NUM(result) = p;
527      //DEN(result) = NULL; // done by omAlloc0Bin
528      //COM(result) = 0; // done by omAlloc0Bin
529      ntTest((number)result);
530      return (number)result;
531    }
532  }
533  return NULL;
534}
535
536
537/// takes over p!
538number ntInit(poly p, const coeffs cf)
539{
540  if (p == NULL) return NULL;
541
542    number g;
543    // TODO/NOTE: the following should not be necessary (due to
544    // Hannes!) as NUM (f) should be over Z!!!
545    CPolyCoeffsEnumerator itr(p);
546
547    n_ClearDenominators(itr, g, ntRing->cf);
548
549    if( !n_GreaterZero(g, ntRing->cf) )
550    {
551      p = p_Neg(p, ntRing); // Ugly :(((
552      g = n_Neg(g, ntRing->cf);
553    }
554
555    // g should be a positive integer now!
556    assume( n_GreaterZero(g, ntRing->cf) );
557
558    fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
559
560    if( !n_IsOne(g, ntRing->cf) )
561    {
562      DEN (f) = p_NSet(g, ntRing);
563//      COM (f) ++; // update COM(f)???
564      assume( DEN (f) != NULL );
565    }
566    else
567    {
568      DEN(f) = NULL;
569      n_Delete(&g, ntRing->cf);
570    }
571
572    NUM(f) = p;
573    COM(f) = 0;
574
575    ntTest((number)f);
576    return (number)f;
577}
578
579int ntInt(number &a, const coeffs cf)
580{
581  ntTest(a);
582  if (IS0(a)) return 0;
583  definiteGcdCancellation(a, cf, FALSE);
584  fraction f = (fraction)a;
585  if (!DENIS1(f)) return 0;
586
587  const poly aAsPoly = NUM(f);
588
589  if(aAsPoly == NULL)
590    return 0;
591
592  if (!p_IsConstant(aAsPoly, ntRing))
593    return 0;
594
595  assume( aAsPoly != NULL );
596
597  return n_Int(p_GetCoeff(aAsPoly, ntRing), ntCoeffs);
598}
599
600/* This method will only consider the numerators of a and b, without
601   cancelling gcd's before.
602   Moreover it may return TRUE only if one or both numerators
603   are zero or if their degrees are equal. Then TRUE is returned iff
604   coeff(numerator(a)) > coeff(numerator(b));
605   In all other cases, FALSE will be returned. */
606BOOLEAN ntGreater(number a, number b, const coeffs cf)
607{
608  ntTest(a);
609  ntTest(b);
610  number aNumCoeff = NULL; int aNumDeg = -1;
611  number bNumCoeff = NULL; int bNumDeg = -1;
612  if (!IS0(a))
613  {
614    fraction fa = (fraction)a;
615    aNumDeg = p_Totaldegree(NUM(fa), ntRing);
616    aNumCoeff = p_GetCoeff(NUM(fa), ntRing);
617  }
618  else return !(ntGreaterZero (b,cf));
619  if (!IS0(b))
620  {
621    fraction fb = (fraction)b;
622    bNumDeg = p_Totaldegree(NUM(fb), ntRing);
623    bNumCoeff = p_GetCoeff(NUM(fb), ntRing);
624  }
625  else return ntGreaterZero(a,cf);
626  if (aNumDeg > bNumDeg) return TRUE;
627  if (aNumDeg < bNumDeg) return FALSE;
628  return n_Greater(aNumCoeff, bNumCoeff, ntCoeffs);
629}
630
631/* this method will only consider the numerator of a, without cancelling
632   the gcd before;
633   returns TRUE iff the leading coefficient of the numerator of a is > 0
634                    or the leading term of the numerator of a is not a
635                    constant */
636BOOLEAN ntGreaterZero(number a, const coeffs cf)
637{
638  ntTest(a);
639  if (IS0(a)) return FALSE;
640  fraction f = (fraction)a;
641  poly g = NUM(f);
642  return (n_GreaterZero(p_GetCoeff(g, ntRing), ntCoeffs) ||
643          (!p_LmIsConstant(g, ntRing)));
644}
645
646void ntCoeffWrite(const coeffs cf, BOOLEAN details)
647{
648  assume( cf != NULL );
649
650  const ring A = cf->extRing;
651
652  assume( A != NULL );
653  assume( A->cf != NULL );
654
655  n_CoeffWrite(A->cf, details);
656
657//  rWrite(A);
658
659  const int P = rVar(A);
660  assume( P > 0 );
661
662  Print("//   %d parameter    : ", P);
663
664  for (int nop=0; nop < P; nop ++)
665    Print("%s ", rRingVar(nop, A));
666
667  assume( A->qideal == NULL );
668
669  PrintS("\n//   minpoly        : 0\n");
670
671/*
672  PrintS("//   Coefficients live in the rational function field\n");
673  Print("//   K(");
674  for (int i = 0; i < rVar(ntRing); i++)
675  {
676    if (i > 0) PrintS(" ");
677    Print("%s", rRingVar(i, ntRing));
678  }
679  PrintS(") with\n");
680  PrintS("//   K: "); n_CoeffWrite(cf->extRing->cf);
681*/
682}
683
684number ntAdd(number a, number b, const coeffs cf)
685{
686  ntTest(a);
687  ntTest(b);
688  if (IS0(a)) return ntCopy(b, cf);
689  if (IS0(b)) return ntCopy(a, cf);
690
691  fraction fa = (fraction)a;
692  fraction fb = (fraction)b;
693
694  poly g = p_Copy(NUM(fa), ntRing);
695  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
696  poly h = p_Copy(NUM(fb), ntRing);
697  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
698  g = p_Add_q(g, h, ntRing);
699
700  if (g == NULL) return NULL;
701
702  poly f;
703  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
704  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
705  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
706  else /* both denom's are != 1 */    f = p_Mult_q(p_Copy(DEN(fa), ntRing),
707                                                   p_Copy(DEN(fb), ntRing),
708                                                   ntRing);
709
710  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
711  NUM(result) = g;
712  DEN(result) = f;
713  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
714  heuristicGcdCancellation((number)result, cf);
715
716//  ntTest((number)result);
717
718  return (number)result;
719}
720
721number ntSub(number a, number b, const coeffs cf)
722{
723  ntTest(a);
724  ntTest(b);
725  if (IS0(a)) return ntNeg(ntCopy(b, cf), cf);
726  if (IS0(b)) return ntCopy(a, cf);
727
728  fraction fa = (fraction)a;
729  fraction fb = (fraction)b;
730
731  poly g = p_Copy(NUM(fa), ntRing);
732  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
733  poly h = p_Copy(NUM(fb), ntRing);
734  if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
735  g = p_Add_q(g, p_Neg(h, ntRing), ntRing);
736
737  if (g == NULL) return NULL;
738
739  poly f;
740  if      (DENIS1(fa) && DENIS1(fb))  f = NULL;
741  else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
742  else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
743  else /* both den's are != 1 */      f = p_Mult_q(p_Copy(DEN(fa), ntRing),
744                                                   p_Copy(DEN(fb), ntRing),
745                                                   ntRing);
746
747  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
748  NUM(result) = g;
749  DEN(result) = f;
750  COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
751  heuristicGcdCancellation((number)result, cf);
752//  ntTest((number)result);
753  return (number)result;
754}
755
756number ntMult(number a, number b, const coeffs cf)
757{
758  ntTest(a); // !!!?
759  ntTest(b); // !!!?
760
761  if (IS0(a) || IS0(b)) return NULL;
762
763  fraction fa = (fraction)a;
764  fraction fb = (fraction)b;
765
766  const poly g = pp_Mult_qq(NUM(fa), NUM(fb), ntRing);
767
768  if (g == NULL) return NULL; // may happen due to zero divisors???
769
770  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
771
772  NUM(result) = g;
773
774  const poly da = DEN(fa);
775  const poly db = DEN(fb);
776
777
778  if (db == NULL)
779  {
780    // b = ? // NULL
781
782    if(da == NULL)
783    { // both fa && fb are ?? // NULL!
784      assume (da == NULL && db == NULL);
785      DEN(result) = NULL;
786      COM(result) = 0;
787    }
788    else
789    {
790      assume (da != NULL && db == NULL);
791      DEN(result) = p_Copy(da, ntRing);
792      COM(result) = COM(fa) + MULT_COMPLEXITY;
793      heuristicGcdCancellation((number)result, cf);
794    }
795  } else
796  { // b = ?? / ??
797    if (da == NULL)
798    { // a == ? // NULL
799      assume( db != NULL && da == NULL);
800      DEN(result) = p_Copy(db, ntRing);
801      COM(result) = COM(fb) + MULT_COMPLEXITY;
802      heuristicGcdCancellation((number)result, cf);
803    }
804    else /* both den's are != 1 */
805    {
806      assume (da != NULL && db != NULL);
807      DEN(result) = pp_Mult_qq(da, db, ntRing);
808      COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
809      heuristicGcdCancellation((number)result, cf);
810    }
811  }
812
813//  ntTest((number)result);
814
815  return (number)result;
816}
817
818number ntDiv(number a, number b, const coeffs cf)
819{
820  ntTest(a);
821  ntTest(b);
822  if (IS0(a)) return NULL;
823  if (IS0(b)) WerrorS(nDivBy0);
824
825  fraction fa = (fraction)a;
826  fraction fb = (fraction)b;
827
828  poly g = p_Copy(NUM(fa), ntRing);
829  if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
830
831  if (g == NULL) return NULL;   /* may happen due to zero divisors */
832
833  poly f = p_Copy(NUM(fb), ntRing);
834  if (!DENIS1(fa)) f = p_Mult_q(f, p_Copy(DEN(fa), ntRing), ntRing);
835
836  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
837  NUM(result) = g;
838  if (!p_IsConstant(f,ntRing) || !n_IsOne(pGetCoeff(f),ntRing->cf))
839    DEN(result) = f;
840  COM(result) = COM(fa) + COM(fb) + MULT_COMPLEXITY;
841  heuristicGcdCancellation((number)result, cf);
842//  ntTest((number)result);
843  return (number)result;
844}
845
846/* 0^0 = 0;
847   for |exp| <= 7 compute power by a simple multiplication loop;
848   for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
849      p^13 = p^1 * p^4 * p^8, where we utilise that
850      p^(2^(k+1)) = p^(2^k) * p^(2^k);
851   intermediate cancellation is controlled by the in-place method
852   heuristicGcdCancellation; see there.
853*/
854void ntPower(number a, int exp, number *b, const coeffs cf)
855{
856  ntTest(a);
857
858  /* special cases first */
859  if (IS0(a))
860  {
861    if (exp >= 0) *b = NULL;
862    else          WerrorS(nDivBy0);
863  }
864  else if (exp ==  0) { *b = ntInit(1, cf); return;}
865  else if (exp ==  1) { *b = ntCopy(a, cf); return;}
866  else if (exp == -1) { *b = ntInvers(a, cf); return;}
867
868  int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
869
870  /* now compute a^expAbs */
871  number pow; number t;
872  if (expAbs <= 7)
873  {
874    pow = ntCopy(a, cf);
875    for (int i = 2; i <= expAbs; i++)
876    {
877      t = ntMult(pow, a, cf);
878      ntDelete(&pow, cf);
879      pow = t;
880      heuristicGcdCancellation(pow, cf);
881    }
882  }
883  else
884  {
885    pow = ntInit(1, cf);
886    number factor = ntCopy(a, cf);
887    while (expAbs != 0)
888    {
889      if (expAbs & 1)
890      {
891        t = ntMult(pow, factor, cf);
892        ntDelete(&pow, cf);
893        pow = t;
894        heuristicGcdCancellation(pow, cf);
895      }
896      expAbs = expAbs / 2;
897      if (expAbs != 0)
898      {
899        t = ntMult(factor, factor, cf);
900        ntDelete(&factor, cf);
901        factor = t;
902        heuristicGcdCancellation(factor, cf);
903      }
904    }
905    ntDelete(&factor, cf);
906  }
907
908  /* invert if original exponent was negative */
909  if (exp < 0)
910  {
911    t = ntInvers(pow, cf);
912    ntDelete(&pow, cf);
913    pow = t;
914  }
915  *b = pow;
916  ntTest(*b);
917}
918
919/* assumes that cf represents the rationals, i.e. Q, and will only
920   be called in that case;
921   assumes furthermore that f != NULL and that the denominator of f != 1;
922   generally speaking, this method removes denominators in the rational
923   coefficients of the numerator and denominator of 'a';
924   more concretely, the following normalizations will be performed,
925   where t^alpha denotes a monomial in the transcendental variables t_k
926   (1) if 'a' is of the form
927          (sum_alpha a_alpha/b_alpha * t^alpha)
928          -------------------------------------
929            (sum_beta c_beta/d_beta * t^beta)
930       with integers a_alpha, b_alpha, c_beta, d_beta, then both the
931       numerator and the denominator will be multiplied by the LCM of
932       the b_alpha's and the d_beta's (if this LCM is != 1),
933   (2) if 'a' is - e.g. after having performed step (1) - of the form
934          (sum_alpha a_alpha * t^alpha)
935          -----------------------------
936            (sum_beta c_beta * t^beta)
937       with integers a_alpha, c_beta, and with a non-constant denominator,
938       then both the numerator and the denominator will be divided by the
939       GCD of the a_alpha's and the c_beta's (if this GCD is != 1),
940   (3) if 'a' is - e.g. after having performed steps (1) and (2) - of the
941       form
942          (sum_alpha a_alpha * t^alpha)
943          -----------------------------
944                        c
945       with integers a_alpha, and c != 1, then 'a' will be replaced by
946       (sum_alpha a_alpha/c * t^alpha);
947   this procedure does not alter COM(f) (this has to be done by the
948   calling procedure);
949   modifies f */
950void handleNestedFractionsOverQ(fraction f, const coeffs cf)
951{
952  assume(nCoeff_is_Q(ntCoeffs));
953  assume(!IS0(f));
954  assume(!DENIS1(f));
955
956  if (!p_IsConstant(DEN(f), ntRing))
957  { /* step (1); see documentation of this procedure above */
958    p_Normalize(NUM(f), ntRing);
959    p_Normalize(DEN(f), ntRing);
960    number lcmOfDenominators = n_Init(1, ntCoeffs);
961    number c; number tmp;
962    poly p = NUM(f);
963    /* careful when using n_Lcm!!! It computes the lcm of the numerator
964       of the 1st argument and the denominator of the 2nd!!! */
965    while (p != NULL)
966    {
967      c = p_GetCoeff(p, ntRing);
968      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
969      n_Delete(&lcmOfDenominators, ntCoeffs);
970      lcmOfDenominators = tmp;
971      pIter(p);
972    }
973    p = DEN(f);
974    while (p != NULL)
975    {
976      c = p_GetCoeff(p, ntRing);
977      tmp = n_Lcm(lcmOfDenominators, c, ntCoeffs);
978      n_Delete(&lcmOfDenominators, ntCoeffs);
979      lcmOfDenominators = tmp;
980      pIter(p);
981    }
982    if (!n_IsOne(lcmOfDenominators, ntCoeffs))
983    { /* multiply NUM(f) and DEN(f) with lcmOfDenominators */
984      NUM(f) = p_Mult_nn(NUM(f), lcmOfDenominators, ntRing);
985      p_Normalize(NUM(f), ntRing);
986      DEN(f) = p_Mult_nn(DEN(f), lcmOfDenominators, ntRing);
987      p_Normalize(DEN(f), ntRing);
988    }
989    n_Delete(&lcmOfDenominators, ntCoeffs);
990    if (!p_IsConstant(DEN(f), ntRing))
991    { /* step (2); see documentation of this procedure above */
992      p = NUM(f);
993      number gcdOfCoefficients = n_Copy(p_GetCoeff(p, ntRing), ntCoeffs);
994      pIter(p);
995      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
996      {
997        c = p_GetCoeff(p, ntRing);
998        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
999        n_Delete(&gcdOfCoefficients, ntCoeffs);
1000        gcdOfCoefficients = tmp;
1001        pIter(p);
1002      }
1003      p = DEN(f);
1004      while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
1005      {
1006        c = p_GetCoeff(p, ntRing);
1007        tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
1008        n_Delete(&gcdOfCoefficients, ntCoeffs);
1009        gcdOfCoefficients = tmp;
1010        pIter(p);
1011      }
1012      if (!n_IsOne(gcdOfCoefficients, ntCoeffs))
1013      { /* divide NUM(f) and DEN(f) by gcdOfCoefficients */
1014        number inverseOfGcdOfCoefficients = n_Invers(gcdOfCoefficients,
1015                                                     ntCoeffs);
1016        NUM(f) = p_Mult_nn(NUM(f), inverseOfGcdOfCoefficients, ntRing);
1017        p_Normalize(NUM(f), ntRing);
1018        DEN(f) = p_Mult_nn(DEN(f), inverseOfGcdOfCoefficients, ntRing);
1019        p_Normalize(DEN(f), ntRing);
1020        n_Delete(&inverseOfGcdOfCoefficients, ntCoeffs);
1021      }
1022      n_Delete(&gcdOfCoefficients, ntCoeffs);
1023    }
1024  }
1025  if (p_IsConstant(DEN(f), ntRing) &&
1026      (!n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs)))
1027  { /* step (3); see documentation of this procedure above */
1028    number inverseOfDen = n_Invers(p_GetCoeff(DEN(f), ntRing), ntCoeffs);
1029    NUM(f) = p_Mult_nn(NUM(f), inverseOfDen, ntRing);
1030    n_Delete(&inverseOfDen, ntCoeffs);
1031    p_Delete(&DEN(f), ntRing);
1032    DEN(f) = NULL;
1033  }
1034
1035  /* Now, due to the above computations, DEN(f) may have become the
1036     1-polynomial which needs to be represented by NULL: */
1037  if ((DEN(f) != NULL) &&
1038      p_IsConstant(DEN(f), ntRing) &&
1039      n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
1040  {
1041    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1042  }
1043
1044  if( DEN(f) != NULL )
1045    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1046    {
1047      NUM(f) = p_Neg(NUM(f), ntRing);
1048      DEN(f) = p_Neg(DEN(f), ntRing);
1049    }
1050
1051  ntTest((number)f); // TODO!
1052}
1053
1054/* modifies a */
1055void heuristicGcdCancellation(number a, const coeffs cf)
1056{
1057//  ntTest(a); // !!!!????
1058  if (IS0(a)) return;
1059
1060  fraction f = (fraction)a;
1061  if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; ntTest(a); return; }
1062
1063  assume( DEN(f) != NULL );
1064
1065  /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
1066  if (p_EqualPolys(NUM(f), DEN(f), ntRing))
1067  { /* numerator and denominator are both != 1 */
1068    p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
1069    p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1070    COM(f) = 0;
1071  } else
1072  {
1073    if (COM(f) > BOUND_COMPLEXITY)
1074      definiteGcdCancellation(a, cf, TRUE);
1075
1076  // TODO: check if it is enough to put the following into definiteGcdCancellation?!
1077  if( DEN(f) != NULL )
1078    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1079    {
1080      NUM(f) = p_Neg(NUM(f), ntRing);
1081      DEN(f) = p_Neg(DEN(f), ntRing);
1082    }
1083  }
1084
1085
1086  ntTest(a); // !!!!????
1087}
1088
1089/// modifies a
1090void definiteGcdCancellation(number a, const coeffs cf,
1091                             BOOLEAN simpleTestsHaveAlreadyBeenPerformed)
1092{
1093  ntTest(a); // !!!!
1094
1095  fraction f = (fraction)a;
1096
1097  if (!simpleTestsHaveAlreadyBeenPerformed)
1098  {
1099    if (IS0(a)) return;
1100    if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
1101
1102    /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
1103    if (p_EqualPolys(NUM(f), DEN(f), ntRing))
1104    { /* numerator and denominator are both != 1 */
1105      p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
1106      p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1107      COM(f) = 0;
1108      ntTest(a); // !!!!
1109      return;
1110    }
1111  }
1112
1113#ifdef HAVE_FACTORY
1114  /* Note that, over Q, singclap_gcd will remove the denominators in all
1115     rational coefficients of pNum and pDen, before starting to compute
1116     the gcd. Thus, we do not need to ensure that the coefficients of
1117     pNum and pDen live in Z; they may well be elements of Q\Z. */
1118  /* singclap_gcd destroys its arguments; we hence need copies: */
1119  poly pGcd = singclap_gcd(p_Copy(NUM(f), ntRing), p_Copy(DEN(f), ntRing), cf->extRing);
1120  if (p_IsConstant(pGcd, ntRing) &&
1121      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
1122  { /* gcd = 1; nothing to cancel;
1123       Suppose the given rational function field is over Q. Although the
1124       gcd is 1, we may have produced fractional coefficients in NUM(f),
1125       DEN(f), or both, due to previous arithmetics. The next call will
1126       remove those nested fractions, in case there are any. */
1127    if (nCoeff_is_Zp(ntCoeffs) && p_IsConstant (DEN (f), ntRing))
1128    {
1129      NUM (f) = p_Div_nn (NUM (f), p_GetCoeff (DEN(f),ntRing), ntRing);
1130      //poly newNum= singclap_pdivide (NUM(f), DEN (f), ntRing);
1131      //p_Delete(&NUM (f), ntRing);
1132      //NUM (f)= newNum;
1133      p_Delete(&DEN (f), ntRing);
1134      DEN (f) = NULL;
1135      COM (f) = 0;
1136    } else if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
1137  }
1138  else
1139  { /* We divide both NUM(f) and DEN(f) by the gcd which is known
1140       to be != 1. */
1141    poly newNum = singclap_pdivide(NUM(f), pGcd, ntRing);
1142    p_Delete(&NUM(f), ntRing);
1143    NUM(f) = newNum;
1144    poly newDen = singclap_pdivide(DEN(f), pGcd, ntRing);
1145    p_Delete(&DEN(f), ntRing);
1146    DEN(f) = newDen;
1147    if (p_IsConstant(DEN(f), ntRing) &&
1148        n_IsOne(p_GetCoeff(DEN(f), ntRing), ntCoeffs))
1149    {
1150      /* DEN(f) = 1 needs to be represented by NULL! */
1151      p_Delete(&DEN(f), ntRing);
1152      newDen = NULL;
1153    }
1154    else
1155    { /* Note that over Q, by cancelling the gcd, we may have produced
1156         fractional coefficients in NUM(f), DEN(f), or both. The next
1157         call will remove those nested fractions, in case there are
1158         any. */
1159      if (nCoeff_is_Q(ntCoeffs)) handleNestedFractionsOverQ(f, cf);
1160    }
1161  }
1162  COM(f) = 0;
1163  p_Delete(&pGcd, ntRing);
1164
1165  if( DEN(f) != NULL )
1166    if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1167    {
1168      NUM(f) = p_Neg(NUM(f), ntRing);
1169      DEN(f) = p_Neg(DEN(f), ntRing);
1170    }
1171#endif /* HAVE_FACTORY */
1172
1173  ntTest(a); // !!!!
1174}
1175
1176// NOTE: modifies a
1177void ntWriteLong(number &a, const coeffs cf)
1178{
1179  ntTest(a);
1180  definiteGcdCancellation(a, cf, FALSE);
1181  if (IS0(a))
1182    StringAppendS("0");
1183  else
1184  {
1185    fraction f = (fraction)a;
1186    // stole logic from napWrite from kernel/longtrans.cc of legacy singular
1187    BOOLEAN omitBrackets = p_IsConstant(NUM(f), ntRing);
1188    if (!omitBrackets) StringAppendS("(");
1189    p_String0Long(NUM(f), ntRing, ntRing);
1190    if (!omitBrackets) StringAppendS(")");
1191    if (!DENIS1(f))
1192    {
1193      StringAppendS("/");
1194      omitBrackets = p_IsConstant(DEN(f), ntRing);
1195      if (!omitBrackets) StringAppendS("(");
1196      p_String0Long(DEN(f), ntRing, ntRing);
1197      if (!omitBrackets) StringAppendS(")");
1198    }
1199  }
1200  ntTest(a); // !!!!
1201}
1202
1203// NOTE: modifies a
1204void ntWriteShort(number &a, const coeffs cf)
1205{
1206  ntTest(a);
1207  definiteGcdCancellation(a, cf, FALSE);
1208  if (IS0(a))
1209    StringAppendS("0");
1210  else
1211  {
1212    fraction f = (fraction)a;
1213    // stole logic from napWrite from kernel/longtrans.cc of legacy singular
1214    BOOLEAN omitBrackets = p_IsConstant(NUM(f), ntRing);
1215    if (!omitBrackets) StringAppendS("(");
1216    p_String0Short(NUM(f), ntRing, ntRing);
1217    if (!omitBrackets) StringAppendS(")");
1218    if (!DENIS1(f))
1219    {
1220      StringAppendS("/");
1221      omitBrackets = p_IsConstant(DEN(f), ntRing);
1222      if (!omitBrackets) StringAppendS("(");
1223      p_String0Short(DEN(f), ntRing, ntRing);
1224      if (!omitBrackets) StringAppendS(")");
1225    }
1226  }
1227  ntTest(a);
1228}
1229
1230const char * ntRead(const char *s, number *a, const coeffs cf)
1231{
1232  poly p;
1233  const char * result = p_Read(s, p, ntRing);
1234  if (p == NULL) *a = NULL;
1235  else *a = ntInit(p, cf);
1236  return result;
1237}
1238
1239void ntNormalize (number &a, const coeffs cf)
1240{
1241  if ((a!=NULL))
1242  {
1243    definiteGcdCancellation(a, cf, FALSE);
1244    if ((DEN(a)!=NULL)
1245    &&(!n_GreaterZero(pGetCoeff(DEN(a)),ntRing->cf)))
1246    {
1247      NUM(a)=p_Neg(NUM(a),ntRing);
1248      DEN(a)=p_Neg(DEN(a),ntRing);
1249    }
1250  }
1251  ntTest(a); // !!!!
1252}
1253
1254/* expects *param to be castable to TransExtInfo */
1255static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
1256{
1257  if (ID != n) return FALSE;
1258  TransExtInfo *e = (TransExtInfo *)param;
1259  /* for rational function fields we expect the underlying
1260     polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
1261     this expectation is based on the assumption that we have properly
1262     registered cf and perform reference counting rather than creating
1263     multiple copies of the same coefficient field/domain/ring */
1264  if (ntRing == e->r)
1265    return TRUE;
1266
1267  // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)...
1268  if( rEqual(ntRing, e->r, TRUE) )
1269  {
1270    rDelete(e->r);
1271    return TRUE;
1272  }
1273
1274  return FALSE;
1275}
1276
1277number ntLcm(number a, number b, const coeffs cf)
1278{
1279  ntTest(a);
1280  ntTest(b);
1281  fraction fb = (fraction)b;
1282  if ((b==NULL)||(DEN(fb)==NULL)) return ntCopy(a,cf);
1283#ifdef HAVE_FACTORY
1284  fraction fa = (fraction)a;
1285  /* singclap_gcd destroys its arguments; we hence need copies: */
1286  poly pa = p_Copy(NUM(fa), ntRing);
1287  poly pb = p_Copy(DEN(fb), ntRing);
1288
1289  /* Note that, over Q, singclap_gcd will remove the denominators in all
1290     rational coefficients of pa and pb, before starting to compute
1291     the gcd. Thus, we do not need to ensure that the coefficients of
1292     pa and pb live in Z; they may well be elements of Q\Z. */
1293  poly pGcd = singclap_gcd(pa, pb, cf->extRing);
1294  if (p_IsConstant(pGcd, ntRing) &&
1295      n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs))
1296  { /* gcd = 1; return pa*pb*/
1297    p_Delete(&pGcd,ntRing);
1298    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1299    NUM(result) = pp_Mult_qq(NUM(fa),DEN(fb),ntRing);
1300
1301    ntTest((number)result); // !!!!
1302
1303    return (number)result;
1304  }
1305
1306
1307  /* return pa*pb/gcd */
1308    poly newNum = singclap_pdivide(NUM(fa), pGcd, ntRing);
1309    p_Delete(&pGcd,ntRing);
1310    fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1311    NUM(result) = p_Mult_q(p_Copy(DEN(fb),ntRing),newNum,ntRing);
1312    ntTest((number)result); // !!!!
1313    return (number)result;
1314
1315#else
1316  Print("// factory needed: transext.cc:ntLcm\n");
1317  return NULL;
1318#endif /* HAVE_FACTORY */
1319  return NULL;
1320}
1321
1322number ntGcd(number a, number b, const coeffs cf)
1323{
1324  ntTest(a);
1325  ntTest(b);
1326  if (a==NULL) return ntCopy(b,cf);
1327  if (b==NULL) return ntCopy(a,cf);
1328#ifdef HAVE_FACTORY
1329  fraction fa = (fraction)a;
1330  fraction fb = (fraction)b;
1331  /* singclap_gcd destroys its arguments; we hence need copies: */
1332  poly pa = p_Copy(NUM(fa), ntRing);
1333  poly pb = p_Copy(NUM(fb), ntRing);
1334
1335  /* Note that, over Q, singclap_gcd will remove the denominators in all
1336     rational coefficients of pa and pb, before starting to compute
1337     the gcd. Thus, we do not need to ensure that the coefficients of
1338     pa and pb live in Z; they may well be elements of Q\Z. */
1339  poly pGcd = singclap_gcd(pa, pb, cf->extRing);
1340  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1341  NUM(result) = pGcd;
1342  ntTest((number)result); // !!!!
1343  return (number)result;
1344#else
1345  Print("// factory needed: transext.cc:ntGcd\n");
1346  return NULL;
1347#endif /* HAVE_FACTORY */
1348}
1349
1350int ntSize(number a, const coeffs cf)
1351{
1352  ntTest(a);
1353  if (IS0(a)) return -1;
1354  /* this has been taken from the old implementation of field extensions,
1355     where we computed the sum of the degrees and the numbers of terms in
1356     the numerator and denominator of a; so we leave it at that, for the
1357     time being */
1358  fraction f = (fraction)a;
1359  poly p = NUM(f);
1360  int noOfTerms = 0;
1361  int numDegree = 0;
1362  while (p != NULL)
1363  {
1364    noOfTerms++;
1365    int d = 0;
1366    for (int i = 1; i <= rVar(ntRing); i++)
1367      d += p_GetExp(p, i, ntRing);
1368    if (d > numDegree) numDegree = d;
1369    pIter(p);
1370  }
1371  int denDegree = 0;
1372  if (!DENIS1(f))
1373  {
1374    p = DEN(f);
1375    while (p != NULL)
1376    {
1377      noOfTerms++;
1378      int d = 0;
1379      for (int i = 1; i <= rVar(ntRing); i++)
1380        d += p_GetExp(p, i, ntRing);
1381      if (d > denDegree) denDegree = d;
1382      pIter(p);
1383    }
1384  }
1385  ntTest(a); // !!!!
1386  return numDegree + denDegree + noOfTerms;
1387}
1388
1389number ntInvers(number a, const coeffs cf)
1390{
1391  ntTest(a);
1392  if (IS0(a))
1393  {
1394    WerrorS(nDivBy0);
1395    return NULL;
1396  }
1397  fraction f = (fraction)a;
1398  assume( f != NULL );
1399
1400  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1401
1402  assume( NUM(f) != NULL );
1403  const poly den = DEN(f);
1404
1405  if (den == NULL)
1406    NUM(result) = p_One(ntRing);
1407  else
1408    NUM(result) = p_Copy(den, ntRing);
1409
1410  if( !NUMIS1(f) )
1411  {
1412    poly num_f=NUM(f);
1413    BOOLEAN neg= !n_GreaterZero(pGetCoeff(num_f),ntRing->cf);
1414    if (neg)
1415    {
1416      num_f=p_Neg(p_Copy(num_f, ntRing), ntRing);
1417      NUM(result)=p_Neg(NUM(result), ntRing);
1418    }
1419    else
1420    {
1421      num_f=p_Copy(num_f, ntRing);
1422    }
1423    DEN(result) = num_f;
1424    COM(result) = COM(f);
1425    if (neg)
1426    {
1427      if (p_IsOne(num_f, ntRing))
1428      {
1429        DEN(result)=NULL;
1430        COM(result) = 0;
1431        p_Delete(&num_f,ntRing);
1432      }
1433    }
1434  }
1435  else
1436  {
1437    DEN(result) = NULL;
1438    COM(result) = 0;
1439  }
1440  ntTest((number)result); // !!!!
1441  return (number)result;
1442}
1443
1444/* assumes that src = Q, dst = Q(t_1, ..., t_s) */
1445number ntMap00(number a, const coeffs src, const coeffs dst)
1446{
1447  if (n_IsZero(a, src)) return NULL;
1448  assume(n_Test(a, src));
1449  assume(src == dst->extRing->cf);
1450  return ntInit(p_NSet(n_Copy(a, src), dst->extRing), dst);
1451}
1452
1453/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
1454number ntMapP0(number a, const coeffs src, const coeffs dst)
1455{
1456  if (n_IsZero(a, src)) return NULL;
1457  assume(n_Test(a, src));
1458  /* mapping via intermediate int: */
1459  int n = n_Int(a, src);
1460  number q = n_Init(n, dst->extRing->cf);
1461  if (n_IsZero(q, dst->extRing->cf))
1462  {
1463    n_Delete(&q, dst->extRing->cf);
1464    return NULL;
1465  }
1466  return ntInit(p_NSet(q, dst->extRing), dst);
1467}
1468
1469/* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or
1470                       src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */
1471number ntCopyMap(number a, const coeffs cf, const coeffs dst)
1472{
1473//  if (n_IsZero(a, cf)) return NULL;
1474
1475  ntTest(a);
1476
1477  if (IS0(a)) return NULL;
1478
1479  const ring rSrc = cf->extRing;
1480  const ring rDst = dst->extRing;
1481
1482  if( rSrc == rDst )
1483    return ntCopy(a, dst); // USUALLY WRONG!
1484
1485  fraction f = (fraction)a;
1486  poly g = prCopyR(NUM(f), rSrc, rDst);
1487
1488  poly h = NULL;
1489
1490  if (!DENIS1(f))
1491     h = prCopyR(DEN(f), rSrc, rDst);
1492
1493  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1494
1495  NUM(result) = g;
1496  DEN(result) = h;
1497  COM(result) = COM(f);
1498  assume(n_Test((number)result, dst));
1499  return (number)result;
1500}
1501
1502number ntCopyAlg(number a, const coeffs cf, const coeffs dst)
1503{
1504  assume( n_Test(a, cf) );
1505  if (n_IsZero(a, cf)) return NULL;
1506
1507  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1508  // DEN(f) = NULL; COM(f) = 0;
1509  NUM(f) = prCopyR((poly)a, cf->extRing, dst->extRing);
1510  assume(n_Test((number)f, dst));
1511  return (number)f;
1512}
1513
1514/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
1515number ntMap0P(number a, const coeffs src, const coeffs dst)
1516{
1517  assume( n_Test(a, src) );
1518  if (n_IsZero(a, src)) return NULL;
1519  int p = rChar(dst->extRing);
1520  number q = nlModP(a, src, dst->extRing->cf);
1521
1522  if (n_IsZero(q, dst->extRing->cf))
1523  {
1524    n_Delete(&q, dst->extRing->cf);
1525    return NULL;
1526  }
1527
1528  poly g = p_NSet(q, dst->extRing);
1529
1530  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1531  NUM(f) = g; // DEN(f) = NULL; COM(f) = 0;
1532  assume(n_Test((number)f, dst));
1533  return (number)f;
1534}
1535
1536/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
1537number ntMapPP(number a, const coeffs src, const coeffs dst)
1538{
1539  assume( n_Test(a, src) );
1540  if (n_IsZero(a, src)) return NULL;
1541  assume(src == dst->extRing->cf);
1542  poly p = p_One(dst->extRing);
1543  p_SetCoeff(p, n_Copy(a, src), dst->extRing);
1544  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1545  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1546  assume(n_Test((number)f, dst));
1547  return (number)f;
1548}
1549
1550/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
1551number ntMapUP(number a, const coeffs src, const coeffs dst)
1552{
1553  assume( n_Test(a, src) );
1554  if (n_IsZero(a, src)) return NULL;
1555  /* mapping via intermediate int: */
1556  int n = n_Int(a, src);
1557  number q = n_Init(n, dst->extRing->cf);
1558  poly p;
1559  if (n_IsZero(q, dst->extRing->cf))
1560  {
1561    n_Delete(&q, dst->extRing->cf);
1562    return NULL;
1563  }
1564  p = p_One(dst->extRing);
1565  p_SetCoeff(p, q, dst->extRing);
1566  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1567  NUM(f) = p; DEN(f) = NULL; COM(f) = 0;
1568  assume(n_Test((number)f, dst));
1569  return (number)f;
1570}
1571
1572nMapFunc ntSetMap(const coeffs src, const coeffs dst)
1573{
1574  /* dst is expected to be a rational function field */
1575  assume(getCoeffType(dst) == ID);
1576
1577  if( src == dst ) return ndCopyMap;
1578
1579  int h = 0; /* the height of the extension tower given by dst */
1580  coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
1581  coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
1582
1583  /* for the time being, we only provide maps if h = 1 and if b is Q or
1584     some field Z/pZ: */
1585  if (h==0)
1586  {
1587    if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst))
1588      return ntMap00;                                 /// Q       -->  Q(T)
1589    if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
1590      return ntMapP0;                                 /// Z/p     -->  Q(T)
1591    if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst))
1592      return ntMap0P;                                 /// Q       --> Z/p(T)
1593    if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
1594    {
1595      if (src->ch == dst->ch) return ntMapPP;         /// Z/p     --> Z/p(T)
1596      else return ntMapUP;                            /// Z/u     --> Z/p(T)
1597    }
1598  }
1599  if (h != 1) return NULL;
1600  if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
1601
1602  /* Let T denote the sequence of transcendental extension variables, i.e.,
1603     K[t_1, ..., t_s] =: K[T];
1604     Let moreover, for any such sequence T, T' denote any subsequence of T
1605     of the form t_1, ..., t_w with w <= s. */
1606
1607  if ((!nCoeff_is_Zp(bSrc)) && (!nCoeff_is_Q(bSrc))) return NULL;
1608
1609  if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst))
1610  {
1611    if (rVar(src->extRing) > rVar(dst->extRing))
1612       return NULL;
1613
1614    for (int i = 0; i < rVar(src->extRing); i++)
1615      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0)
1616         return NULL;
1617
1618    if (src->type==n_transExt)
1619       return ntCopyMap;          /// Q(T')   --> Q(T)
1620    else
1621       return ntCopyAlg;
1622  }
1623
1624  if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst))
1625  {
1626    if (rVar(src->extRing) > rVar(dst->extRing))
1627       return NULL;
1628
1629    for (int i = 0; i < rVar(src->extRing); i++)
1630      if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0)
1631         return NULL;
1632
1633    if (src->type==n_transExt)
1634       return ntCopyMap;         /// Z/p(T') --> Z/p(T)
1635    else
1636       return ntCopyAlg;
1637  }
1638
1639  return NULL;                                 /// default
1640}
1641#if 0
1642nMapFunc ntSetMap_T(const coeffs src, const coeffs dst)
1643{
1644  nMapFunc n=ntSetMap(src,dst);
1645  if (n==ntCopyAlg) printf("n=ntCopyAlg\n");
1646  else if (n==ntCopyMap) printf("n=ntCopyMap\n");
1647  else if (n==ntMapUP) printf("n=ntMapUP\n");
1648  else if (n==ntMap0P) printf("n=ntMap0P\n");
1649  else if (n==ntMapP0) printf("n=ntMapP0\n");
1650  else if (n==ntMap00) printf("n=ntMap00\n");
1651  else if (n==NULL) printf("n=NULL\n");
1652  else printf("n=?\n");
1653  return n;
1654}
1655#endif
1656
1657void ntKillChar(coeffs cf)
1658{
1659  if ((--cf->extRing->ref) == 0)
1660    rDelete(cf->extRing);
1661}
1662#ifdef HAVE_FACTORY
1663number ntConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
1664{
1665  if (n.isZero()) return NULL;
1666  poly p=convFactoryPSingP(n,ntRing);
1667  fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1668  NUM(result) = p;
1669  //DEN(result) = NULL; // done by omAlloc0Bin
1670  //COM(result) = 0; // done by omAlloc0Bin
1671  ntTest((number)result);
1672  return (number)result;
1673}
1674CanonicalForm ntConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs cf )
1675{
1676  ntTest(n);
1677  if (IS0(n)) return CanonicalForm(0);
1678
1679  fraction f = (fraction)n;
1680  return convSingPFactoryP(NUM(f),ntRing);
1681}
1682#endif
1683
1684static int ntParDeg(number a, const coeffs cf)
1685{
1686  ntTest(a);
1687  if (IS0(a)) return -1;
1688  fraction fa = (fraction)a;
1689  return cf->extRing->pFDeg(NUM(fa),cf->extRing);
1690}
1691
1692/// return the specified parameter as a number in the given trans.ext.
1693static number ntParameter(const int iParameter, const coeffs cf)
1694{
1695  assume(getCoeffType(cf) == ID);
1696
1697  const ring R = cf->extRing;
1698  assume( R != NULL );
1699  assume( 0 < iParameter && iParameter <= rVar(R) );
1700
1701  poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
1702
1703//  return (number) p;
1704
1705  fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
1706  NUM(f) = p;
1707  DEN(f) = NULL;
1708  COM(f) = 0;
1709
1710  ntTest((number)f);
1711
1712  return (number)f;
1713}
1714
1715/// if m == var(i)/1 => return i,
1716int ntIsParam(number m, const coeffs cf)
1717{
1718  ntTest(m);
1719  assume(getCoeffType(cf) == ID);
1720
1721  const ring R = cf->extRing;
1722  assume( R != NULL );
1723
1724  fraction f = (fraction)m;
1725
1726  if( DEN(f) != NULL )
1727    return 0;
1728
1729  return p_Var( NUM(f), R );
1730}
1731
1732struct NTNumConverter
1733{
1734  static inline poly convert(const number& n)
1735  {
1736    // suitable for trans. ext. numbers that are fractions of polys
1737    return NUM((fraction)n); // return the numerator
1738  }
1739};
1740
1741
1742static void ntClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1743{
1744  assume(cf != NULL);
1745  assume(getCoeffType(cf) == ID);
1746  // all coeffs are given by fractions of polynomails over integers!!!
1747  // without denominators!!!
1748
1749  const ring   R = cf->extRing;
1750  assume(R != NULL);
1751  const coeffs Q = R->cf;
1752  assume(Q != NULL);
1753  assume(nCoeff_is_Q(Q));
1754
1755
1756  numberCollectionEnumerator.Reset();
1757
1758  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
1759  {
1760    c = ntInit(1, cf);
1761    return;
1762  }
1763
1764  // all coeffs are given by integers after returning from this routine
1765
1766  // part 1, collect product of all denominators /gcds
1767  poly cand = NULL;
1768
1769  do
1770  {
1771    number &n = numberCollectionEnumerator.Current();
1772
1773    ntNormalize(n, cf);
1774
1775    fraction f = (fraction)n;
1776
1777    assume( f != NULL );
1778
1779    const poly den = DEN(f);
1780
1781    assume( den == NULL ); // ?? / 1 ?
1782
1783    const poly num = NUM(f);
1784
1785    if( cand == NULL )
1786      cand = p_Copy(num, R);
1787    else
1788      cand = singclap_gcd(cand, p_Copy(num, R), R); // gcd(cand, num)
1789
1790    if( p_IsConstant(cand, R) )
1791      break;
1792  }
1793  while( numberCollectionEnumerator.MoveNext() ) ;
1794
1795
1796  // part2: all coeffs = all coeffs * cand
1797  if( cand != NULL )
1798  {
1799  if( !p_IsConstant(cand, R) )
1800  {
1801    c = ntInit(cand, cf);
1802    numberCollectionEnumerator.Reset();
1803    while (numberCollectionEnumerator.MoveNext() )
1804    {
1805      number &n = numberCollectionEnumerator.Current();
1806      const number t = ntDiv(n, c, cf); // TODO: rewrite!?
1807      ntDelete(&n, cf);
1808      n = t;
1809    }
1810  } // else NUM (result) = p_One(R);
1811  else { p_Delete(&cand, R); cand = NULL; }
1812  }
1813
1814  // Quick and dirty fix for constant content clearing: consider numerators???
1815  CRecursivePolyCoeffsEnumerator<NTNumConverter> itr(numberCollectionEnumerator); // recursively treat the NUM(numbers) as polys!
1816  number cc;
1817
1818  n_ClearContent(itr, cc, Q);
1819  number g = ntInit(p_NSet(cc, R), cf);
1820
1821  if( cand != NULL )
1822  {
1823    number gg = ntMult(g, c, cf);
1824    ntDelete(&g, cf);
1825    ntDelete(&c, cf); c = gg;
1826  } else
1827    c = g;
1828  ntTest(c);
1829}
1830
1831static void ntClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
1832{
1833  assume(cf != NULL);
1834  assume(getCoeffType(cf) == ID); // both over Q(a) and Zp(a)!
1835  // all coeffs are given by fractions of polynomails over integers!!!
1836
1837  numberCollectionEnumerator.Reset();
1838
1839  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
1840  {
1841    c = ntInit(1, cf);
1842    return;
1843  }
1844
1845  // all coeffs are given by integers after returning from this routine
1846
1847  // part 1, collect product of all denominators /gcds
1848  poly cand = NULL;
1849
1850  const ring R = cf->extRing;
1851  assume(R != NULL);
1852
1853  const coeffs Q = R->cf;
1854  assume(Q != NULL);
1855//  assume(nCoeff_is_Q(Q));
1856
1857  do
1858  {
1859    number &n = numberCollectionEnumerator.Current();
1860
1861    ntNormalize(n, cf);
1862
1863    fraction f = (fraction)ntGetDenom (n, cf);
1864
1865    assume( f != NULL );
1866
1867    const poly den = NUM(f);
1868
1869    if( den == NULL ) // ?? / 1 ?
1870      continue;
1871
1872    if( cand == NULL )
1873      cand = p_Copy(den, R);
1874    else
1875    {
1876      // cand === LCM( cand, den )!!!!
1877      // NOTE: maybe it's better to make the product and clearcontent afterwards!?
1878      // TODO: move the following to factory?
1879      poly gcd = singclap_gcd(p_Copy(cand, R), p_Copy(den, R), R); // gcd(cand, den) is monic no mater leading coeffs! :((((
1880      if (nCoeff_is_Q (Q))
1881      {
1882        number LcGcd= n_Gcd (p_GetCoeff (cand, R), p_GetCoeff(den, R), Q);
1883        gcd = p_Mult_nn(gcd, LcGcd, R);
1884        n_Delete(&LcGcd,Q);
1885      }
1886//      assume( n_IsOne(pGetCoeff(gcd), Q) ); // TODO: this may be wrong...
1887      cand = p_Mult_q(cand, p_Copy(den, R), R); // cand *= den
1888      const poly t = singclap_pdivide( cand, gcd, R ); // cand' * den / gcd(cand', den)
1889      p_Delete(&cand, R);
1890      p_Delete(&gcd, R);
1891      cand = t;
1892    }
1893  }
1894  while( numberCollectionEnumerator.MoveNext() );
1895
1896  if( cand == NULL )
1897  {
1898    c = ntInit(1, cf);
1899    return;
1900  }
1901
1902  c = ntInit(cand, cf);
1903
1904  numberCollectionEnumerator.Reset();
1905
1906  number d = NULL;
1907
1908  while (numberCollectionEnumerator.MoveNext() )
1909  {
1910    number &n = numberCollectionEnumerator.Current();
1911    number t = ntMult(n, c, cf); // TODO: rewrite!?
1912    ntDelete(&n, cf);
1913
1914    ntNormalize(t, cf); // TODO: needed?
1915    n = t;
1916
1917    fraction f = (fraction)t;
1918    assume( f != NULL );
1919
1920    const poly den = DEN(f);
1921
1922    if( den != NULL ) // ?? / ?? ?
1923    {
1924      assume( p_IsConstant(den, R) );
1925      assume( pNext(den) == NULL );
1926
1927      if( d == NULL )
1928        d = n_Copy(pGetCoeff(den), Q);
1929      else
1930      {
1931        number g = n_Lcm(d, pGetCoeff(den), Q);
1932        n_Delete(&d, Q); d = g;
1933      }
1934    }
1935  }
1936
1937  if( d != NULL )
1938  {
1939    numberCollectionEnumerator.Reset();
1940    while (numberCollectionEnumerator.MoveNext() )
1941    {
1942      number &n = numberCollectionEnumerator.Current();
1943      fraction f = (fraction)n;
1944
1945      assume( f != NULL );
1946
1947      const poly den = DEN(f);
1948
1949      if( den == NULL ) // ?? / 1 ?
1950        NUM(f) = p_Mult_nn(NUM(f), d, R);
1951      else
1952      {
1953        assume( p_IsConstant(den, R) );
1954        assume( pNext(den) == NULL );
1955
1956        number ddd = n_Div(d, pGetCoeff(den), Q); // but be an integer now!!!
1957        NUM(f) = p_Mult_nn(NUM(f), ddd, R);
1958        n_Delete(&ddd, Q);
1959
1960        p_Delete(&DEN(f), R);
1961        DEN(f) = NULL; // TODO: check if this is needed!?
1962      }
1963
1964      assume( DEN(f) == NULL );
1965    }
1966
1967    NUM(c) = p_Mult_nn(NUM(c), d, R);
1968    n_Delete(&d, Q);
1969  }
1970
1971
1972  ntTest(c);
1973}
1974
1975BOOLEAN ntInitChar(coeffs cf, void * infoStruct)
1976{
1977
1978  assume( infoStruct != NULL );
1979
1980  TransExtInfo *e = (TransExtInfo *)infoStruct;
1981
1982  assume( e->r                != NULL);      // extRing;
1983  assume( e->r->cf            != NULL);      // extRing->cf;
1984  assume( e->r->qideal == NULL );
1985
1986  assume( cf != NULL );
1987  assume(getCoeffType(cf) == ID);                // coeff type;
1988
1989  ring R = e->r;
1990  assume(R != NULL);
1991
1992  R->ref ++; // increase the ref.counter for the ground poly. ring!
1993
1994  cf->extRing           = R;
1995  /* propagate characteristic up so that it becomes
1996     directly accessible in cf: */
1997  cf->ch = R->cf->ch;
1998  cf->factoryVarOffset = R->cf->factoryVarOffset + rVar(R);
1999
2000  cf->cfGreaterZero  = ntGreaterZero;
2001  cf->cfGreater      = ntGreater;
2002  cf->cfEqual        = ntEqual;
2003  cf->cfIsZero       = ntIsZero;
2004  cf->cfIsOne        = ntIsOne;
2005  cf->cfIsMOne       = ntIsMOne;
2006  cf->cfInit         = ntInit;
2007  cf->cfInit_bigint  = ntInit_bigint;
2008  cf->cfInt          = ntInt;
2009  cf->cfNeg          = ntNeg;
2010  cf->cfAdd          = ntAdd;
2011  cf->cfSub          = ntSub;
2012  cf->cfMult         = ntMult;
2013  cf->cfDiv          = ntDiv;
2014  cf->cfExactDiv     = ntDiv;
2015  cf->cfPower        = ntPower;
2016  cf->cfCopy         = ntCopy;
2017  cf->cfWriteLong    = ntWriteLong;
2018  cf->cfRead         = ntRead;
2019  cf->cfNormalize    = ntNormalize;
2020  cf->cfDelete       = ntDelete;
2021  cf->cfSetMap       = ntSetMap;
2022  cf->cfGetDenom     = ntGetDenom;
2023  cf->cfGetNumerator = ntGetNumerator;
2024  cf->cfRePart       = ntCopy;
2025  cf->cfImPart       = ntImPart;
2026  cf->cfCoeffWrite   = ntCoeffWrite;
2027#ifdef LDEBUG
2028  cf->cfDBTest       = ntDBTest;
2029#endif
2030  cf->cfGcd          = ntGcd;
2031  cf->cfLcm          = ntLcm;
2032  cf->cfSize         = ntSize;
2033  cf->nCoeffIsEqual  = ntCoeffIsEqual;
2034  cf->cfInvers       = ntInvers;
2035  cf->cfIntDiv       = ntDiv;
2036  cf->cfKillChar     = ntKillChar;
2037
2038  if( rCanShortOut(ntRing) )
2039    cf->cfWriteShort = ntWriteShort;
2040  else
2041    cf->cfWriteShort = ntWriteLong;
2042
2043#ifndef HAVE_FACTORY
2044  PrintS("// Warning: The 'factory' module is not available.\n");
2045  PrintS("//          Hence gcd's cannot be cancelled in any\n");
2046  PrintS("//          computed fraction!\n");
2047#else
2048  cf->convFactoryNSingN =ntConvFactoryNSingN;
2049  cf->convSingNFactoryN =ntConvSingNFactoryN;
2050#endif
2051  cf->cfParDeg = ntParDeg;
2052
2053  cf->iNumberOfParameters = rVar(R);
2054  cf->pParameterNames = R->names;
2055  cf->cfParameter = ntParameter;
2056
2057  if( nCoeff_is_Q(R->cf) )
2058    cf->cfClearContent = ntClearContent;
2059
2060  cf->cfClearDenominators = ntClearDenominators;
2061
2062  return FALSE;
2063}
Note: See TracBrowser for help on using the repository browser.