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

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