source: git/libpolys/polys/ext_fields/transext.cc @ 0736b6

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