My Project
p_Mult_q.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /***************************************************************
5  * File: p_Mult_q.cc
6  * Purpose: multiplication of polynomials
7  * Author: obachman (Olaf Bachmann)
8  * Created: 8/00
9  *******************************************************************/
10 
11 #include "misc/auxiliary.h"
12 #include "factory/factory.h"
13 
14 #include "misc/options.h"
15 
16 #include "coeffs/numbers.h"
18 #include "polys/kbuckets.h"
19 #include "polys/clapsing.h"
20 
25 #include "polys/flintconv.h"
26 #include "polys/flint_mpoly.h"
27 
28 #include "p_Mult_q.h"
29 
30 
31 BOOLEAN pqLength(poly p, poly q, int &lp, int &lq, const int min)
32 {
33  int l = 0;
34 
35  do
36  {
37  if (p == NULL)
38  {
39  lp = l;
40  if (l < min)
41  {
42  if (q != NULL)
43  lq = l+1;
44  else
45  lq = l;
46  return FALSE;
47  }
48  lq = l + pLength(q);
49  return TRUE;
50  }
51  pIter(p);
52  if (q == NULL)
53  {
54  lq = l;
55  if (l < min)
56  {
57  lp = l+1;
58  return FALSE;
59  }
60  lp = l + 1 + pLength(p);
61  return TRUE;
62  }
63  pIter(q);
64  l++;
65  }
66  while (1);
67 }
68 
69 static void pqLengthApprox(poly p, poly q, int &lp, int &lq, const int min)
70 {
71  int l = 0;
72 
73  do
74  {
75  if (p == NULL)
76  {
77  lp=l;
78  lq=l+(q!=NULL);
79  return;
80  }
81  if (q == NULL) /* && p!=NULL */
82  {
83  lp=l+1;
84  lq=l;
85  return;
86  }
87  if (l>min) /* && p,q!=NULL */
88  {
89  lp=l; lq=l;
90  return;
91  }
92  pIter(p);
93  pIter(q);
94  l++;
95  }
96  while (1);
97 }
98 
99 
100 static poly _p_Mult_q_Bucket(poly p, const int lp,
101  poly q, const int lq,
102  const int copy, const ring r)
103 {
104  assume(p != NULL && pNext(p) != NULL && q != NULL && pNext(q) != NULL);
107  assume(lp >= 1 && lq >= 1);
108  p_Test(p, r);
109  p_Test(q, r);
110 
111  poly res = pp_Mult_mm(p,q,r); // holds initially q1*p
112  poly qq = pNext(q); // we iter of this
113  poly qn = pp_Mult_mm(qq, p,r); // holds p1*qi
114  poly pp = pNext(p); // used for Lm(qq)*pp
115  poly rr = res; // last monom which is surely not NULL
116  poly rn = pNext(res); // pNext(rr)
117  number n, n1;
118 
119  kBucket_pt bucket = kBucketCreate(r);
120 
121  // initialize bucket
122  kBucketInit(bucket, pNext(rn), lp - 2);
123  pNext(rn) = NULL;
124 
125  // now the main loop
126  Top:
127  if (rn == NULL) goto Smaller;
128  p_LmCmpAction(rn, qn, r, goto Equal, goto Greater, goto Smaller);
129 
130  Greater:
131  // rn > qn, so iter
132  rr = rn;
133  pNext(rn) = kBucketExtractLm(bucket);
134  pIter(rn);
135  goto Top;
136 
137  // rn < qn, append qn to rr, and compute next Lm(qq)*pp
138  Smaller:
139  pNext(rr) = qn;
140  rr = qn;
141  pIter(qn);
142  Work: // compute res + Lm(qq)*pp
143  if (rn == NULL)
144  {
145  pNext(rr) = pp_Mult_mm(pp, qq, r);
146  kBucketInit(bucket, pNext(pNext(rr)), lp - 2);
147  pNext(pNext(rr)) = NULL;
148  }
149  else
150  {
151  kBucketSetLm(bucket, rn);
152  kBucket_Plus_mm_Mult_pp(bucket, qq, pp, lp - 1);
153  pNext(rr) = kBucketExtractLm(bucket);
154  }
155 
156  pIter(qq);
157  if (qq == NULL) goto Finish;
158  rn = pNext(rr);
159  goto Top;
160 
161  Equal:
162  n1 = pGetCoeff(rn);
163  n = n_Add(n1, pGetCoeff(qn), r->cf);
164  n_Delete(&n1, r->cf);
165  if (n_IsZero(n, r->cf))
166  {
167  n_Delete(&n, r->cf);
168  p_LmFree(rn, r);
169  }
170  else
171  {
172  pSetCoeff0(rn, n);
173  rr = rn;
174  }
175  rn = kBucketExtractLm(bucket);
176  n_Delete(&pGetCoeff(qn),r->cf);
177  qn = p_LmFreeAndNext(qn, r);
178  goto Work;
179 
180  Finish:
181  assume(rr != NULL && pNext(rr) != NULL);
182  pNext(pNext(rr)) = kBucketClear(bucket);
183  kBucketDestroy(&bucket);
184 
185  if (!copy)
186  {
187  p_Delete(&p, r);
188  p_Delete(&q, r);
189  }
190  p_Test(res, r);
191  return res;
192 }
193 
194 #ifdef HAVE_RINGS
195 static poly _p_Mult_q_Normal_ZeroDiv(poly p, poly q, const int copy, const ring r)
196 {
197  assume(p != NULL && pNext(p) != NULL && q != NULL && pNext(q) != NULL);
199  p_Test(p, r);
200  p_Test(q, r);
201 
202  poly res = pp_Mult_mm(p,q,r); // holds initially q1*p
203  poly qq = pNext(q); // we iter of this
204 
205  while (qq != NULL)
206  {
207  res = p_Plus_mm_Mult_qq(res, qq, p, r);
208  pIter(qq);
209  }
210 
211  if (!copy)
212  {
213  p_Delete(&p, r);
214  p_Delete(&q, r);
215  }
216 
217  p_Test(res, r);
218 
219  return res;
220 }
221 #endif
222 
223 static poly _p_Mult_q_Normal(poly p, poly q, const int copy, const ring r)
224 {
225  assume(r != NULL);
226  assume(p != NULL && pNext(p) != NULL && q != NULL && pNext(q) != NULL);
227 #ifdef HAVE_RINGS
228  assume(nCoeff_is_Domain(r->cf));
229 #endif
230  pAssume1(! p_HaveCommonMonoms(p, q, r));
231  p_Test(p, r);
232  p_Test(q, r);
233 
234  poly res = pp_Mult_mm(p,q,r); // holds initially q1*p
235  poly qq = pNext(q); // we iter of this
236  poly qn = pp_Mult_mm(qq, p,r); // holds p1*qi
237  poly pp = pNext(p); // used for Lm(qq)*pp
238  poly rr = res; // last monom which is surely not NULL
239  poly rn = pNext(res); // pNext(rr)
240  number n, n1;
241 
242  // now the main loop
243  Top:
244  if (rn == NULL) goto Smaller;
245  p_LmCmpAction(rn, qn, r, goto Equal, goto Greater, goto Smaller);
246 
247  Greater:
248  // rn > qn, so iter
249  rr = rn;
250  pIter(rn);
251  goto Top;
252 
253  // rn < qn, append qn to rr, and compute next Lm(qq)*pp
254  Smaller:
255  pNext(rr) = qn;
256  rr = qn;
257  pIter(qn);
258 
259  Work: // compute res + Lm(qq)*pp
260  if (rn == NULL)
261  pNext(rr) = pp_Mult_mm(pp, qq, r);
262  else
263  {
264  pNext(rr) = p_Plus_mm_Mult_qq(rn, qq, pp, r);
265  }
266 
267  pIter(qq);
268  if (qq == NULL) goto Finish;
269  rn = pNext(rr);
270  goto Top;
271 
272  Equal:
273  n1 = pGetCoeff(rn);
274  n = n_Add(n1, pGetCoeff(qn), r->cf);
275  n_Delete(&n1, r->cf);
276  if (n_IsZero(n, r->cf))
277  {
278  n_Delete(&n, r->cf);
279  rn = p_LmFreeAndNext(rn, r);
280  }
281  else
282  {
283  pSetCoeff0(rn, n);
284  rr = rn;
285  pIter(rn);
286  }
287  n_Delete(&pGetCoeff(qn),r->cf);
288  qn = p_LmFreeAndNext(qn, r);
289  goto Work;
290 
291  Finish:
292  if (!copy)
293  {
294  p_Delete(&p, r);
295  p_Delete(&q, r);
296  }
297  p_Test(res, r);
298  return res;
299 }
300 
301 
302 // Use factory if min(pLength(p), pLength(q)) >= MIN_LENGTH_FACTORY (>MIN_LENGTH_BUCKET)
303 // Not thoroughly tested what is best
304 #define MIN_LENGTH_FACTORY 200
305 #define MIN_LENGTH_FACTORY_QQ 60
306 #define MIN_FLINT_QQ 10
307 #define MIN_FLINT_Zp 20
308 #define MIN_FLINT_Z 10
309 
310 /// Returns: p * q,
311 /// Destroys: if !copy then p, q
312 /// Assumes: pLength(p) >= 2 pLength(q) >=2, !rIsPluralRing(r)
313 poly _p_Mult_q(poly p, poly q, const int copy, const ring r)
314 {
315  assume(r != NULL);
316 #ifdef HAVE_RINGS
317  if (!nCoeff_is_Domain(r->cf))
318  return _p_Mult_q_Normal_ZeroDiv(p, q, copy, r);
319 #endif
320  int lp, lq, l;
321  poly pt;
322 
323  // MIN_LENGTH_FACTORY must be >= MIN_LENGTH_FACTORY_QQ, MIN_FLINT_QQ, MIN_FLINT_Zp 20
325 
326  if (lp < lq)
327  {
328  pt = p;
329  p = q;
330  q = pt;
331  l = lp;
332  lp = lq;
333  lq = l;
334  }
335  BOOLEAN pure_polys=(p_GetComp(p,r)==0) && (p_GetComp(q,r)==0);
336  #ifdef HAVE_FLINT
337  #if __FLINT_RELEASE >= 20503
338  if (lq>MIN_FLINT_QQ)
339  {
340  fmpq_mpoly_ctx_t ctx;
341  if (pure_polys && rField_is_Q(r) && !convSingRFlintR(ctx,r))
342  {
343  // lq is a lower bound for the length of p and q
344  poly res=Flint_Mult_MP(p,lq,q,lq,ctx,r);
345  if (!copy)
346  {
347  p_Delete(&p,r);
348  p_Delete(&q,r);
349  }
350  return res;
351  }
352  }
353  if (lq>MIN_FLINT_Zp)
354  {
355  nmod_mpoly_ctx_t ctx;
356  if (pure_polys && rField_is_Zp(r) && !convSingRFlintR(ctx,r))
357  {
358  // lq is a lower bound for the length of p and q
359  poly res=Flint_Mult_MP(p,lq,q,lq,ctx,r);
360  if (!copy)
361  {
362  p_Delete(&p,r);
363  p_Delete(&q,r);
364  }
365  return res;
366  }
367  }
368  if (lq>MIN_FLINT_Z)
369  {
370  fmpz_mpoly_ctx_t ctx;
371  if (pure_polys && rField_is_Z(r) && !convSingRFlintR(ctx,r))
372  {
373  // lq is a lower bound for the length of p and q
374  poly res=Flint_Mult_MP(p,lq,q,lq,ctx,r);
375  if (!copy)
376  {
377  p_Delete(&p,r);
378  p_Delete(&q,r);
379  }
380  return res;
381  }
382  }
383  #endif
384  #endif
386  return _p_Mult_q_Normal(p, q, copy, r);
387  else if (pure_polys
388  && (((lq >= MIN_LENGTH_FACTORY)
389  && (r->cf->convSingNFactoryN!=ndConvSingNFactoryN))
390  || ((lq >= MIN_LENGTH_FACTORY_QQ)
391  && rField_is_Q(r))))
392  {
393  poly h=singclap_pmult(p,q,r);
394  if (!copy)
395  {
396  p_Delete(&p,r);
397  p_Delete(&q,r);
398  }
399  return h;
400  }
401  else
402  {
403  lp=pLength(p);
404  lq=pLength(q);
405  return _p_Mult_q_Bucket(p, lp, q, lq, copy, r);
406  }
407 }
All the auxiliary stuff.
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
int l
Definition: cfEzgcd.cc:100
int p
Definition: cfModGcd.cc:4078
poly singclap_pmult(poly f, poly g, const ring r)
Definition: clapsing.cc:577
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of 'a' and 'b', i.e., a+b
Definition: coeffs.h:650
static FORCE_INLINE BOOLEAN nCoeff_is_Domain(const coeffs r)
returns TRUE, if r is a field or r has no zero divisors (i.e is a domain)
Definition: coeffs.h:739
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
Definition: coeffs.h:464
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:455
CanonicalForm res
Definition: facAbsFact.cc:60
CFArray copy(const CFList &list)
write elements of list into an array
static int min(int a, int b)
Definition: fast_mult.cc:268
static BOOLEAN Equal(number a, number b, const coeffs)
Definition: flintcf_Q.cc:324
This file is work in progress and currently not part of the official Singular.
static bool Greater(mono_type m1, mono_type m2)
STATIC_VAR Poly * h
Definition: janet.cc:971
void kBucketClear(kBucket_pt bucket, poly *p, int *length)
Definition: kbuckets.cc:521
void kBucketDestroy(kBucket_pt *bucket_pt)
Definition: kbuckets.cc:216
void kBucketInit(kBucket_pt bucket, poly lm, int length)
Definition: kbuckets.cc:493
poly kBucketExtractLm(kBucket_pt bucket)
Definition: kbuckets.cc:511
kBucket_pt kBucketCreate(const ring bucket_ring)
Creation/Destruction of buckets.
Definition: kbuckets.cc:209
void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
Bpoly == Bpoly + m*p; where m is a monom Does not destroy p and m assume (l <= 0 || pLength(p) == l)
Definition: kbuckets.cc:827
void kBucketSetLm(kBucket_pt bucket, poly lm)
#define assume(x)
Definition: mod2.h:389
#define p_GetComp(p, r)
Definition: monomials.h:64
#define pIter(p)
Definition: monomials.h:37
#define pNext(p)
Definition: monomials.h:36
#define pAssume1(cond)
Definition: monomials.h:171
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
#define pSetCoeff0(p, n)
Definition: monomials.h:59
Definition: lq.h:40
CanonicalForm ndConvSingNFactoryN(number, BOOLEAN, const coeffs)
Definition: numbers.cc:303
#define NULL
Definition: omList.c:12
#define TEST_OPT_NOT_BUCKETS
Definition: options.h:105
static void pqLengthApprox(poly p, poly q, int &lp, int &lq, const int min)
Definition: p_Mult_q.cc:69
#define MIN_LENGTH_FACTORY
Definition: p_Mult_q.cc:304
poly _p_Mult_q(poly p, poly q, const int copy, const ring r)
Returns: p * q, Destroys: if !copy then p, q Assumes: pLength(p) >= 2 pLength(q) >=2,...
Definition: p_Mult_q.cc:313
#define MIN_FLINT_Z
Definition: p_Mult_q.cc:308
BOOLEAN pqLength(poly p, poly q, int &lp, int &lq, const int min)
return TRUE and lp == pLength(p), lq == pLength(q), if min(pLength(p), pLength(q)) >= min FALSE if mi...
Definition: p_Mult_q.cc:31
#define MIN_FLINT_QQ
Definition: p_Mult_q.cc:306
static poly _p_Mult_q_Normal(poly p, poly q, const int copy, const ring r)
Definition: p_Mult_q.cc:223
#define MIN_LENGTH_FACTORY_QQ
Definition: p_Mult_q.cc:305
static poly _p_Mult_q_Bucket(poly p, const int lp, poly q, const int lq, const int copy, const ring r)
Definition: p_Mult_q.cc:100
static poly _p_Mult_q_Normal_ZeroDiv(poly p, poly q, const int copy, const ring r)
Definition: p_Mult_q.cc:195
#define MIN_FLINT_Zp
Definition: p_Mult_q.cc:307
#define MIN_LENGTH_BUCKET
Definition: p_Mult_q.h:21
static poly pp_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:1033
#define p_LmCmpAction(p, q, r, actionE, actionG, actionS)
Definition: p_polys.h:1721
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:903
static unsigned pLength(poly a)
Definition: p_polys.h:191
static poly p_LmFreeAndNext(poly p, ring)
Definition: p_polys.h:713
static void p_LmFree(poly p, ring)
Definition: p_polys.h:685
static poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq, const ring r)
Definition: p_polys.h:1185
BOOLEAN pHaveCommonMonoms(poly p, poly q)
Definition: pDebug.cc:175
#define p_Test(p, r)
Definition: p_polys.h:162
static BOOLEAN rField_is_Z(const ring r)
Definition: ring.h:510
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:501
static BOOLEAN rField_is_Domain(const ring r)
Definition: ring.h:488
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:507
#define rField_is_Ring(R)
Definition: ring.h:486