My Project
Loading...
Searching...
No Matches
transext.cc
Go to the documentation of this file.
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.
24* For the special case of K = Q (i.e., when computing over the
25* rationals), this definite cancellation procedure will also take
26* care of nested fractions: If there are fractional coefficients
27* in the numerator or denominator of a number, then this number
28* is being replaced by a quotient of two polynomials over Z, or
29* - if the denominator is a constant - by a polynomial over Q.
30*
31* TODO: the description above needs a major update!!!
32*/
33#define TRANSEXT_PRIVATES
34
35#include "misc/auxiliary.h"
36
37#include "factory/factory.h"
38
39#include "reporter/reporter.h"
40
41#include "coeffs/coeffs.h"
42#include "coeffs/numbers.h"
43
44#include "coeffs/longrat.h"
45
48#include "polys/simpleideals.h"
49
50#include "polys/clapsing.h"
51#include "polys/clapconv.h"
52
53#include "polys/prCopy.h"
54#include "transext.h"
55#include "algext.h"
56
58
59
60/* constants for controlling the complexity of numbers */
61#define ADD_COMPLEXITY 1 /**< complexity increase due to + and - */
62#define MULT_COMPLEXITY 2 /**< complexity increase due to * and / */
63#define DIFF_COMPLEXITY 2 /**< complexity increase due to diff */
64#define BOUND_COMPLEXITY 10 /**< maximum complexity of a number */
65
66/// TRUE iff num. represents 1
67#define NUMIS1(f) (p_IsOne(NUM(f), cf->extRing))
68
69#define COM(f) (f)->complexity
70
71
72#ifdef LDEBUG
73static BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs r);
74#endif
75
76#define ntTest(a) n_Test(a, cf)
77
78/* polynomial ring in which the numerators and denominators of our
79 numbers live */
80#define ntRing cf->extRing
81
82/* coeffs object in which the coefficients of our numbers live;
83 * methods attached to ntCoeffs may be used to compute with the
84 * coefficients of our numbers, e.g., use ntCoeffs->nAdd to add
85 * coefficients of our numbers */
86#define ntCoeffs cf->extRing->cf
87
88
89VAR omBin fractionObjectBin = omGetSpecBin(sizeof(fractionObject));
90
91/// forward declarations
92static void heuristicGcdCancellation(number a, const coeffs cf);
93static void definiteGcdCancellation(number a, const coeffs cf,
94 BOOLEAN simpleTestsHaveAlreadyBeenPerformed);
95
96/* test routine, usualy disabled *
97 * if want to activate it, activate also the calls to check_N *
98 *
99void check_normalized(number t,const coeffs cf, const char *f, int l)
100{
101 if (IS0(t)) return;
102 if(rField_is_Q(ntRing))
103 {
104 poly pp=NUM(t);
105 while(pp!=NULL)
106 {
107 if (((SR_HDL(pGetCoeff(pp)) & SR_INT)==0)&&(SR_HDL(pGetCoeff(pp))!=NULL))
108 {
109 if (pGetCoeff(pp)->s==0)
110 {
111 Print("NUM not normalized in %s:%d\n",f,l);
112 p_Normalize(pp,ntRing);
113 }
114 else if (pGetCoeff(pp)->s==1)
115 Print("NUM is rational in %s:%d\n",f,l);
116 }
117 pIter(pp);
118 }
119 pp=DEN(t);
120 while(pp!=NULL)
121 {
122 if (((SR_HDL(pGetCoeff(pp)) & SR_INT)==0)&&(SR_HDL(pGetCoeff(pp))!=NULL))
123 {
124 if (pGetCoeff(pp)->s==0)
125 {
126 Print("NUM not normalized in %s:%d\n",f,l);
127 p_Normalize(pp,ntRing);
128 }
129 else if (pGetCoeff(pp)->s==1)
130 Print("DEN is rational in %s:%d\n",f,l);
131 }
132 pIter(pp);
133 }
134 }
135}
136#define check_N(A,B) check_normalized(A,B,__FILE__,__LINE__)
137*/
138
139#ifdef LDEBUG
140static BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs cf)
141{
143
144 if (IS0(a)) return TRUE;
145
146 const fraction t = (fraction)a;
147
148 //check_N(a,cf);
149 const poly num = NUM(t);
150 assume(num != NULL); ///< t != 0 ==> numerator(t) != 0
151
152 p_Test(num, ntRing);
153
155 for( poly p = num; p != NULL; pIter(p) )
157 {
158 Print("ERROR in %s:%d: non-integer Q coeff in num. poly\n",f,l);
159 Print("TERM: "); p_wrp(p, ntRing); PrintLn();
160 return FALSE;
161 }
162
163 const poly den = DEN(t);
164
165 if (den != NULL) // !DENIS1(f)
166 {
167 p_Test(den, ntRing);
168
170 for( poly p = den; p != NULL; pIter(p) )
172 {
173 Print("ERROR in %s:%d: non-integer Q coeff in den. poly\n",f,l);
174 Print("TERM: "); p_wrp(p, ntRing); PrintLn();
175 return FALSE;
176 }
177
179 {
180 if( p_IsConstant(den, ntRing) )
181 {
182 Print("ERROR in %s:%d: constant den. poly / Zp\n",f,l);
183 PrintS("NUM: "); p_Write(num, ntRing);
184 PrintS("DEN: "); p_Write(den, ntRing);
185 return FALSE;
186 }
187
188 if( !n_IsOne(pGetCoeff(den), ntCoeffs) )
189 {
190 Print("ERROR in %s:%d: non-monic den. poly / Zp\n",f,l);
191 PrintS("NUM: "); p_Write(num, ntRing);
192 PrintS("DEN: "); p_Write(den, ntRing);
193 return FALSE;
194 }
195 }
196
197 if (COM(t)==0)
198 {
199 poly gcd = singclap_gcd_r( num, den, ntRing );
200 if(gcd!=NULL)
201 {
202 if( !p_IsOne(gcd, ntRing) )
203 {
204 Print("ERROR in %s:%d: 1 != GCD between num. & den. poly\n",f,l);
205 PrintS("GCD: "); p_Write(gcd, ntRing);
206 PrintS("NUM: "); p_Write(num, ntRing);
207 PrintS("DEN: "); p_Write(den, ntRing);
208 return FALSE;
209 }
210 p_Delete( &gcd, ntRing );
211 }
212 }
213 return TRUE;
214
216 {
217 Print("?/1 in %s:%d\n",f,l);
218 return FALSE;
219 }
221 {
222 Print("negative sign of DEN. of a fraction in %s:%d\n",f,l);
223 return FALSE;
224 }
225 // test that den is over integers!?
226 }
227 else
228 {
229 return TRUE;
230
231 // num != NULL // den == NULL
232// if( COM(t) != 0 )
233// {
234// Print("?//NULL with non-zero complexity: %d in %s:%d\n", COM(t), f, l);
235// return FALSE;
236// }
237 // test that nume is over integers!?
238 }
240 {
241 poly p=num; // !=NULL
242 do
243 {
244 number n=pGetCoeff(p);
245 n_Test(n,ntCoeffs);
246 if ((!(SR_HDL(n) & SR_INT))&&(n->s==0))
247 /* not normalized, just do for the following test*/
248 {
250 n=pGetCoeff(p);
251 }
252 if (!(SR_HDL(n) & SR_INT))
253 {
254 if (n->s<2)
255 Print("rational coeff in num: %s:%d\n",f,l);
256 }
257 pIter(p);
258 } while(p!=NULL);
259 p=den;
260 while(p!=NULL)
261 {
262 number n=pGetCoeff(p);
263 if (!(SR_HDL(n) & SR_INT))
264 {
265 if (n->s!=3)
266 Print("rational coeff in den.:%s:%d\n",f,l);
267 }
268 pIter(p);
269 }
270 }
271 return TRUE;
272}
273#endif
274
275poly gcd_over_Q ( poly f, poly g, const ring r)
276{
277 poly res;
278 f=p_Copy(f,r);
279 p_Cleardenom(f, r);
280 g=p_Copy(g,r);
281 p_Cleardenom(g, r);
283 p_Delete(&f, r);
284 p_Delete(&g, r);
285 return res;
286}
287
288/* returns the bottom field in this field extension tower; if the tower
289 is flat, i.e., if there is no extension, then r itself is returned;
290 as a side-effect, the counter 'height' is filled with the height of
291 the extension tower (in case the tower is flat, 'height' is zero) */
292static coeffs nCoeff_bottom(const coeffs r, int &height)
293{
294 assume(r != NULL);
295 coeffs cf = r;
296 height = 0;
297 while (nCoeff_is_Extension(cf))
298 {
299 assume(cf->extRing != NULL); assume(cf->extRing->cf != NULL);
300 cf = cf->extRing->cf;
301 height++;
302 }
303 return cf;
304}
305
306static BOOLEAN ntIsZero(number a, const coeffs cf)
307{
308 //check_N(a,cf);
309 ntTest(a); // !!!
310 return (IS0(a));
311}
312
313static void ntDelete(number * a, const coeffs cf)
314{
315 //check_N(*a,cf);
316 ntTest(*a); // !!!
317
318 fraction f = (fraction)(*a);
319 if (IS0(f)) return;
320 p_Delete(&NUM(f), ntRing);
321 if (!DENIS1(f)) p_Delete(&DEN(f), ntRing);
323 *a = NULL;
324}
325
326static BOOLEAN ntEqual(number a, number b, const coeffs cf)
327{
328 //check_N(a,cf);
329 //check_N(b,cf);
330 ntTest(a);
331 ntTest(b);
332
333 /// simple tests
334 if (a == b) return TRUE;
335 if ((IS0(a)) && (!IS0(b))) return FALSE;
336 if ((IS0(b)) && (!IS0(a))) return FALSE;
337
338 /// cheap test if gcd's have been cancelled in both numbers
339 fraction fa = (fraction)a;
340 fraction fb = (fraction)b;
341 if ((COM(fa) == 1) && (COM(fb) == 1))
342 {
343 poly f = p_Add_q(p_Copy(NUM(fa), ntRing),
344 p_Neg(p_Copy(NUM(fb), ntRing), ntRing),
345 ntRing);
346 if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
347 if (DENIS1(fa) && DENIS1(fb)) return TRUE;
348 if (DENIS1(fa) && !DENIS1(fb)) return FALSE;
349 if (!DENIS1(fa) && DENIS1(fb)) return FALSE;
350 f = p_Add_q(p_Copy(DEN(fa), ntRing),
351 p_Neg(p_Copy(DEN(fb), ntRing), ntRing),
352 ntRing);
353 if (f != NULL) { p_Delete(&f, ntRing); return FALSE; }
354 return TRUE;
355 }
356
357 /* default: the more expensive multiplication test
358 a/b = c/d <==> a*d = b*c */
359 poly f = p_Copy(NUM(fa), ntRing);
360 if (!DENIS1(fb)) f = p_Mult_q(f, p_Copy(DEN(fb), ntRing), ntRing);
361 poly g = p_Copy(NUM(fb), ntRing);
362 if (!DENIS1(fa)) g = p_Mult_q(g, p_Copy(DEN(fa), ntRing), ntRing);
363 poly h = p_Add_q(f, p_Neg(g, ntRing), ntRing);
364 if (h == NULL) return TRUE;
365 else
366 {
367 p_Delete(&h, ntRing);
368 return FALSE;
369 }
370}
371
372static number ntCopy(number a, const coeffs cf)
373{
374 //check_N(a,cf);
375 ntTest(a); // !!!
376 if (IS0(a)) return NULL;
377 fraction f = (fraction)a;
378 poly g = NUM(f);
379 poly h = NULL;
380 h =DEN(f);
381 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
382 NUM(result) = p_Copy(g,cf->extRing);
383 DEN(result) = p_Copy(h,cf->extRing);
384 COM(result) = COM(f);
385 ntTest((number)result);
386 return (number)result;
387}
388
389/* assumes that cf represents the rationals, i.e. Q, and will only
390 be called in that case;
391 assumes furthermore that f != NULL and that the denominator of f != 1;
392 generally speaking, this method removes denominators in the rational
393 coefficients of the numerator and denominator of 'a';
394 more concretely, the following normalizations will be performed,
395 where t^alpha denotes a monomial in the transcendental variables t_k
396 (1) if 'a' is of the form
397 (sum_alpha a_alpha/b_alpha * t^alpha)
398 -------------------------------------
399 (sum_beta c_beta/d_beta * t^beta)
400 with integers a_alpha, b_alpha, c_beta, d_beta, then both the
401 numerator and the denominator will be multiplied by the LCM of
402 the b_alpha's and the d_beta's (if this LCM is != 1),
403 (2) if 'a' is - e.g. after having performed step (1) - of the form
404 (sum_alpha a_alpha * t^alpha)
405 -----------------------------
406 (sum_beta c_beta * t^beta)
407 with integers a_alpha, c_beta, and with a non-constant denominator,
408 then both the numerator and the denominator will be divided by the
409 GCD of the a_alpha's and the c_beta's (if this GCD is != 1),
410 this procedure does not alter COM(f) (this has to be done by the
411 calling procedure);
412 modifies f */
413static void handleNestedFractionsOverQ(fraction f, const coeffs cf)
414{
416 assume(!IS0(f));
417 assume(!DENIS1(f));
418
419 { /* step (1); see documentation of this procedure above */
420 number lcmOfDenominators = n_Init(1, ntCoeffs);
421 number c; number tmp;
422 poly p = NUM(f);
423 /* careful when using n_NormalizeHelper!!! It computes the lcm of the numerator
424 of the 1st argument and the denominator of the 2nd!!! */
425 while (p != NULL)
426 {
427 c = p_GetCoeff(p, ntRing);
428 tmp = n_NormalizeHelper(lcmOfDenominators, c, ntCoeffs);
429 n_Delete(&lcmOfDenominators, ntCoeffs);
430 lcmOfDenominators = tmp;
431 pIter(p);
432 }
433 p = DEN(f);
434 while (p != NULL)
435 {
436 c = p_GetCoeff(p, ntRing);
437 tmp = n_NormalizeHelper(lcmOfDenominators, c, ntCoeffs);
438 n_Delete(&lcmOfDenominators, ntCoeffs);
439 lcmOfDenominators = tmp;
440 pIter(p);
441 }
442 if (!n_IsOne(lcmOfDenominators, ntCoeffs))
443 { /* multiply NUM(f) and DEN(f) with lcmOfDenominators */
444 NUM(f) = __p_Mult_nn(NUM(f), lcmOfDenominators, ntRing);
445 p_Normalize(NUM(f), ntRing);
446 DEN(f) = __p_Mult_nn(DEN(f), lcmOfDenominators, ntRing);
447 p_Normalize(DEN(f), ntRing);
448 }
449 n_Delete(&lcmOfDenominators, ntCoeffs);
450 if (DEN(f)!=NULL)
451 { /* step (2); see documentation of this procedure above */
452 p = NUM(f);
453 number gcdOfCoefficients = n_Copy(p_GetCoeff(p, ntRing), ntCoeffs);
454 pIter(p);
455 while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
456 {
457 c = p_GetCoeff(p, ntRing);
458 tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
459 n_Delete(&gcdOfCoefficients, ntCoeffs);
460 gcdOfCoefficients = tmp;
461 pIter(p);
462 }
463 p = DEN(f);
464 while ((p != NULL) && (!n_IsOne(gcdOfCoefficients, ntCoeffs)))
465 {
466 c = p_GetCoeff(p, ntRing);
467 tmp = n_Gcd(c, gcdOfCoefficients, ntCoeffs);
468 n_Delete(&gcdOfCoefficients, ntCoeffs);
469 gcdOfCoefficients = tmp;
470 pIter(p);
471 }
472 if (!n_IsOne(gcdOfCoefficients, ntCoeffs))
473 { /* divide NUM(f) and DEN(f) by gcdOfCoefficients */
474 number inverseOfGcdOfCoefficients = n_Invers(gcdOfCoefficients,
475 ntCoeffs);
476 NUM(f) = __p_Mult_nn(NUM(f), inverseOfGcdOfCoefficients, ntRing);
477 p_Normalize(NUM(f), ntRing);
478 DEN(f) = __p_Mult_nn(DEN(f), inverseOfGcdOfCoefficients, ntRing);
479 p_Normalize(DEN(f), ntRing);
480 n_Delete(&inverseOfGcdOfCoefficients, ntCoeffs);
481 }
482 n_Delete(&gcdOfCoefficients, ntCoeffs);
483 }
484 }
485
486 /* Now, due to the above computations, DEN(f) may have become the
487 1-polynomial which needs to be represented by NULL: */
488 if ((DEN(f) != NULL) &&
489 p_IsConstant(DEN(f), ntRing) &&
491 {
492 p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
493 }
494
495 if( DEN(f) != NULL )
496 if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
497 {
498 NUM(f) = p_Neg(NUM(f), ntRing);
499 DEN(f) = p_Neg(DEN(f), ntRing);
500 }
502 ntTest((number)f); // TODO!
503}
504
505/// TODO: normalization of a!?
506static number ntGetNumerator(number &a, const coeffs cf)
507{
508 //check_N(a,cf);
509 ntTest(a);
510 if (IS0(a)) return NULL;
511
513
514 fraction f = (fraction)a;
515 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
516
517 const BOOLEAN denis1= DENIS1 (f);
518
519 if (getCoeffType (ntCoeffs) == n_Q && !denis1)
521
522 if (getCoeffType (ntCoeffs) == n_Q && denis1)
523 {
524 assume( DEN (f) == NULL );
525
526 number g;
527 // TODO/NOTE: the following should not be necessary (due to
528 // Hannes!) as NUM (f) should be over Z!!!
529 CPolyCoeffsEnumerator itr(NUM(f));
530
531
533
534 if( !n_GreaterZero(g, ntCoeffs) )
535 {
536 NUM (f) = p_Neg(NUM (f), ntRing);
537 g = n_InpNeg(g, ntCoeffs);
538 }
539
540 // g should be a positive integer now!
542
543 if( !n_IsOne(g, ntCoeffs) )
544 {
545 DEN (f) = p_NSet(g, ntRing);
546 COM (f) ++;
547 assume( DEN (f) != NULL );
548 }
549 else
551
552 ntTest(a);
553 }
554
555 // Call ntNormalize instead of above?!?
556
557 NUM (result) = p_Copy (NUM (f), ntRing); // ???
558 //DEN (result) = NULL; // done by ..Alloc0..
559 //COM (result) = 0; // done by ..Alloc0..
560
561 ntTest((number)result);
562 //check_N((number)result,cf);
563 return (number)result;
564}
565
566/// TODO: normalization of a!?
567static number ntGetDenom(number &a, const coeffs cf)
568{
569 //check_N(a,cf);
570 ntTest(a);
571
572 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
573 //DEN (result)= NULL; // done by ..Alloc0..
574 //COM (result)= 0; // done by ..Alloc0..
575
576 if (IS0(a))
577 {
578 NUM (result) = p_One(ntRing);
579 return (number)result;
580 }
581
583
584 fraction f = (fraction)a;
585
586 assume( !IS0(f) );
587
588 const BOOLEAN denis1 = DENIS1 (f);
589
590 if( denis1 && (getCoeffType (ntCoeffs) != n_Q) ) // */1 or 0
591 {
592 NUM (result)= p_One(ntRing);
593 ntTest((number)result);
594 return (number)result;
595 }
596
597 if (!denis1) // */* / Q
598 {
599 assume( DEN (f) != NULL );
600
601 if (getCoeffType (ntCoeffs) == n_Q)
603
604 ntTest(a);
605
606 if( DEN (f) != NULL ) // is it ?? // 1 now???
607 {
608 assume( !p_IsOne(DEN (f), ntRing) );
609
610 NUM (result) = p_Copy (DEN (f), ntRing);
611 ntTest((number)result);
612 return (number)result;
613 }
614// NUM (result) = p_One(ntRing); // NOTE: just in order to be sure...
615 }
616
617 // */1 / Q
619 assume( DEN (f) == NULL );
620
621 number g;
622// poly num= p_Copy (NUM (f), ntRing); // ???
623
624
625 // TODO/NOTE: the following should not be necessary (due to
626 // Hannes!) as NUM (f) should be over Z!!!
627 CPolyCoeffsEnumerator itr(NUM(f));
628
629 n_ClearDenominators(itr, g, ntCoeffs); // may return -1 :(((
630
631 if( !n_GreaterZero(g, ntCoeffs) )
632 {
633// NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
634// g = n_InpNeg(g, ntCoeffs);
635 NUM (f) = p_Neg(NUM (f), ntRing); // Ugly :(((
636 g = n_InpNeg(g, ntCoeffs);
637 }
638
639 // g should be a positive integer now!
641
642 if( !n_IsOne(g, ntCoeffs) )
643 {
645 assume( !n_IsOne(g, ntCoeffs) );
646
647 DEN (f) = p_NSet(g, ntRing); // update COM(f)???
648 assume( DEN (f) != NULL );
649 COM (f) ++;
650
651 NUM (result)= p_Copy (DEN (f), ntRing);
652 }
653 else
654 { // common denom == 1?
655 NUM (result)= p_NSet(g, ntRing); // p_Copy (DEN (f), ntRing);
656// n_Delete(&g, ntCoeffs);
657 }
658
659// if (!p_IsConstant (num, ntRing) && pNext(num) != NULL)
660// else
661// g= p_GetAllDenom (num, ntRing);
662// result= (fraction) ntSetMap (ntCoeffs, cf) (g, ntCoeffs, cf);
663
664 ntTest((number)result);
665 //check_N((number)result,cf);
666 return (number)result;
667}
668
669static BOOLEAN ntIsOne(number a, const coeffs cf)
670{
671 //check_N(a,cf);
672 ntTest(a); // !!!
674 fraction f = (fraction)a;
675 return (f!=NULL) && DENIS1(f) && NUMIS1(f);
676}
677
678static BOOLEAN ntIsMOne(number a, const coeffs cf)
679{
680 //check_N(a,cf);
681 ntTest(a);
683 fraction f = (fraction)a;
684 if ((f==NULL) || (!DENIS1(f))) return FALSE;
685 poly g = NUM(f);
686 if (!p_IsConstant(g, ntRing)) return FALSE;
688}
689
690/// this is in-place, modifies a
691static number ntNeg(number a, const coeffs cf)
692{
693 //check_N(a,cf);
694 ntTest(a);
695 if (!IS0(a))
696 {
697 fraction f = (fraction)a;
698 NUM(f) = p_Neg(NUM(f), ntRing);
699 }
700 ntTest(a);
701 return a;
702}
703
704number ntInit(long i, const coeffs cf)
705{
706 if (i != 0)
707 {
708 poly p=p_ISet(i, ntRing);
709 if (p!=NULL)
710 {
711 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
712 NUM(result) = p;
713 //DEN(result) = NULL; // done by omAlloc0Bin
714 //COM(result) = 0; // done by omAlloc0Bin
715 ntTest((number)result);
716 //check_N((number)result,cf);
717 return (number)result;
718 }
719 }
720 return NULL;
721}
722
723
724// takes over p!
725number ntInit(poly p, const coeffs cf)
726{
727 if (p == NULL) return NULL;
728
729 p_Test( p, ntRing);
730 fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
731
733 {
734 number g;
735 // the following is necessary because
736 // NUM (f) should be over Z,
737 // while p may be over Q
739
741
742 if( !n_GreaterZero(g, ntCoeffs) )
743 {
744 p = p_Neg(p, ntRing);
745 g = n_InpNeg(g, ntCoeffs);
746 }
747
748 // g should be a positive integer now!
750
751 if( !n_IsOne(g, ntCoeffs) )
752 {
753 DEN (f) = p_NSet(g, ntRing);
754 p_Normalize(DEN(f), ntRing);
755 assume( DEN (f) != NULL );
756 }
757 else
758 {
759 //DEN(f) = NULL; // done by omAlloc0
761 }
762 }
763
765 NUM(f) = p;
766 //COM(f) = 0; // done by omAlloc0
767
768 //check_N((number)f,cf);
769 ntTest((number)f);
770 return (number)f;
771}
772
773static long ntInt(number &a, const coeffs cf)
774{
775 //check_N(a,cf);
776 ntTest(a);
777 if (IS0(a)) return 0;
779 fraction f = (fraction)a;
780 if (!DENIS1(f)) return 0;
781
782 const poly aAsPoly = NUM(f);
783
784 if(aAsPoly == NULL)
785 return 0;
786
787 if (!p_IsConstant(aAsPoly, ntRing))
788 return 0;
789
790 assume( aAsPoly != NULL );
791
792 return n_Int(p_GetCoeff(aAsPoly, ntRing), ntCoeffs);
793}
794
795/* return FALSE, if a is a constant <=0 */
796static BOOLEAN ntGreaterZero(number a, const coeffs cf)
797{
798 //check_N(a,cf);
799 ntTest(a);
800 if (IS0(a)) return FALSE;
801 fraction f = (fraction)a;
802 poly g = NUM(f);
804}
805
806static BOOLEAN ntGreater(number a, number b, const coeffs cf)
807{
808 //check_N(a,cf);
809 //check_N(b,cf);
810 ntTest(a);
811 ntTest(b);
812 if (IS0(a))
813 {
814 if (IS0(b)) return FALSE;
815 fraction fb = (fraction)b;
816 return (!n_GreaterZero(pGetCoeff(NUM(fb)), ntCoeffs));
817 }
818 if (IS0(b))
819 {
820 fraction fa = (fraction)a;
821 return (n_GreaterZero(pGetCoeff(NUM(fa)), ntCoeffs));
822 }
823 // now: a!=0, b!=0
824 fraction fa = (fraction)a;
825 number aNumCoeff = p_GetCoeff(NUM(fa), ntRing);
826 int aNumDeg = p_Totaldegree(NUM(fa), ntRing);
827 number aDenCoeff = NULL; int aDenDeg = 0;
828 if (DEN(fa)!=NULL)
829 {
830 aDenCoeff=p_GetCoeff(DEN(fa),ntRing);
831 aDenDeg = p_Totaldegree(DEN(fa), ntRing);
832 }
833 fraction fb = (fraction)b;
834 number bNumCoeff = p_GetCoeff(NUM(fb), ntRing);
835 int bNumDeg = p_Totaldegree(NUM(fb), ntRing);
836 number bDenCoeff = NULL; int bDenDeg = 0;
837 if (DEN(fb)!=NULL)
838 {
839 bDenCoeff=p_GetCoeff(DEN(fb),ntRing);
840 bDenDeg = p_Totaldegree(DEN(fb), ntRing);
841 }
842 if (aNumDeg-aDenDeg > bNumDeg-bDenDeg) return TRUE;
843 if (aNumDeg-aDenDeg < bNumDeg-bDenDeg) return FALSE;
844 number aa;
845 number bb;
846 if (bDenCoeff==NULL) aa=n_Copy(aNumCoeff,ntCoeffs);
847 else aa=n_Mult(aNumCoeff,bDenCoeff,ntCoeffs);
848 if (aDenCoeff==NULL) bb=n_Copy(bNumCoeff,ntCoeffs);
849 else bb=n_Mult(bNumCoeff,aDenCoeff,ntCoeffs);
850 BOOLEAN rr= n_Greater(aa, bb, ntCoeffs);
851 n_Delete(&aa,ntCoeffs);
852 n_Delete(&bb,ntCoeffs);
853 return rr;
854}
855
856static void ntCoeffWrite(const coeffs cf, BOOLEAN details)
857{
858 assume( cf != NULL );
859
860 const ring A = cf->extRing;
861
862 assume( A != NULL );
863 assume( A->cf != NULL );
864
865 n_CoeffWrite(A->cf, details);
866
867// rWrite(A);
868
869 const int P = rVar(A);
870 assume( P > 0 );
871
872 PrintS("(");
873
874 for (int nop=0; nop < P; nop ++)
875 {
876 Print("%s", rRingVar(nop, A));
877 if (nop!=P-1) PrintS(", ");
878 }
879
880 PrintS(")");
881
882 assume( A->qideal == NULL );
883
884/*
885 PrintS("// Coefficients live in the rational function field\n");
886 Print("// K(");
887 for (int i = 0; i < rVar(ntRing); i++)
888 {
889 if (i > 0) PrintS(" ");
890 Print("%s", rRingVar(i, ntRing));
891 }
892 PrintS(") with\n");
893 PrintS("// K: "); n_CoeffWrite(cf->extRing->cf);
894*/
895}
896
897number ntDiff(number a, number d, const coeffs cf)
898{
899 //check_N(a,cf);
900 //check_N(d,cf);
901 ntTest(a);
902 ntTest(d);
903
904 if (IS0(d))
905 {
906 WerrorS("ringvar expected");
907 return NULL;
908 }
909 fraction t = (fraction) d;
910 if (!DENIS1(t))
911 {
912 WerrorS("expected differentiation by a variable");
913 return NULL;
914 }
915 int k=p_Var(NUM(t),ntRing);
916 if (k==0)
917 {
918 WerrorS("expected differentiation by a variable");
919 return NULL;
920 }
921
922 if (IS0(a)) return ntCopy(a, cf);
923
924 fraction fa = (fraction)a;
925 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
926 if (DENIS1(fa))
927 {
928 NUM(result) = p_Diff(NUM(fa),k,ntRing);
929 //DEN(result) = NULL; // done by ..Alloc0..
930 if (NUM(result)==NULL)
931 {
933 return(NULL);
934 }
936 //check_N((number)result,cf);
937 ntTest((number)result);
938 return (number)result;
939 }
940
941 poly fg = p_Mult_q(p_Copy(DEN(fa),ntRing),p_Diff(NUM(fa),k,ntRing),ntRing);
942 poly gf = p_Mult_q(p_Copy(NUM(fa),ntRing),p_Diff(DEN(fa),k,ntRing),ntRing);
943 NUM(result) = p_Sub(fg,gf,ntRing);
944 if (NUM(result)==NULL) return(NULL);
945 DEN(result) = pp_Mult_qq(DEN(fa), DEN(fa), ntRing);
948
949 //check_N((number)result,cf);
950 ntTest((number)result);
951 return (number)result;
952}
953
954static number ntAdd(number a, number b, const coeffs cf)
955{
956 //check_N(a,cf);
957 //check_N(b,cf);
958 ntTest(a);
959 ntTest(b);
960 if (IS0(a)) return ntCopy(b, cf);
961 if (IS0(b)) return ntCopy(a, cf);
962
963 fraction fa = (fraction)a;
964 fraction fb = (fraction)b;
965
966 poly g = p_Copy(NUM(fa), ntRing);
967 if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
968 poly h = p_Copy(NUM(fb), ntRing);
969 if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
970 g = p_Add_q(g, h, ntRing);
971
972 if (g == NULL) return NULL;
973
974 poly f;
975 if (DENIS1(fa) && DENIS1(fb)) f = NULL;
976 else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
977 else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
978 else /* both denom's are != 1 */ f = p_Mult_q(p_Copy(DEN(fa), ntRing),
979 p_Copy(DEN(fb), ntRing),
980 ntRing);
981
982 fraction result = (fraction)omAllocBin(fractionObjectBin);
983 NUM(result) = g;
984 DEN(result) = f;
987
988// ntTest((number)result);
989
990 //check_N((number)result,cf);
991 ntTest((number)result);
992 return (number)result;
993}
994
995static void ntInpAdd(number &a, number b, const coeffs cf)
996{
997 ntTest(a);
998 ntTest(b);
999 if (IS0(b)) return;
1000 if (IS0(a)) { a=ntCopy(b, cf); return;}
1001
1002 fraction fa = (fraction)a;
1003 fraction fb = (fraction)b;
1004
1005 if (DENIS1(fa) && DENIS1(fb))
1006 {
1007 poly g = p_Add_q(NUM(fa), p_Copy(NUM(fb),ntRing), ntRing);
1008 NUM(fa)=g;
1009 if (g==NULL) { omFreeBin((ADDRESS)a, fractionObjectBin); a=NULL; }
1010 else heuristicGcdCancellation((number)fa, cf);
1011 return;
1012 }
1013 poly g = NUM(fa);
1014 if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
1015 poly h = p_Copy(NUM(fb), ntRing);
1016 if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
1017 g = p_Add_q(g, h, ntRing);
1018 if (g==NULL) { omFreeBin((ADDRESS)a, fractionObjectBin); a=NULL; return;}
1019 poly f;
1020 if (!DENIS1(fa) && DENIS1(fb)) f = DEN(fa);
1021 else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
1022 else /* both denom's are != 1 */ f = p_Mult_q(DEN(fa),
1023 p_Copy(DEN(fb), ntRing),
1024 ntRing);
1025
1026 NUM(fa) = g;
1027 DEN(fa) = f;
1028 COM(fa) += COM(fb) + ADD_COMPLEXITY;
1030
1031 ntTest(a);
1032}
1033
1034static number ntSub(number a, number b, const coeffs cf)
1035{
1036 //check_N(a,cf);
1037 //check_N(b,cf);
1038 ntTest(a);
1039 ntTest(b);
1040 if (IS0(a)) return ntNeg(ntCopy(b, cf), cf);
1041 if (IS0(b)) return ntCopy(a, cf);
1042
1043 fraction fa = (fraction)a;
1044 fraction fb = (fraction)b;
1045
1046 poly g = p_Copy(NUM(fa), ntRing);
1047 if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
1048 poly h = p_Copy(NUM(fb), ntRing);
1049 if (!DENIS1(fa)) h = p_Mult_q(h, p_Copy(DEN(fa), ntRing), ntRing);
1050 g = p_Add_q(g, p_Neg(h, ntRing), ntRing);
1051
1052 if (g == NULL) return NULL;
1053
1054 poly f;
1055 if (DENIS1(fa) && DENIS1(fb)) f = NULL;
1056 else if (!DENIS1(fa) && DENIS1(fb)) f = p_Copy(DEN(fa), ntRing);
1057 else if (DENIS1(fa) && !DENIS1(fb)) f = p_Copy(DEN(fb), ntRing);
1058 else /* both den's are != 1 */ f = p_Mult_q(p_Copy(DEN(fa), ntRing),
1059 p_Copy(DEN(fb), ntRing),
1060 ntRing);
1061
1062 fraction result = (fraction)omAllocBin(fractionObjectBin);
1063 NUM(result) = g;
1064 DEN(result) = f;
1065 COM(result) = COM(fa) + COM(fb) + ADD_COMPLEXITY;
1067// ntTest((number)result);
1068 //check_N((number)result,cf);
1069 ntTest((number)result);
1070 return (number)result;
1071}
1072
1073static number ntMult(number a, number b, const coeffs cf)
1074{
1075 //check_N(a,cf);
1076 //check_N(b,cf);
1077 ntTest(a); // !!!?
1078 ntTest(b); // !!!?
1079
1080 if (IS0(a) || IS0(b)) return NULL;
1081
1082 fraction fa = (fraction)a;
1083 fraction fb = (fraction)b;
1084
1085 const poly g = pp_Mult_qq(NUM(fa), NUM(fb), ntRing);
1086
1087 if (g == NULL) return NULL; // may happen due to zero divisors???
1088
1089 fraction result = (fraction)omAllocBin(fractionObjectBin);
1090
1091 NUM(result) = g;
1092
1093 const poly da = DEN(fa);
1094 const poly db = DEN(fb);
1095
1096
1097 //check_N((number)result,cf);
1098 if (db == NULL)
1099 {
1100 // b = ? // NULL
1101
1102 if(da == NULL)
1103 { // both fa && fb are ?? // NULL!
1104 DEN(result) = NULL;
1105 COM(result) = 0;
1107 }
1108 else
1109 {
1110 DEN(result) = p_Copy(da, ntRing);
1113 //check_N((number)result,cf);
1114 }
1115 }
1116 else
1117 { // b = ?? / ??
1118 if (da == NULL)
1119 { // a == ? // NULL
1120 DEN(result) = p_Copy(db, ntRing);
1123 //check_N((number)result,cf);
1124 }
1125 else /* both den's are != 1 */
1126 {
1127 DEN(result) = pp_Mult_qq(da, db, ntRing);
1130 //check_N((number)result,cf);
1131 }
1132 }
1133
1134// ntTest((number)result);
1135
1136 //check_N((number)result,cf);
1137 ntTest((number)result);
1138 return (number)result;
1139}
1140
1141static void ntInpMult(number &a, number b, const coeffs cf)
1142{
1143 ntTest(a); // !!!?
1144 ntTest(b); // !!!?
1145
1146 if (IS0(a)) return;
1147 if (IS0(b)) { omFreeBin((ADDRESS)a, fractionObjectBin); a=NULL; return; }
1148
1149 fraction fa = (fraction)a;
1150 fraction fb = (fraction)b;
1151
1152 const poly g = p_Mult_q(NUM(fa), p_Copy(NUM(fb),ntRing), ntRing);
1153
1154 if (g == NULL) // may happen due to zero divisors???
1155 { omFreeBin((ADDRESS)a, fractionObjectBin); a=NULL; return; }
1156
1157 NUM(fa) = g;
1158
1159 const poly da = DEN(fa);
1160 const poly db = DEN(fb);
1161
1162
1163 //check_N((number)result,cf);
1164 if (db == NULL)
1165 {
1166 // b = ? // NULL
1167
1168 if(da == NULL)
1169 { // both fa && fb are ?? // NULL!
1170 DEN(fa) = NULL;
1171 COM(fa) = 0;
1173 }
1174 else
1175 {
1176 //DEN(fa) = da;
1178 heuristicGcdCancellation((number)fa, cf);
1179 }
1180 }
1181 else
1182 { // b = ?? / ??
1183 if (da == NULL)
1184 { // a == ? // NULL
1185 DEN(fa) = p_Copy(db, ntRing);
1186 COM(fa) = COM(fb) + MULT_COMPLEXITY;
1187 heuristicGcdCancellation((number)fa, cf);
1188 }
1189 else /* both den's are != 1 */
1190 {
1191 DEN(fa) = p_Mult_q(da, p_Copy(db,ntRing), ntRing);
1192 COM(fa) += COM(fb) + MULT_COMPLEXITY;
1193 heuristicGcdCancellation((number)fa, cf);
1194 }
1195 }
1196
1197 ntTest((number)fa);
1198}
1199
1200static void ntNormalizeDen(fraction result, const ring R)
1201{
1202 if ((nCoeff_has_simple_inverse(R->cf))
1203 && (result!=NULL)
1204 && (DEN(result)!=NULL))
1205 {
1206 poly n=DEN(result);
1207 if (!n_IsOne(pGetCoeff(n),R->cf))
1208 {
1209 number inv=n_Invers(pGetCoeff(n),R->cf);
1210 DEN(result)=__p_Mult_nn(n,inv,R);
1211 NUM(result)=__p_Mult_nn(NUM(result),inv,R);
1212 n_Delete(&inv,R->cf);
1213 if (p_IsOne(DEN(result), R))
1214 {
1215 n=DEN(result);
1216 DEN(result)=NULL;
1217 COM(result) = 0;
1218 p_Delete(&n,R);
1219 }
1220 }
1221 }
1222}
1223
1224static number ntDiv(number a, number b, const coeffs cf)
1225{
1226 //check_N(a,cf);
1227 //check_N(b,cf);
1228 ntTest(a);
1229 ntTest(b);
1230 if (IS0(a)) return NULL;
1231 if (IS0(b)) WerrorS(nDivBy0);
1232
1233 fraction fa = (fraction)a;
1234 fraction fb = (fraction)b;
1235
1236 poly g = p_Copy(NUM(fa), ntRing);
1237 if (!DENIS1(fb)) g = p_Mult_q(g, p_Copy(DEN(fb), ntRing), ntRing);
1238
1239 if (g == NULL) return NULL; /* may happen due to zero divisors */
1240
1241 poly f = p_Copy(NUM(fb), ntRing);
1242 if (!DENIS1(fa)) f = p_Mult_q(f, p_Copy(DEN(fa), ntRing), ntRing);
1243
1244 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1245 NUM(result) = g;
1247 {
1248 g=p_Neg(g,ntRing);
1249 f=p_Neg(f,ntRing);
1250 NUM(result) = g;
1251 }
1253 {
1254 DEN(result) = f;
1255 }
1256 else
1257 {
1258 p_Delete(&f, ntRing);
1259 }
1261// definiteGcdCancellation((number)result, cf,FALSE);
1263// ntTest((number)result);
1264 //check_N((number)result,cf);
1266 ntTest((number)result);
1267 return (number)result;
1268}
1269
1270static number ntInvers(number a, const coeffs cf)
1271{
1272 //check_N(a,cf);
1273 ntTest(a);
1274 if (IS0(a))
1275 {
1277 return NULL;
1278 }
1279 fraction f = (fraction)a;
1280 assume( f != NULL );
1281
1282 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1283
1284 assume( NUM(f) != NULL );
1285 const poly den = DEN(f);
1286
1287 if (den == NULL)
1288 NUM(result) = p_One(ntRing);
1289 else
1290 NUM(result) = p_Copy(den, ntRing);
1291
1292 if( !NUMIS1(f) )
1293 {
1294 poly num_f=NUM(f);
1296 if (neg)
1297 {
1298 num_f=p_Neg(p_Copy(num_f, ntRing), ntRing);
1299 NUM(result)=p_Neg(NUM(result), ntRing);
1300 }
1301 else
1302 {
1303 num_f=p_Copy(num_f, ntRing);
1304 }
1305 DEN(result) = num_f;
1306 COM(result) = COM(f);
1307 if (neg)
1308 {
1309 if (p_IsOne(num_f, ntRing))
1310 {
1311 DEN(result)=NULL;
1312 //COM(result) = 0;
1313 p_Delete(&num_f,ntRing);
1314 }
1315 }
1316 }
1317 //else// Alloc0
1318 //{
1319 // DEN(result) = NULL;
1320 // COM(result) = 0;
1321 //}
1323 ntTest((number)result); // !!!!
1324 //check_N((number)result,cf);
1325 return (number)result;
1326}
1327
1328/* 0^0 = 0;
1329 for |exp| <= 7 compute power by a simple multiplication loop;
1330 for |exp| >= 8 compute power along binary presentation of |exp|, e.g.
1331 p^13 = p^1 * p^4 * p^8, where we utilise that
1332 p^(2^(k+1)) = p^(2^k) * p^(2^k);
1333 intermediate cancellation is controlled by the in-place method
1334 heuristicGcdCancellation; see there.
1335*/
1336static void ntPower(number a, int exp, number *b, const coeffs cf)
1337{
1338 ntTest(a);
1339
1340 /* special cases first */
1341 if (IS0(a))
1342 {
1343 if (exp >= 0) *b = NULL;
1344 else WerrorS(nDivBy0);
1345 }
1346 else if (exp == 0) { *b = ntInit(1, cf); return;}
1347 else if (exp == 1) { *b = ntCopy(a, cf); return;}
1348 else if (exp == -1) { *b = ntInvers(a, cf); return;}
1349
1350 int expAbs = exp; if (expAbs < 0) expAbs = -expAbs;
1351
1352 /* now compute a^expAbs */
1353 number pow; number t;
1354 if (expAbs <= 7)
1355 {
1356 pow = ntCopy(a, cf);
1357 for (int i = 2; i <= expAbs; i++)
1358 {
1359 ntInpMult(pow, a, cf);
1361 }
1362 }
1363 else
1364 {
1365 pow = ntInit(1, cf);
1366 number factor = ntCopy(a, cf);
1367 while (expAbs != 0)
1368 {
1369 if (expAbs & 1)
1370 {
1373 }
1374 expAbs = expAbs / 2;
1375 if (expAbs != 0)
1376 {
1379 }
1380 }
1381 ntDelete(&factor, cf);
1382 }
1383
1384 /* invert if original exponent was negative */
1385 if (exp < 0)
1386 {
1387 t = ntInvers(pow, cf);
1388 ntDelete(&pow, cf);
1389 pow = t;
1390 }
1391 *b = pow;
1392 ntTest(*b);
1393 //check_N(*b,cf);
1394}
1395
1396/* modifies a */
1397/* this is an intermediate simplification routine - not a comple "normalize" */
1398static void heuristicGcdCancellation(number a, const coeffs cf)
1399{
1400 if (IS0(a)) return;
1401
1402 fraction f = (fraction)a;
1403 p_Normalize(NUM(f),ntRing);
1404 if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; return; }
1405
1406 assume( DEN(f) != NULL );
1407 p_Normalize(DEN(f),ntRing);
1408
1409 /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
1410 if (p_EqualPolys(NUM(f), DEN(f), ntRing))
1411 { /* numerator and denominator are both != 1 */
1412 p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
1413 p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1414 COM(f) = 0;
1415 }
1416 else
1417 {
1418 if (COM(f) > BOUND_COMPLEXITY)
1420
1421 // TODO: check if it is enough to put the following into definiteGcdCancellation?!
1422 if( DEN(f) != NULL )
1423 {
1424 if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1425 {
1426 NUM(f) = p_Neg(NUM(f), ntRing);
1427 DEN(f) = p_Neg(DEN(f), ntRing);
1428 }
1429 if (ntCoeffs->has_simple_Inverse)
1430 {
1431 if (!n_IsOne(pGetCoeff(DEN(f)),ntCoeffs))
1432 {
1433 number inv=n_Invers(pGetCoeff(DEN(f)),ntCoeffs);
1434 DEN(f)=__p_Mult_nn(DEN(f),inv,ntRing);
1435 NUM(f)=__p_Mult_nn(NUM(f),inv,ntRing);
1436 }
1437 if(p_LmIsConstant(DEN(f),ntRing))
1438 {
1439 p_Delete(&DEN(f),ntRing);
1440 COM(f)=0;
1441 }
1442 }
1443 if ((DEN(f)!=NULL)
1444 && (pNext(DEN(f))==NULL))
1445 {
1446 poly den_f=DEN(f);
1447 poly h=NUM(f);
1448 loop
1449 {
1450 if (h==NULL)
1451 {
1452 h=NUM(f);
1453 do
1454 {
1455 p_ExpVectorDiff(h,h,den_f,ntRing);
1456 pIter(h);
1457 } while(h!=NULL);
1458 p_ExpVectorDiff(den_f,den_f,den_f,ntRing);
1459 break;
1460 }
1461 int i=0;
1462 do
1463 {
1464 i++;
1465 if (p_GetExp(den_f,i,ntRing) > p_GetExp(h,i,ntRing)) return;
1466 } while(i<ntRing->N);
1467 pIter(h);
1468 }
1469 }
1470 }
1471 }
1472 if ((DEN(f)!=NULL)
1473 && (pNext(DEN(f))==NULL)
1474 && (p_LmIsConstantComp(DEN(f),ntRing))
1475 && (n_IsOne(pGetCoeff(DEN(f)),ntCoeffs)))
1476 {
1477 p_Delete(&DEN(f),ntRing);
1478 COM(f)=0;
1479 }
1480}
1481
1482/// modifies a
1483static void definiteGcdCancellation(number a, const coeffs cf,
1484 BOOLEAN simpleTestsHaveAlreadyBeenPerformed)
1485{
1486// ntTest(a); // !!!!
1487
1488 fraction f = (fraction)a;
1489
1490 if (IS0(a)) return;
1491 if (COM(f)==0) return;
1492 if (DENIS1(f) || NUMIS1(f)) { COM(f) = 0; ntTest(a); return; }
1493 if (!simpleTestsHaveAlreadyBeenPerformed)
1494 {
1495
1496 /* check whether NUM(f) = DEN(f), and - if so - replace 'a' by 1 */
1497 if (p_EqualPolys(NUM(f), DEN(f), ntRing))
1498 { /* numerator and denominator are both != 1 */
1499 p_Delete(&NUM(f), ntRing); NUM(f) = p_ISet(1, ntRing);
1500 p_Delete(&DEN(f), ntRing); DEN(f) = NULL;
1501 COM(f) = 0;
1502 ntTest(a);
1503 return;
1504 }
1505 }
1506 /*if (rField_is_Q(ntRing))
1507 {
1508 number c=n_Copy(pGetCoeff(NUM(f)),ntCoeffs);
1509 poly p=pNext(NUM(f));
1510 while((p!=NULL)&&(!n_IsOne(c,ntCoeffs)))
1511 {
1512 number cc=n_Gcd(c,pGetCoeff(p),ntCoeffs);
1513 n_Delete(&c,ntCoeffs);
1514 c=cc;
1515 pIter(p);
1516 };
1517 p=DEN(f);
1518 while((p!=NULL)&&(!n_IsOne(c,ntCoeffs)))
1519 {
1520 number cc=n_Gcd(c,pGetCoeff(p),ntCoeffs);
1521 n_Delete(&c,ntCoeffs);
1522 c=cc;
1523 pIter(p);
1524 };
1525 if(!n_IsOne(c,ntCoeffs))
1526 {
1527 p=NUM(f);
1528 do
1529 {
1530 number cc=n_Div(pGetCoeff(p),c,ntCoeffs);
1531 n_Normalize(cc,ntCoeffs);
1532 p_SetCoeff(p,cc,ntRing);
1533 pIter(p);
1534 } while(p!=NULL);
1535 p=DEN(f);
1536 do
1537 {
1538 number cc=n_Div(pGetCoeff(p),c,ntCoeffs);
1539 n_Normalize(cc,ntCoeffs);
1540 p_SetCoeff(p,cc,ntRing);
1541 pIter(p);
1542 } while(p!=NULL);
1543 n_Delete(&c,ntCoeffs);
1544 if(pNext(DEN(f))==NULL)
1545 {
1546 if (p_IsOne(DEN(f),ntRing))
1547 {
1548 p_LmDelete(&DEN(f),ntRing);
1549 COM(f)=0;
1550 return;
1551 }
1552 else
1553 {
1554 return;
1555 }
1556 }
1557 }
1558 }*/
1559
1560 /* here we assume: NUM(f), DEN(f) !=NULL, in Z_a reqp. Z/p_a */
1561 //StringSetS("");ntWriteLong(a,cf);
1562 poly pGcd = singclap_gcd_and_divide(NUM(f), DEN(f), ntRing);
1563 //PrintS("gcd= ");p_wrp(pGcd,ntRing);PrintLn();
1564 if (p_IsConstant(pGcd, ntRing)
1565 && n_IsOne(p_GetCoeff(pGcd, ntRing), ntCoeffs)
1566 )
1567 { /* gcd = 1; nothing to cancel;
1568 Suppose the given rational function field is over Q. Although the
1569 gcd is 1, we may have produced fractional coefficients in NUM(f),
1570 DEN(f), or both, due to previous arithmetics. The next call will
1571 remove those nested fractions, in case there are any. */
1573 {
1574 number d=p_GetCoeff (DEN(f),ntRing);
1575 BOOLEAN d_not_1=FALSE;
1576 if (!n_IsOne(d,ntCoeffs))
1577 {
1578 NUM (f) = p_Div_nn (NUM (f), d, ntRing);
1579 d_not_1=TRUE;
1580 }
1581 if (p_IsConstant (DEN (f), ntRing))
1582 {
1583 p_Delete(&DEN (f), ntRing);
1584 DEN (f) = NULL;
1585 }
1586 else if (d_not_1)
1587 {
1588 DEN (f) = p_Div_nn (DEN (f), d, ntRing);
1589 }
1591 }
1592 else
1593 { /* We divide both NUM(f) and DEN(f) by the gcd which is known
1594 to be != 1. */
1595 if (p_IsConstant(DEN(f), ntRing) &&
1597 {
1598 /* DEN(f) = 1 needs to be represented by NULL! */
1599 p_Delete(&DEN(f), ntRing);
1600 DEN(f) = NULL;
1601 }
1602 else
1603 {
1605 {
1606 NUM (f) = p_Div_nn (NUM (f), p_GetCoeff (DEN(f),ntRing), ntRing);
1607 if (p_IsConstant (DEN (f), ntRing))
1608 {
1609 p_Delete(&DEN (f), ntRing);
1610 DEN (f) = NULL;
1611 }
1612 else
1613 {
1614 p_Norm (DEN (f),ntRing);
1615 }
1616 }
1617 }
1618 }
1619 p_Delete(&pGcd, ntRing);
1620// StringAppendS(" -> ");ntWriteLong(a,cf);StringAppendS("\n");{ char* s = StringEndS(); Print("%s", s); omFree(s); }
1621 COM(f) = 0;
1622
1623 if( DEN(f) != NULL )
1624 {
1625 if( !n_GreaterZero(pGetCoeff(DEN(f)), ntCoeffs) )
1626 {
1627 NUM(f) = p_Neg(NUM(f), ntRing);
1628 DEN(f) = p_Neg(DEN(f), ntRing);
1629 if (p_IsConstant(DEN(f), ntRing) &&
1631 {
1632 /* DEN(f) = 1 needs to be represented by NULL! */
1633 p_Delete(&DEN(f), ntRing);
1634 DEN (f) = NULL;
1635 }
1636 }
1637 }
1638 ntTest(a); // !!!!
1639}
1640
1641static void ntWriteLong(number a, const coeffs cf)
1642{
1643 ntTest(a);
1644 if (IS0(a))
1645 StringAppendS("0");
1646 else
1647 {
1648 fraction f = (fraction)a;
1649 // stole logic from napWrite from kernel/longtrans.cc of legacy singular
1650 BOOLEAN omitBrackets = p_IsConstant(NUM(f), ntRing);
1651 if (!omitBrackets) StringAppendS("(");
1652 p_String0Long(NUM(f), ntRing, ntRing);
1653 if (!omitBrackets) StringAppendS(")");
1654 if (!DENIS1(f))
1655 {
1656 StringAppendS("/");
1657 omitBrackets = p_IsConstant(DEN(f), ntRing);
1658 if (!omitBrackets) StringAppendS("(");
1659 p_String0Long(DEN(f), ntRing, ntRing);
1660 if (!omitBrackets) StringAppendS(")");
1661 }
1662 }
1663 ntTest(a); // !!!!
1664}
1665
1666static void ntWriteShort(number a, const coeffs cf)
1667{
1668 ntTest(a);
1669 if (IS0(a))
1670 StringAppendS("0");
1671 else
1672 {
1673 fraction f = (fraction)a;
1674 // stole logic from napWrite from kernel/longtrans.cc of legacy singular
1675 BOOLEAN omitBrackets = p_IsConstant(NUM(f), ntRing);
1676 if (!omitBrackets) StringAppendS("(");
1677 p_String0Short(NUM(f), ntRing, ntRing);
1678 if (!omitBrackets) StringAppendS(")");
1679 if (!DENIS1(f))
1680 {
1681 StringAppendS("/");
1682 omitBrackets = p_IsConstant(DEN(f), ntRing);
1683 if (!omitBrackets) StringAppendS("(");
1684 p_String0Short(DEN(f), ntRing, ntRing);
1685 if (!omitBrackets) StringAppendS(")");
1686 }
1687 }
1688 ntTest(a);
1689}
1690
1691static const char * ntRead(const char *s, number *a, const coeffs cf)
1692{
1693 poly p;
1694 const char * result = p_Read(s, p, ntRing);
1695 if (p == NULL) *a = NULL;
1696 else *a = ntInit(p, cf);
1697 ntTest(*a);
1698 return result;
1699}
1700
1701static void ntNormalize (number &a, const coeffs cf)
1702{
1703 if ( /*(*/ a!=NULL /*)*/ )
1704 {
1705 //PrintS("num=");p_wrp(NUM(a),ntRing);
1706 //PrintS(" den=");p_wrp(DEN(a),ntRing);PrintLn();
1707 if (COM((fraction)a)>0) definiteGcdCancellation(a, cf, FALSE);
1708 if ((DEN((fraction)a)!=NULL)
1709 &&(!n_GreaterZero(pGetCoeff(DEN((fraction)a)),ntCoeffs)))
1710 {
1711 NUM((fraction)a)=p_Neg(NUM((fraction)a),ntRing);
1712 DEN((fraction)a)=p_Neg(DEN((fraction)a),ntRing);
1713 }
1714 }
1715 ntNormalizeDen((fraction)a,ntRing);
1716 ntTest(a); // !!!!
1717}
1718
1719static number ntExactDiv(number a, number b, const coeffs cf)
1720{
1721 number r=ntDiv(a,b,cf);
1722 ntNormalize(r,cf);
1723 return r;
1724}
1725
1726/* expects *param to be castable to TransExtInfo */
1727static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param)
1728{
1729 if (n_transExt != n) return FALSE;
1730 TransExtInfo *e = (TransExtInfo *)param;
1731 /* for rational function fields we expect the underlying
1732 polynomial rings to be IDENTICAL, i.e. the SAME OBJECT;
1733 this expectation is based on the assumption that we have properly
1734 registered cf and perform reference counting rather than creating
1735 multiple copies of the same coefficient field/domain/ring */
1736 if (ntRing == e->r)
1737 return TRUE;
1738
1739 // NOTE: Q(a)[x] && Q(a)[y] should better share the _same_ Q(a)...
1740 if( rEqual(ntRing, e->r, TRUE) )
1741 {
1742 rDelete(e->r);
1743 return TRUE;
1744 }
1745
1746 return FALSE;
1747}
1748
1749static number ntNormalizeHelper(number a, number b, const coeffs cf)
1750{
1751 ntTest(a);
1752 ntTest(b);
1753 fraction fb = (fraction)b;
1754 if ((b==NULL)||(DEN(fb)==NULL)) return ntCopy(a,cf);
1755 fraction fa = (fraction)a;
1756
1757 poly pGcd;
1758 if (nCoeff_is_Q(ntCoeffs))
1759 {
1760 poly pa = NUM(fa);
1761 poly pb = DEN(fb);
1763 {
1764 pGcd = p_Copy(pa,ntRing);
1765 p_SetCoeff (pGcd, n_Gcd (pGetCoeff(pGcd), pGetCoeff(pb), ntCoeffs), ntRing);
1766 }
1767 else
1768 {
1769 number contentpa, contentpb, tmp;
1770
1771 contentpb= n_Copy(p_GetCoeff(pb, ntRing),ntCoeffs);
1772 pIter(pb);
1773 while (pb != NULL)
1774 {
1775 tmp = n_SubringGcd(contentpb, p_GetCoeff(pb, ntRing) , ntCoeffs);
1776 n_Delete(&contentpb, ntCoeffs);
1777 contentpb = tmp;
1778 pIter(pb);
1779 }
1780
1781 contentpa= n_Copy(p_GetCoeff(pa, ntRing),ntCoeffs);
1782 pIter(pa);
1783 while (pa != NULL)
1784 {
1785 tmp = n_SubringGcd(contentpa, p_GetCoeff(pa, ntRing), ntCoeffs);
1786 n_Delete(&contentpa, ntCoeffs);
1787 contentpa = tmp;
1788 pIter(pa);
1789 }
1790
1791 tmp= n_SubringGcd (contentpb, contentpa, ntCoeffs);
1792 n_Delete(&contentpa, ntCoeffs);
1793 n_Delete(&contentpb, ntCoeffs);
1794 contentpa= tmp;
1795
1796 pGcd = gcd_over_Q(NUM(fa), DEN(fb), ntRing);
1797 pGcd= __p_Mult_nn (pGcd, contentpa, ntRing);
1798 n_Delete(&contentpa, ntCoeffs);
1799 }
1800 }
1801 else
1802 pGcd = singclap_gcd_r(NUM(fa), DEN(fb), ntRing);
1803
1804 /* Note that, over Q, singclap_gcd will remove the denominators in all
1805 rational coefficients of pa and pb, before starting to compute
1806 the gcd. Thus, we do not need to ensure that the coefficients of
1807 pa and pb live in Z; they may well be elements of Q\Z. */
1808
1809 if (p_IsConstant(pGcd, ntRing) &&
1811 { /* gcd = 1; return pa*pb*/
1812 p_Delete(&pGcd,ntRing);
1813 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1814 NUM(result) = pp_Mult_qq(NUM(fa),DEN(fb),ntRing);
1815
1816 ntTest((number)result); // !!!!
1817
1818 return (number)result;
1819 }
1820
1821 /* return pa*pb/gcd */
1822 poly newNum = singclap_pdivide(NUM(fa), pGcd, ntRing);
1823 p_Delete(&pGcd,ntRing);
1824 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1825 NUM(result) = p_Mult_q(p_Copy(DEN(fb),ntRing),newNum,ntRing);
1826 ntTest((number)result); // !!!!
1827 return (number)result;
1828}
1829
1830static number ntGcd(number a, number b, const coeffs cf)
1831{
1832 ntTest(a);
1833 ntTest(b);
1834 if (a==NULL) return ntCopy(b,cf);
1835 if (b==NULL) return ntCopy(a,cf);
1836 fraction fa = (fraction)a;
1837 fraction fb = (fraction)b;
1838
1839
1840 poly pGcd;
1841 if (nCoeff_is_Q(ntCoeffs))
1842 {
1843 poly pa = NUM(fa);
1844 poly pb = NUM(fb);
1846 {
1847 pGcd = p_Copy(pa,ntRing);
1849 }
1850 else
1851 {
1852 number contentpa, contentpb, tmp;
1853
1854 contentpb= n_Copy(p_GetCoeff(pb, ntRing),ntCoeffs);
1855 pIter(pb);
1856 while (pb != NULL)
1857 {
1858 tmp = n_SubringGcd(contentpb, p_GetCoeff(pb, ntRing) , ntCoeffs);
1859 n_Delete(&contentpb, ntCoeffs);
1860 contentpb = tmp;
1861 pIter(pb);
1862 }
1863
1864 contentpa= n_Copy(p_GetCoeff(pa, ntRing),ntCoeffs);
1865 pIter(pa);
1866 while (pa != NULL)
1867 {
1868 tmp = n_SubringGcd(contentpa, p_GetCoeff(pa, ntRing), ntCoeffs);
1869 n_Delete(&contentpa, ntCoeffs);
1870 contentpa = tmp;
1871 pIter(pa);
1872 }
1873
1874 tmp= n_SubringGcd (contentpb, contentpa, ntCoeffs);
1875 n_Delete(&contentpa, ntCoeffs);
1876 n_Delete(&contentpb, ntCoeffs);
1877 contentpa= tmp;
1878
1879 pGcd = gcd_over_Q(NUM(fa), NUM(fb), ntRing);
1880 pGcd= __p_Mult_nn (pGcd, contentpa, ntRing);
1881 n_Delete(&contentpa, ntCoeffs);
1882 }
1883 }
1884 else
1885 pGcd = singclap_gcd_r(NUM(fa), NUM(fb), ntRing);
1886 /* Note that, over Q, singclap_gcd will remove the denominators in all
1887 rational coefficients of pa and pb, before starting to compute
1888 the gcd. Thus, we do not need to ensure that the coefficients of
1889 pa and pb live in Z; they may well be elements of Q\Z. */
1890
1891 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
1892 NUM(result) = pGcd;
1893 ntTest((number)result); // !!!!
1894 return (number)result;
1895}
1896//number ntGcd_dummy(number a, number b, const coeffs cf)
1897//{
1898// extern char my_yylinebuf[80];
1899// Print("ntGcd in >>%s<<\n",my_yylinebuf);
1900// return ntGcd(a,b,cf);
1901//}
1902
1903static int ntSize(number a, const coeffs cf)
1904{
1905 ntTest(a);
1906 if (IS0(a)) return 0;
1907 fraction f = (fraction)a;
1908 poly p = NUM(f);
1909 unsigned long noOfTerms = 0;
1910 unsigned long numDegree = 0;
1911 if (p!=NULL)
1912 {
1913 numDegree = p_Totaldegree(p,ntRing);
1914 noOfTerms = pLength(p);
1915 }
1916 unsigned long denDegree = 0;
1917 if (!DENIS1(f))
1918 {
1919 denDegree = p_Totaldegree(DEN(f),ntRing);
1920 noOfTerms += pLength(DEN(f));
1921 }
1922 ntTest(a); // !!!!
1923 // avoid int overflow:
1924 unsigned long t= ((numDegree + denDegree)*(numDegree + denDegree) + 1) * noOfTerms; // must be >0
1925 if (t>INT_MAX) return INT_MAX;
1926 else return (int)t;
1927}
1928
1929/* assumes that src = Q or Z, dst = Q(t_1, ..., t_s) */
1930static number ntMap00(number a, const coeffs src, const coeffs dst)
1931{
1932 n_Test(a, src);
1933
1934 if (n_IsZero(a, src)) return NULL;
1935 assume(src->rep == dst->extRing->cf->rep);
1936 if ((SR_HDL(a) & SR_INT) || (a->s==3))
1937 {
1938 number res=ntInit(p_NSet(n_Copy(a, src), dst->extRing), dst);
1939 n_Test(res, dst);
1940 return res;
1941 }
1942 number nn=n_GetDenom(a,src);
1943 number zz=n_GetNumerator(a,src);
1944 number res=ntInit(p_NSet(zz,dst->extRing), dst);
1945 fraction ff=(fraction)res;
1946 if (n_IsOne(nn,src)) DEN(ff)=NULL;
1947 else DEN(ff)=p_NSet(nn,dst->extRing);
1948
1949 n_Test((number)ff,dst);
1950 //check_N((number)ff,dst);
1951 return (number)ff;
1952}
1953
1954static number ntMapZ0(number a, const coeffs src, const coeffs dst)
1955{
1956 n_Test(a, src);
1957 if (n_IsZero(a, src)) return NULL;
1958 nMapFunc nMap=n_SetMap(src,dst->extRing->cf);
1959 poly p=p_NSet(nMap(a, src,dst->extRing->cf), dst->extRing);
1960 if (n_IsZero(pGetCoeff(p),dst->extRing->cf))
1961 p_Delete(&p,dst->extRing);
1962 number res=ntInit(p, dst);
1963 n_Test(res,dst);
1964 return res;
1965}
1966
1967/* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */
1968static number ntMapP0(number a, const coeffs src, const coeffs dst)
1969{
1970 n_Test(a, src);
1971 if (n_IsZero(a, src)) return NULL;
1972 /* mapping via intermediate int: */
1973 int n = n_Int(a, src);
1974 number q = n_Init(n, dst->extRing->cf);
1975 if (n_IsZero(q, dst->extRing->cf))
1976 {
1977 n_Delete(&q, dst->extRing->cf);
1978 return NULL;
1979 }
1980 return ntInit(p_NSet(q, dst->extRing), dst);
1981}
1982
1983 /* assumes that either src = K(t_1, ..., t_s), dst = K(t_1, ..., t_s) */
1984static number ntCopyMap(number a, const coeffs cf, const coeffs dst)
1985{
1986 ntTest(a);
1987 if (IS0(a)) return NULL;
1988
1989 const ring rSrc = cf->extRing;
1990 const ring rDst = dst->extRing;
1991
1992 if( rSrc == rDst )
1993 return ntCopy(a, dst); // USUALLY WRONG!
1994
1995 fraction f = (fraction)a;
1996 poly g = prCopyR(NUM(f), rSrc, rDst);
1997
1998 poly h = NULL;
1999
2000 if (!DENIS1(f))
2001 h = prCopyR(DEN(f), rSrc, rDst);
2002
2003 fraction result = (fraction)omAllocBin(fractionObjectBin);
2004
2005 NUM(result) = g;
2006 DEN(result) = h;
2007 COM(result) = COM(f);
2008 //check_N((number)result,dst);
2009 n_Test((number)result, dst);
2010 return (number)result;
2011}
2012
2013static number ntGenMap(number a, const coeffs cf, const coeffs dst)
2014{
2015 ntTest(a);
2016 if (IS0(a)) return NULL;
2017
2018 const ring rSrc = cf->extRing;
2019 const ring rDst = dst->extRing;
2020
2021 const nMapFunc nMap=n_SetMap(rSrc->cf,rDst->cf);
2022 fraction f = (fraction)a;
2023 poly g = prMapR(NUM(f), nMap, rSrc, rDst);
2024 /* g may contain summands with coeff 0 */
2025 poly hh=g;
2026 poly prev=NULL;
2027 while(hh!=NULL)
2028 {
2029 if (n_IsZero(pGetCoeff(hh),rDst->cf))
2030 {
2031 if (prev==NULL)
2032 {
2033 g=p_LmFreeAndNext(g,rDst);
2034 hh=g;
2035 }
2036 else
2037 {
2038 prev->next=p_LmFreeAndNext(prev->next,rDst);
2039 hh=prev->next;
2040 }
2041 }
2042 else
2043 {
2044 prev=hh;
2045 pIter(hh);
2046 }
2047 }
2048 if (g==NULL) return NULL;
2049
2050 poly h = NULL;
2051
2052 if (!DENIS1(f))
2053 {
2054 h = prMapR(DEN(f), nMap, rSrc, rDst);
2055 /* h may contain summands with coeff 0 */
2056 hh=h;
2057 prev=NULL;
2058 while(hh!=NULL)
2059 {
2060 if (n_IsZero(pGetCoeff(hh),rDst->cf))
2061 {
2062 if (prev==NULL)
2063 {
2064 h=p_LmFreeAndNext(h,rDst);
2065 hh=h;
2066 }
2067 else
2068 {
2069 prev->next=p_LmFreeAndNext(prev->next,rDst);
2070 hh=prev->next;
2071 }
2072 }
2073 else
2074 {
2075 prev=hh;
2076 pIter(hh);
2077 }
2078 }
2079 if (h==NULL) WerrorS("mapping to */0");
2080 }
2081
2082 fraction result = (fraction)omAllocBin(fractionObjectBin);
2083
2084 NUM(result) = g;
2085 DEN(result) = h;
2086 COM(result) = COM(f);
2087 //check_N((number)result,dst);
2088 n_Test((number)result, dst);
2089 return (number)result;
2090}
2091
2092static number ntCopyAlg(number a, const coeffs cf, const coeffs dst)
2093{
2094 n_Test(a, cf) ;
2095 if (n_IsZero(a, cf)) return NULL;
2096 return ntInit(prCopyR((poly)a, cf->extRing, dst->extRing),dst);
2097}
2098
2099static number ntGenAlg(number a, const coeffs cf, const coeffs dst)
2100{
2101 n_Test(a, cf) ;
2102 if (n_IsZero(a, cf)) return NULL;
2103
2104 const nMapFunc nMap=n_SetMap(cf->extRing->cf,dst->extRing->cf);
2105 return ntInit(prMapR((poly)a, nMap, cf->extRing, dst->extRing),dst);
2106}
2107
2108/* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */
2109static number ntMap0P(number a, const coeffs src, const coeffs dst)
2110{
2111 n_Test(a, src) ;
2112 if (n_IsZero(a, src)) return NULL;
2113 // int p = rChar(dst->extRing);
2114
2115 number q = nlModP(a, src, dst->extRing->cf); // FIXME? TODO? // extern number nlModP(number q, const coeffs Q, const coeffs Zp); // Map q \in QQ \to Zp
2116
2117 if (n_IsZero(q, dst->extRing->cf))
2118 {
2119 n_Delete(&q, dst->extRing->cf);
2120 return NULL;
2121 }
2122
2123 poly g = p_NSet(q, dst->extRing);
2124
2125 fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
2126 NUM(f) = g; // DEN(f) = NULL; COM(f) = 0;
2127 n_Test((number)f, dst);
2128 //check_N((number)f,dst);
2129 return (number)f;
2130}
2131
2132/* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */
2133static number ntMapPP(number a, const coeffs src, const coeffs dst)
2134{
2135 n_Test(a, src) ;
2136 if (n_IsZero(a, src)) return NULL;
2137 assume(src == dst->extRing->cf);
2138 poly p = p_One(dst->extRing);
2139 p_SetCoeff(p, n_Copy(a, src), dst->extRing);
2140 fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
2141 NUM(f) = p; // DEN(f) = NULL; COM(f) = 0;
2142 n_Test((number)f, dst);
2143 //check_N((number)f,dst);
2144 return (number)f;
2145}
2146
2147/* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */
2148static number ntMapUP(number a, const coeffs src, const coeffs dst)
2149{
2150 n_Test(a, src) ;
2151 if (n_IsZero(a, src)) return NULL;
2152 /* mapping via intermediate int: */
2153 int n = n_Int(a, src);
2154 number q = n_Init(n, dst->extRing->cf);
2155 poly p;
2156 if (n_IsZero(q, dst->extRing->cf))
2157 {
2158 n_Delete(&q, dst->extRing->cf);
2159 return NULL;
2160 }
2161 p = p_One(dst->extRing);
2162 p_SetCoeff(p, q, dst->extRing);
2163 fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
2164 NUM(f) = p; // DEN(f) = NULL; COM(f) = 0;
2165 n_Test((number)f, dst);
2166 //check_N((number)f,dst);
2167 return (number)f;
2168}
2169
2170nMapFunc ntSetMap(const coeffs src, const coeffs dst)
2171{
2172 /* dst is expected to be a rational function field */
2174
2175 int h = 0; /* the height of the extension tower given by dst */
2176 coeffs bDst = nCoeff_bottom(dst, h); /* the bottom field in the tower dst */
2177 coeffs bSrc = nCoeff_bottom(src, h); /* the bottom field in the tower src */
2178
2179 /* for the time being, we only provide maps if h = 1 and if b is Q or
2180 some field Z/pZ: */
2181 if (h==0)
2182 {
2183 if (((src->rep==n_rep_gap_rat) || (src->rep==n_rep_gap_gmp))
2184 && (nCoeff_is_Q(dst->extRing->cf) || nCoeff_is_Z(dst->extRing->cf)))
2185 return ntMap00; /// Q or Z --> Q(T)
2186 if (src->rep==n_rep_gmp)
2187 return ntMapZ0; /// Z --> K(T)
2188 if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst))
2189 return ntMapP0; /// Z/p --> Q(T)
2190 if (nCoeff_is_Q_or_BI(src) && nCoeff_is_Zp(bDst))
2191 return ntMap0P; /// Q --> Z/p(T)
2192 if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst))
2193 {
2194 if (src->ch == dst->ch) return ntMapPP; /// Z/p --> Z/p(T)
2195 else return ntMapUP; /// Z/u --> Z/p(T)
2196 }
2197 if (nCoeff_is_Zn(src) && nCoeff_is_Zn(bDst))
2198 {
2199 if (mpz_cmp(src->modNumber,bDst->modNumber)==0) return ntMapPP; /// Z/p --> Z/p(T)
2200 }
2201 }
2202 if (h != 1) return NULL;
2203 //if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL;
2204
2205 /* Let T denote the sequence of transcendental extension variables, i.e.,
2206 K[t_1, ..., t_s] =: K[T];
2207 Let moreover, for any such sequence T, T' denote any subsequence of T
2208 of the form t_1, ..., t_w with w <= s. */
2209
2210 if (rVar(src->extRing) > rVar(dst->extRing))
2211 return NULL;
2212
2213 for (int i = 0; i < rVar(src->extRing); i++)
2214 if (strcmp(rRingVar(i, src->extRing), rRingVar(i, dst->extRing)) != 0)
2215 return NULL;
2216
2217 if (src->type==n_transExt)
2218 {
2219 if (src->extRing->cf==dst->extRing->cf)
2220 return ntCopyMap; /// K(T') --> K(T)
2221 else
2222 return ntGenMap; /// K(T') --> K'(T)
2223 }
2224 else
2225 {
2226 if (src->extRing->cf==dst->extRing->cf)
2227 return ntCopyAlg; /// K(T') --> K(T)
2228 else
2229 return ntGenAlg; /// K(T') --> K'(T)
2230 }
2231
2232 return NULL; /// default
2233}
2234#if 0
2235nMapFunc ntSetMap_T(const coeffs src, const coeffs dst)
2236{
2237 nMapFunc n=ntSetMap(src,dst);
2238 if (n==ntCopyAlg) printf("n=ntCopyAlg\n");
2239 else if (n==ntCopyMap) printf("n=ntCopyMap\n");
2240 else if (n==ntMapUP) printf("n=ntMapUP\n");
2241 else if (n==ntMap0P) printf("n=ntMap0P\n");
2242 else if (n==ntMapP0) printf("n=ntMapP0\n");
2243 else if (n==ntMap00) printf("n=ntMap00\n");
2244 else if (n==NULL) printf("n=NULL\n");
2245 else printf("n=?\n");
2246 return n;
2247}
2248#endif
2249
2251{
2252 rDecRefCnt(cf->extRing);
2253 if (cf->extRing->ref < 0)
2254 rDelete(cf->extRing);
2255}
2256static number ntConvFactoryNSingN( const CanonicalForm n, const coeffs cf)
2257{
2258 if (n.isZero()) return NULL;
2259 poly p=convFactoryPSingP(n,ntRing);
2261 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
2262 NUM(result) = p;
2263 //DEN(result) = NULL; // done by omAlloc0Bin
2264 //COM(result) = 0; // done by omAlloc0Bin
2265 ntTest((number)result);
2266 return (number)result;
2267}
2268static CanonicalForm ntConvSingNFactoryN( number n, BOOLEAN /*setChar*/, const coeffs cf )
2269{
2270 ntTest(n);
2271 if (IS0(n)) return CanonicalForm(0);
2272
2273 fraction f = (fraction)n;
2274 return convSingPFactoryP(NUM(f),ntRing);
2275}
2276
2277static int ntParDeg(number a, const coeffs cf)
2278{
2279 ntTest(a);
2280 if (IS0(a)) return -1;
2281 fraction fa = (fraction)a;
2282 return cf->extRing->pFDeg(NUM(fa),cf->extRing);
2283}
2284
2285/// return the specified parameter as a number in the given trans.ext.
2286static number ntParameter(const int iParameter, const coeffs cf)
2287{
2289
2290 const ring R = cf->extRing;
2291 assume( R != NULL );
2292 assume( 0 < iParameter && iParameter <= rVar(R) );
2293
2294 poly p = p_One(R); p_SetExp(p, iParameter, 1, R); p_Setm(p, R);
2295 p_Test(p,R);
2296
2297 fraction f = (fraction)omAlloc0Bin(fractionObjectBin);
2298 NUM(f) = p;
2299 //DEN(f) = NULL;
2300 //COM(f) = 0;
2301
2302 ntTest((number)f);
2303
2304 return (number)f;
2305}
2306
2307/// if m == var(i)/1 => return i,
2308int ntIsParam(number m, const coeffs cf)
2309{
2310 ntTest(m);
2312
2313 const ring R = cf->extRing;
2314 assume( R != NULL );
2315
2316 fraction f = (fraction)m;
2317
2318 if( DEN(f) != NULL )
2319 return 0;
2320
2321 return p_Var( NUM(f), R );
2322}
2323
2325{
2326 static inline poly convert(const number& n)
2327 {
2328 // suitable for trans. ext. numbers that are fractions of polys
2329 return NUM((fraction)n); // return the numerator
2330 }
2331};
2332
2333
2334static void ntClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
2335{
2336 assume(cf != NULL);
2338 // all coeffs are given by fractions of polynomails over integers!!!
2339 // without denominators!!!
2340
2341 const ring R = cf->extRing;
2342 assume(R != NULL);
2343 const coeffs Q = R->cf;
2344 assume(Q != NULL);
2346
2347
2348 numberCollectionEnumerator.Reset();
2349
2350 if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2351 {
2352 c = ntInit(1, cf);
2353 return;
2354 }
2355
2356 // all coeffs are given by integers after returning from this routine
2357
2358 // part 1, collect product of all denominators /gcds
2359 poly cand = NULL;
2360
2361 do
2362 {
2363 number &n = numberCollectionEnumerator.Current();
2364
2365 ntNormalize(n, cf);
2366
2367 fraction f = (fraction)n;
2368
2369 assume( f != NULL );
2370
2371 const poly den = DEN(f);
2372
2373 assume( den == NULL ); // ?? / 1 ?
2374
2375 const poly num = NUM(f);
2376
2377 if( cand == NULL )
2378 cand = p_Copy(num, R);
2379 else
2380 {
2381 poly tmp = singclap_gcd_r(cand, num, R); // gcd(cand, num)
2382 p_Delete(&cand,R);
2383 cand=tmp;
2384 }
2385
2386 if( p_IsConstant(cand, R) )
2387 break;
2388 }
2389 while( numberCollectionEnumerator.MoveNext() ) ;
2390
2391
2392 // part2: all coeffs = all coeffs * cand
2393 if( cand != NULL )
2394 {
2395 if( !p_IsConstant(cand, R) )
2396 {
2397 c = ntInit(cand, cf);
2398 numberCollectionEnumerator.Reset();
2399 while (numberCollectionEnumerator.MoveNext() )
2400 {
2401 number &n = numberCollectionEnumerator.Current();
2402 const number t = ntDiv(n, c, cf); // TODO: rewrite!?
2403 ntDelete(&n, cf);
2404 n = t;
2405 }
2406 } // else NUM (result) = p_One(R);
2407 else { p_Delete(&cand, R); cand = NULL; }
2408 }
2409
2410 // Quick and dirty fix for constant content clearing: consider numerators???
2411 CRecursivePolyCoeffsEnumerator<NTNumConverter> itr(numberCollectionEnumerator); // recursively treat the NUM(numbers) as polys!
2412 number cc;
2413
2414 n_ClearContent(itr, cc, Q);
2415 number g = ntInit(p_NSet(cc, R), cf);
2416
2417 if( cand != NULL )
2418 {
2419 number gg = ntMult(g, c, cf);
2420 ntDelete(&g, cf);
2421 ntDelete(&c, cf); c = gg;
2422 } else
2423 c = g;
2424 ntTest(c);
2425}
2426
2427static void ntClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
2428{
2429 assume(cf != NULL);
2430 assume(getCoeffType(cf) == n_transExt); // both over Q(a) and Zp(a)!
2431 // all coeffs are given by fractions of polynomails over integers!!!
2432
2433 numberCollectionEnumerator.Reset();
2434
2435 if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
2436 {
2437 c = ntInit(1, cf);
2438 return;
2439 }
2440
2441 // all coeffs are given by integers after returning from this routine
2442
2443 // part 1, collect product of all denominators /gcds
2444 poly cand = NULL;
2445
2446 const ring R = cf->extRing;
2447 assume(R != NULL);
2448
2449 const coeffs Q = R->cf;
2450 assume(Q != NULL);
2451// assume(nCoeff_is_Q(Q));
2452
2453 do
2454 {
2455 number &n = numberCollectionEnumerator.Current();
2456
2457 ntNormalize(n, cf);
2458
2459 fraction f = (fraction)ntGetDenom (n, cf);
2460
2461 assume( f != NULL );
2462
2463 const poly den = NUM(f);
2464
2465 if( den == NULL ) // ?? / 1 ?
2466 continue;
2467
2468 if( cand == NULL )
2469 cand = p_Copy(den, R);
2470 else
2471 {
2472 // cand === LCM( cand, den )!!!!
2473 // NOTE: maybe it's better to make the product and clearcontent afterwards!?
2474 // TODO: move the following to factory?
2475 poly gcd = singclap_gcd_r(cand, den, R); // gcd(cand, den) is monic no mater leading coeffs! :((((
2476 if (nCoeff_is_Q (Q))
2477 {
2478 number LcGcd= n_SubringGcd (p_GetCoeff (cand, R), p_GetCoeff(den, R), Q);
2479 gcd = __p_Mult_nn(gcd, LcGcd, R);
2480 n_Delete(&LcGcd,Q);
2481 }
2482// assume( n_IsOne(pGetCoeff(gcd), Q) ); // TODO: this may be wrong...
2483 cand = p_Mult_q(cand, p_Copy(den, R), R); // cand *= den
2484 const poly t = singclap_pdivide( cand, gcd, R ); // cand' * den / gcd(cand', den)
2485 p_Delete(&cand, R);
2486 p_Delete(&gcd, R);
2487 cand = t;
2488 }
2489 }
2490 while( numberCollectionEnumerator.MoveNext() );
2491
2492 if( cand == NULL )
2493 {
2494 c = ntInit(1, cf);
2495 return;
2496 }
2497
2498 c = ntInit(cand, cf);
2499
2500 numberCollectionEnumerator.Reset();
2501
2502 number d = NULL;
2503
2504 while (numberCollectionEnumerator.MoveNext() )
2505 {
2506 number &n = numberCollectionEnumerator.Current();
2507 number t = ntMult(n, c, cf); // TODO: rewrite!?
2508 ntDelete(&n, cf);
2509
2510 ntNormalize(t, cf); // TODO: needed?
2511 n = t;
2512
2513 fraction f = (fraction)t;
2514 assume( f != NULL );
2515
2516 const poly den = DEN(f);
2517
2518 if( den != NULL ) // ?? / ?? ?
2519 {
2520 assume( p_IsConstant(den, R) );
2521 assume( pNext(den) == NULL );
2522
2523 if( d == NULL )
2524 d = n_Copy(pGetCoeff(den), Q);
2525 else
2526 {
2527 number g = n_NormalizeHelper(d, pGetCoeff(den), Q);
2528 n_Delete(&d, Q); d = g;
2529 }
2530 }
2531 }
2532
2533 if( d != NULL )
2534 {
2535 numberCollectionEnumerator.Reset();
2536 while (numberCollectionEnumerator.MoveNext() )
2537 {
2538 number &n = numberCollectionEnumerator.Current();
2539 fraction f = (fraction)n;
2540
2541 assume( f != NULL );
2542
2543 const poly den = DEN(f);
2544
2545 if( den == NULL ) // ?? / 1 ?
2546 NUM(f) = __p_Mult_nn(NUM(f), d, R);
2547 else
2548 {
2549 assume( p_IsConstant(den, R) );
2550 assume( pNext(den) == NULL );
2551
2552 number ddd = n_Div(d, pGetCoeff(den), Q); // but be an integer now!!!
2553 NUM(f) = __p_Mult_nn(NUM(f), ddd, R);
2554 n_Delete(&ddd, Q);
2555
2556 p_Delete(&DEN(f), R);
2557 DEN(f) = NULL; // TODO: check if this is needed!?
2558 }
2559
2560 assume( DEN(f) == NULL );
2561 }
2562
2563 NUM((fraction)c) = __p_Mult_nn(NUM((fraction)c), d, R);
2564 n_Delete(&d, Q);
2565 }
2566
2567
2568 ntTest(c);
2569}
2570
2571static number ntChineseRemainder(number *x, number *q,int rl, BOOLEAN /*sym*/,CFArray &inv_cache,const coeffs cf)
2572{
2573 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
2574
2575 poly *P=(poly*)omAlloc(rl*sizeof(poly*));
2576 number *X=(number *)omAlloc(rl*sizeof(number));
2577
2578 int i;
2579
2580 for(i=0;i<rl;i++) P[i]=p_Copy(NUM((fraction)(x[i])),cf->extRing);
2581 NUM(result)=p_ChineseRemainder(P,X,q,rl,inv_cache,cf->extRing);
2582
2583 for(i=0;i<rl;i++)
2584 {
2585 P[i]=p_Copy(DEN((fraction)(x[i])),cf->extRing);
2586 if (P[i]==NULL) P[i]=p_One(cf->extRing);
2587 }
2588 DEN(result)=p_ChineseRemainder(P,X,q,rl,inv_cache,cf->extRing);
2589
2590 omFreeSize(X,rl*sizeof(number));
2591 omFreeSize(P,rl*sizeof(poly*));
2592 if (p_IsConstant(DEN(result), ntRing)
2593 && n_IsOne(pGetCoeff(DEN(result)), ntCoeffs))
2594 {
2595 p_Delete(&DEN(result),ntRing);
2596 }
2597 ntTest((number)result);
2598 return ((number)result);
2599}
2600
2601static number ntFarey(number p, number n, const coeffs cf)
2602{
2603 // n is really a bigint
2604 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
2605 NUM(result)=p_Farey(p_Copy(NUM((fraction)p),cf->extRing),n,cf->extRing);
2606 DEN(result)=p_Farey(p_Copy(DEN((fraction)p),cf->extRing),n,cf->extRing);
2607 ntTest((number)result);
2608 return ((number)result);
2609}
2610
2611static number ntInitMPZ(mpz_t m, const coeffs r)
2612{
2613 fraction result = (fraction)omAlloc0Bin(fractionObjectBin);
2614 number n=n_InitMPZ(m,r->extRing->cf);
2615 NUM(result)=p_NSet(n,r->extRing);
2616 return ((number)result);
2617}
2618
2619static void ntMPZ(mpz_t m, number &n, const coeffs r)
2620{
2621 mpz_init(m);
2622 if (n!=NULL)
2623 {
2624 fraction nn=(fraction)n;
2625 if (DENIS1(nn))
2626 {
2627 if (p_IsConstant(NUM(nn),r->extRing))
2628 {
2629 n_MPZ(m,pGetCoeff(NUM(nn)),r->extRing->cf);
2630 return;
2631 }
2632 }
2633 }
2634}
2635
2636BOOLEAN ntInitChar(coeffs cf, void * infoStruct)
2637{
2638
2639 assume( infoStruct != NULL );
2640
2641 TransExtInfo *e = (TransExtInfo *)infoStruct;
2642
2643 assume( e->r != NULL); // extRing;
2644 assume( e->r->cf != NULL); // extRing->cf;
2645 assume( e->r->qideal == NULL );
2646
2647 assume( cf != NULL );
2648 assume(getCoeffType(cf) == n_transExt); // coeff type;
2649
2650 ring R = e->r;
2651 assume(R != NULL);
2652
2653 rIncRefCnt(R); // increase the ref.counter for the ground poly. ring!
2654
2655 cf->extRing = R;
2656 /* propagate characteristic up so that it becomes
2657 directly accessible in cf: */
2658 cf->ch = R->cf->ch;
2659
2660 cf->is_field=TRUE;
2661 cf->is_domain=TRUE;
2662 cf->rep=n_rep_rat_fct;
2663
2664 cf->factoryVarOffset = R->cf->factoryVarOffset + rVar(R);
2665
2666 cf->cfCoeffName = naCoeffName; // FIXME? TODO? // extern char* naCoeffString(const coeffs r);
2667
2668 cf->cfGreaterZero = ntGreaterZero;
2669 cf->cfGreater = ntGreater;
2670 cf->cfEqual = ntEqual;
2671 cf->cfIsZero = ntIsZero;
2672 cf->cfIsOne = ntIsOne;
2673 cf->cfIsMOne = ntIsMOne;
2674 cf->cfInit = ntInit;
2675 cf->cfFarey = ntFarey;
2676 cf->cfChineseRemainder = ntChineseRemainder;
2677 cf->cfInt = ntInt;
2678 cf->cfAdd = ntAdd;
2679 cf->cfInpAdd = ntInpAdd;
2680 cf->cfInpNeg = ntNeg;
2681 cf->cfSub = ntSub;
2682 cf->cfMult = ntMult;
2683 cf->cfInpMult = ntInpMult;
2684 cf->cfDiv = ntDiv;
2685 cf->cfExactDiv = ntExactDiv;
2686 cf->cfPower = ntPower;
2687 cf->cfCopy = ntCopy;
2688 cf->cfWriteLong = ntWriteLong;
2689 cf->cfRead = ntRead;
2690 cf->cfNormalize = ntNormalize;
2691 cf->cfDelete = ntDelete;
2692 cf->cfSetMap = ntSetMap;
2693 cf->cfGetDenom = ntGetDenom;
2694 cf->cfGetNumerator = ntGetNumerator;
2695 //cf->cfRePart = ntCopy;
2696 //cf->cfImPart = ntImPart;
2697 cf->cfCoeffWrite = ntCoeffWrite;
2698#ifdef LDEBUG
2699 cf->cfDBTest = ntDBTest;
2700#endif
2701 //cf->cfGcd = ntGcd_dummy;
2702 cf->cfSubringGcd = ntGcd;
2703 cf->cfNormalizeHelper = ntNormalizeHelper;
2704 cf->cfSize = ntSize;
2705 cf->nCoeffIsEqual = ntCoeffIsEqual;
2706 cf->cfInvers = ntInvers;
2707 cf->cfKillChar = ntKillChar;
2708 cf->cfInitMPZ = ntInitMPZ;
2709 cf->cfMPZ = ntMPZ;
2710
2711 if( rCanShortOut(ntRing) )
2712 cf->cfWriteShort = ntWriteShort;
2713 else
2714 cf->cfWriteShort = ntWriteLong;
2715
2716 cf->convFactoryNSingN =ntConvFactoryNSingN;
2717 cf->convSingNFactoryN =ntConvSingNFactoryN;
2718 cf->cfParDeg = ntParDeg;
2719
2720 cf->iNumberOfParameters = rVar(R);
2721 cf->pParameterNames = (const char**)R->names;
2722 cf->cfParameter = ntParameter;
2723 cf->has_simple_Inverse= FALSE;
2724 /* cf->has_simple_Alloc= FALSE; */
2725
2726
2727 if( nCoeff_is_Q(R->cf) )
2728 cf->cfClearContent = ntClearContent;
2729
2730 cf->cfClearDenominators = ntClearDenominators;
2731
2732 return FALSE;
2733}
2734
2736template class IEnumerator<snumber*>;
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
Concrete implementation of enumerators over polynomials.
char * naCoeffName(const coeffs r)
Definition: algext.cc:1345
All the auxiliary stuff.
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
void * ADDRESS
Definition: auxiliary.h:119
CanonicalForm num(const CanonicalForm &f)
CanonicalForm den(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
Definition: cfModGcd.cc:69
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
g
Definition: cfModGcd.cc:4090
CanonicalForm cf
Definition: cfModGcd.cc:4083
CanonicalForm b
Definition: cfModGcd.cc:4103
FILE * f
Definition: checklibs.c:9
CanonicalForm convSingPFactoryP(poly p, const ring r)
Definition: clapconv.cc:136
poly convFactoryPSingP(const CanonicalForm &f, const ring r)
Definition: clapconv.cc:40
poly singclap_pdivide(poly f, poly g, const ring r)
Definition: clapsing.cc:624
poly singclap_gcd_and_divide(poly &f, poly &g, const ring r)
clears denominators of f and g, divides by gcd(f,g)
Definition: clapsing.cc:170
poly singclap_gcd_r(poly f, poly g, const ring r)
Definition: clapsing.cc:68
This is a polynomial enumerator for simple iteration over coefficients of polynomials.
go into polynomials over an alg. extension recursively
factory's main class
Definition: canonicalform.h:86
CF_NO_INLINE bool isZero() const
virtual reference Current()=0
Gets the current element in the collection (read and write).
virtual void Reset()=0
Sets the enumerator to its initial position: -1, which is before the first element in the collection.
virtual bool MoveNext()=0
Advances the enumerator to the next element of the collection. returns true if the enumerator was suc...
Templated enumerator interface for simple iteration over a generic collection of T's.
Definition: Enumerator.h:125
Coefficient rings, fields and other domains suitable for Singular polynomials.
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of 'a' and 'b', i.e., a*b
Definition: coeffs.h:633
static FORCE_INLINE long n_Int(number &n, const coeffs r)
conversion of n to an int; 0 if not possible in Z/pZ: the representing int lying in (-p/2 ....
Definition: coeffs.h:544
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:448
static FORCE_INLINE number n_NormalizeHelper(number a, number b, const coeffs r)
assume that r is a quotient field (otherwise, return 1) for arguments (a1/a2,b1/b2) return (lcm(a1,...
Definition: coeffs.h:692
static FORCE_INLINE number n_GetDenom(number &n, const coeffs r)
return the denominator of n (if elements of r are by nature not fractional, result is 1)
Definition: coeffs.h:600
static FORCE_INLINE void n_CoeffWrite(const coeffs r, BOOLEAN details=TRUE)
output the coeff description
Definition: coeffs.h:716
static FORCE_INLINE BOOLEAN nCoeff_is_Z(const coeffs r)
Definition: coeffs.h:813
static FORCE_INLINE BOOLEAN nCoeff_is_Extension(const coeffs r)
Definition: coeffs.h:843
#define n_Test(a, r)
BOOLEAN n_Test(number a, const coeffs r)
Definition: coeffs.h:709
n_coeffType
Definition: coeffs.h:27
@ n_Q
rational (GMP) numbers
Definition: coeffs.h:30
@ n_Zp
\F{p < 2^31}
Definition: coeffs.h:29
@ n_transExt
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition: coeffs.h:38
static FORCE_INLINE number n_Gcd(number a, number b, const coeffs r)
in Z: return the gcd of 'a' and 'b' in Z/nZ, Z/2^kZ: computed as in the case Z in Z/pZ,...
Definition: coeffs.h:661
static FORCE_INLINE number n_Invers(number a, const coeffs r)
return the multiplicative inverse of 'a'; raise an error if 'a' is not invertible
Definition: coeffs.h:561
static FORCE_INLINE BOOLEAN n_GreaterZero(number n, const coeffs r)
ordered fields: TRUE iff 'n' is positive; in Z/pZ: TRUE iff 0 < m <= roundedBelow(p/2),...
Definition: coeffs.h:491
static FORCE_INLINE void n_MPZ(mpz_t result, number &n, const coeffs r)
conversion of n to a GMP integer; 0 if not possible
Definition: coeffs.h:548
static FORCE_INLINE BOOLEAN n_IsMOne(number n, const coeffs r)
TRUE iff 'n' represents the additive inverse of the one element, i.e. -1.
Definition: coeffs.h:469
static FORCE_INLINE nMapFunc n_SetMap(const coeffs src, const coeffs dst)
set the mapping function pointers for translating numbers from src to dst
Definition: coeffs.h:697
static FORCE_INLINE number n_InpNeg(number n, const coeffs r)
in-place negation of n MUST BE USED: n = n_InpNeg(n) (no copy is returned)
Definition: coeffs.h:554
static FORCE_INLINE BOOLEAN nCoeff_has_simple_inverse(const coeffs r)
TRUE, if the computation of the inverse is fast, i.e. prefer leading coeff. 1 over content.
Definition: coeffs.h:899
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition: coeffs.h:612
static FORCE_INLINE BOOLEAN nCoeff_is_Q(const coeffs r)
Definition: coeffs.h:803
static FORCE_INLINE BOOLEAN n_Greater(number a, number b, const coeffs r)
ordered fields: TRUE iff 'a' is larger than 'b'; in Z/pZ: TRUE iff la > lb, where la and lb are the l...
Definition: coeffs.h:508
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
Definition: coeffs.h:461
static FORCE_INLINE void n_ClearDenominators(ICoeffsEnumerator &numberCollectionEnumerator, number &d, const coeffs r)
(inplace) Clears denominators on a collection of numbers number d is the LCM of all the coefficient d...
Definition: coeffs.h:932
static FORCE_INLINE BOOLEAN nCoeff_is_Q_or_BI(const coeffs r)
Definition: coeffs.h:826
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition: coeffs.h:422
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:452
static FORCE_INLINE BOOLEAN nCoeff_is_Zn(const coeffs r)
Definition: coeffs.h:823
static FORCE_INLINE number n_InitMPZ(mpz_t n, const coeffs r)
conversion of a GMP integer to number
Definition: coeffs.h:539
static FORCE_INLINE BOOLEAN nCoeff_is_Zp(const coeffs r)
Definition: coeffs.h:797
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:535
static FORCE_INLINE void n_ClearContent(ICoeffsEnumerator &numberCollectionEnumerator, number &c, const coeffs r)
Computes the content and (inplace) divides it out on a collection of numbers number c is the content ...
Definition: coeffs.h:925
static FORCE_INLINE number n_GetNumerator(number &n, const coeffs r)
return the numerator of n (if elements of r are by nature not fractional, result is n)
Definition: coeffs.h:605
@ n_rep_gap_rat
(number), see longrat.h
Definition: coeffs.h:111
@ n_rep_gap_gmp
(), see rinteger.h, new impl.
Definition: coeffs.h:112
@ n_rep_gmp
(mpz_ptr), see rmodulon,h
Definition: coeffs.h:115
@ n_rep_rat_fct
(fraction), see transext.h
Definition: coeffs.h:114
static FORCE_INLINE number n_SubringGcd(number a, number b, const coeffs r)
Definition: coeffs.h:663
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:73
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:575
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition: coeffs.h:465
static BOOLEAN fa(leftv res, leftv args)
Definition: cohomo.cc:3764
static BOOLEAN pa(leftv res, leftv args)
Definition: cohomo.cc:3723
static BOOLEAN fb(leftv res, leftv args)
Definition: cohomo.cc:3711
static BOOLEAN pb(leftv res, leftv args)
Definition: cohomo.cc:3747
#define Print
Definition: emacs.cc:80
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
CanonicalForm res
Definition: facAbsFact.cc:60
CanonicalForm factor
Definition: facAbsFact.cc:97
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define VAR
Definition: globaldefs.h:5
STATIC_VAR Poly * h
Definition: janet.cc:971
poly p_ChineseRemainder(poly *xx, mpz_ptr *x, mpz_ptr *q, int rl, mpz_ptr *C, const ring R)
number nlModP(number q, const coeffs, const coeffs Zp)
Definition: longrat.cc:1577
#define SR_INT
Definition: longrat.h:67
static FORCE_INLINE BOOLEAN nlIsInteger(number q, const coeffs r)
Definition: longrat.h:94
#define assume(x)
Definition: mod2.h:389
#define pIter(p)
Definition: monomials.h:37
#define pNext(p)
Definition: monomials.h:36
#define p_GetCoeff(p, r)
Definition: monomials.h:50
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:357
The main handler for Singular numbers which are suitable for Singular polynomials.
const char *const nDivBy0
Definition: numbers.h:89
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
#define omAlloc0Bin(bin)
Definition: omAllocDecl.h:206
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:259
#define omGetSpecBin(size)
Definition: omBin.h:11
#define NULL
Definition: omList.c:12
omBin_t * omBin
Definition: omStructs.h:12
poly p_Diff(poly a, int k, const ring r)
Definition: p_polys.cc:1898
poly p_Farey(poly p, number N, const ring r)
Definition: p_polys.cc:54
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1297
poly p_Div_nn(poly p, const number n, const ring r)
Definition: p_polys.cc:1505
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3813
void p_Norm(poly p1, const ring r)
Definition: p_polys.cc:3719
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2845
int p_Var(poly m, const ring r)
Definition: p_polys.cc:4656
poly p_One(const ring r)
Definition: p_polys.cc:1313
poly p_Sub(poly p1, poly p2, const ring r)
Definition: p_polys.cc:1990
poly p_NSet(number n, const ring r)
returns the poly representing the number n, destroys n
Definition: p_polys.cc:1473
const char * p_Read(const char *st, poly &rc, const ring r)
Definition: p_polys.cc:1370
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition: p_polys.cc:4512
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1105
static int pLength(poly a)
Definition: p_polys.h:188
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:934
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1112
void p_Write(poly p, ring lmRing, ring tailRing)
Definition: polys0.cc:342
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent @Note: VarOffset encodes the position in p->exp
Definition: p_polys.h:486
static void p_ExpVectorDiff(poly pr, poly p1, poly p2, const ring r)
Definition: p_polys.h:1472
void p_String0Long(const poly p, ring lmRing, ring tailRing)
print p in a long way
Definition: polys0.cc:203
void p_String0Short(const poly p, ring lmRing, ring tailRing)
print p in a short way, if possible
Definition: polys0.cc:184
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:231
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:410
static BOOLEAN p_LmIsConstantComp(const poly p, const ring r)
Definition: p_polys.h:1004
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent @Note: the integer VarOffset encodes:
Definition: p_polys.h:467
static BOOLEAN p_LmIsConstant(const poly p, const ring r)
Definition: p_polys.h:1021
static BOOLEAN p_IsOne(const poly p, const ring R)
either poly(1) or gen(k)?!
Definition: p_polys.h:1969
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition: p_polys.h:1962
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:899
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1149
static poly p_LmFreeAndNext(poly p, ring)
Definition: p_polys.h:709
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:844
static long p_Totaldegree(poly p, const ring r)
Definition: p_polys.h:1505
#define p_Test(p, r)
Definition: p_polys.h:159
#define __p_Mult_nn(p, n, r)
Definition: p_polys.h:969
void p_wrp(poly p, ring lmRing, ring tailRing)
Definition: polys0.cc:373
poly prMapR(poly src, nMapFunc nMap, ring src_r, ring dest_r)
Definition: prCopy.cc:45
poly prCopyR(poly p, ring src_r, ring dest_r)
Definition: prCopy.cc:34
void StringAppendS(const char *st)
Definition: reporter.cc:107
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
void rDelete(ring r)
unconditionally deletes fields in r
Definition: ring.cc:450
BOOLEAN rEqual(ring r1, ring r2, BOOLEAN qr)
returns TRUE, if r1 equals r2 FALSE, otherwise Equality is determined componentwise,...
Definition: ring.cc:1746
static ring rIncRefCnt(ring r)
Definition: ring.h:837
static char * rRingVar(short i, const ring r)
Definition: ring.h:577
static void rDecRefCnt(ring r)
Definition: ring.h:838
static BOOLEAN rCanShortOut(const ring r)
Definition: ring.h:586
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:592
#define R
Definition: sirandom.c:27
#define A
Definition: sirandom.c:24
#define Q
Definition: sirandom.c:26
static poly convert(const number &n)
Definition: transext.cc:2326
#define loop
Definition: structs.h:75
#define SR_HDL(A)
Definition: tgb.cc:35
static int ntSize(number a, const coeffs cf)
Definition: transext.cc:1903
static BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs r)
Definition: transext.cc:140
static number ntParameter(const int iParameter, const coeffs cf)
return the specified parameter as a number in the given trans.ext.
Definition: transext.cc:2286
static void ntWriteLong(number a, const coeffs cf)
Definition: transext.cc:1641
static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void *param)
Definition: transext.cc:1727
static void heuristicGcdCancellation(number a, const coeffs cf)
forward declarations
Definition: transext.cc:1398
static const char * ntRead(const char *s, number *a, const coeffs cf)
Definition: transext.cc:1691
number ntDiff(number a, number d, const coeffs cf)
Definition: transext.cc:897
nMapFunc ntSetMap(const coeffs src, const coeffs dst)
Get a mapping function from src into the domain of this type (n_transExt)
Definition: transext.cc:2170
static number ntSub(number a, number b, const coeffs cf)
Definition: transext.cc:1034
static number ntCopyAlg(number a, const coeffs cf, const coeffs dst)
Definition: transext.cc:2092
static long ntInt(number &a, const coeffs cf)
Definition: transext.cc:773
static number ntMapPP(number a, const coeffs src, const coeffs dst)
Definition: transext.cc:2133
static void ntDelete(number *a, const coeffs cf)
Definition: transext.cc:313
static BOOLEAN ntIsOne(number a, const coeffs cf)
Definition: transext.cc:669
static void ntNormalizeDen(fraction result, const ring R)
Definition: transext.cc:1200
poly gcd_over_Q(poly f, poly g, const ring r)
helper routine for calling singclap_gcd_r
Definition: transext.cc:275
BOOLEAN ntInitChar(coeffs cf, void *infoStruct)
Initialize the coeffs object.
Definition: transext.cc:2636
static BOOLEAN ntIsZero(number a, const coeffs cf)
Definition: transext.cc:306
static number ntGetNumerator(number &a, const coeffs cf)
TODO: normalization of a!?
Definition: transext.cc:506
static CanonicalForm ntConvSingNFactoryN(number n, BOOLEAN, const coeffs cf)
Definition: transext.cc:2268
static void ntClearDenominators(ICoeffsEnumerator &numberCollectionEnumerator, number &c, const coeffs cf)
Definition: transext.cc:2427
#define ntRing
Definition: transext.cc:80
static BOOLEAN ntIsMOne(number a, const coeffs cf)
Definition: transext.cc:678
static number ntMapZ0(number a, const coeffs src, const coeffs dst)
Definition: transext.cc:1954
static void ntMPZ(mpz_t m, number &n, const coeffs r)
Definition: transext.cc:2619
static void ntCoeffWrite(const coeffs cf, BOOLEAN details)
Definition: transext.cc:856
static number ntMap0P(number a, const coeffs src, const coeffs dst)
Definition: transext.cc:2109
#define NUMIS1(f)
TRUE iff num. represents 1.
Definition: transext.cc:67
static number ntGcd(number a, number b, const coeffs cf)
Definition: transext.cc:1830
static number ntConvFactoryNSingN(const CanonicalForm n, const coeffs cf)
Definition: transext.cc:2256
static void ntClearContent(ICoeffsEnumerator &numberCollectionEnumerator, number &c, const coeffs cf)
Definition: transext.cc:2334
#define ntCoeffs
Definition: transext.cc:86
static number ntInvers(number a, const coeffs cf)
Definition: transext.cc:1270
static number ntNormalizeHelper(number a, number b, const coeffs cf)
Definition: transext.cc:1749
static void ntPower(number a, int exp, number *b, const coeffs cf)
Definition: transext.cc:1336
number ntInit(long i, const coeffs cf)
Definition: transext.cc:704
#define COM(f)
Definition: transext.cc:69
static number ntMapP0(number a, const coeffs src, const coeffs dst)
Definition: transext.cc:1968
static number ntGetDenom(number &a, const coeffs cf)
TODO: normalization of a!?
Definition: transext.cc:567
static number ntMapUP(number a, const coeffs src, const coeffs dst)
Definition: transext.cc:2148
static void ntInpMult(number &a, number b, const coeffs cf)
Definition: transext.cc:1141
#define MULT_COMPLEXITY
complexity increase due to * and /
Definition: transext.cc:62
static void ntKillChar(coeffs cf)
Definition: transext.cc:2250
static void ntInpAdd(number &a, number b, const coeffs cf)
Definition: transext.cc:995
static void ntWriteShort(number a, const coeffs cf)
Definition: transext.cc:1666
static number ntDiv(number a, number b, const coeffs cf)
Definition: transext.cc:1224
static number ntCopy(number a, const coeffs cf)
Definition: transext.cc:372
static int ntParDeg(number a, const coeffs cf)
Definition: transext.cc:2277
static BOOLEAN ntGreaterZero(number a, const coeffs cf)
Definition: transext.cc:796
static number ntChineseRemainder(number *x, number *q, int rl, BOOLEAN, CFArray &inv_cache, const coeffs cf)
Definition: transext.cc:2571
static number ntGenMap(number a, const coeffs cf, const coeffs dst)
Definition: transext.cc:2013
int ntIsParam(number m, const coeffs cf)
if m == var(i)/1 => return i,
Definition: transext.cc:2308
static number ntMap00(number a, const coeffs src, const coeffs dst)
Definition: transext.cc:1930
static number ntMult(number a, number b, const coeffs cf)
Definition: transext.cc:1073
#define ADD_COMPLEXITY
complexity increase due to + and -
Definition: transext.cc:61
#define ntTest(a)
Definition: transext.cc:76
static number ntGenAlg(number a, const coeffs cf, const coeffs dst)
Definition: transext.cc:2099
static number ntNeg(number a, const coeffs cf)
this is in-place, modifies a
Definition: transext.cc:691
static void definiteGcdCancellation(number a, const coeffs cf, BOOLEAN simpleTestsHaveAlreadyBeenPerformed)
modifies a
Definition: transext.cc:1483
static number ntAdd(number a, number b, const coeffs cf)
Definition: transext.cc:954
static number ntFarey(number p, number n, const coeffs cf)
Definition: transext.cc:2601
static number ntCopyMap(number a, const coeffs cf, const coeffs dst)
Definition: transext.cc:1984
static void ntNormalize(number &a, const coeffs cf)
Definition: transext.cc:1701
static BOOLEAN ntEqual(number a, number b, const coeffs cf)
Definition: transext.cc:326
VAR omBin fractionObjectBin
Definition: transext.cc:89
static coeffs nCoeff_bottom(const coeffs r, int &height)
Definition: transext.cc:292
#define DIFF_COMPLEXITY
complexity increase due to diff
Definition: transext.cc:63
static BOOLEAN ntGreater(number a, number b, const coeffs cf)
Definition: transext.cc:806
static void handleNestedFractionsOverQ(fraction f, const coeffs cf)
Definition: transext.cc:413
static number ntExactDiv(number a, number b, const coeffs cf)
Definition: transext.cc:1719
#define BOUND_COMPLEXITY
maximum complexity of a number
Definition: transext.cc:64
static number ntInitMPZ(mpz_t m, const coeffs r)
Definition: transext.cc:2611
struct for passing initialization parameters to naInitChar
Definition: transext.h:88
int gcd(int a, int b)
Definition: walkSupport.cc:836