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

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