source: git/libpolys/polys/ext_fields/transext.cc @ 6b32a6b

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