source: git/libpolys/polys/ext_fields/transext.cc @ 737bbc2

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