My Project
Loading...
Searching...
No Matches
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
31BOOLEAN 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
69static 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
100static 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
195static 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
223static 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)
313poly _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 && ((r->cf->extRing==NULL)||(r->cf->extRing->qideal!=NULL))
389 /* exclude trans. extensions: may contain rat.funct as cf */
390 && (((lq >= MIN_LENGTH_FACTORY)
391 && (r->cf->convSingNFactoryN!=ndConvSingNFactoryN))
393 && rField_is_Q(r))))
394 {
395 poly h=singclap_pmult(p,q,r);
396 if (!copy)
397 {
398 p_Delete(&p,r);
399 p_Delete(&q,r);
400 }
401 return h;
402 }
403 else
404 {
405 lp=pLength(p);
406 lq=pLength(q);
407 return _p_Mult_q_Bucket(p, lp, q, lq, copy, r);
408 }
409}
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:647
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:736
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_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:452
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:815
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
#define pSetCoeff0(p, n)
Definition: monomials.h:59
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
Definition: lq.h:40
CanonicalForm ndConvSingNFactoryN(number, BOOLEAN, const coeffs)
Definition: numbers.cc:313
#define NULL
Definition: omList.c:12
#define TEST_OPT_NOT_BUCKETS
Definition: options.h:106
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 int pLength(poly a)
Definition: p_polys.h:188
static poly pp_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:1029
#define p_LmCmpAction(p, q, r, actionE, actionG, actionS)
Definition: p_polys.h:1717
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:899
static poly p_LmFreeAndNext(poly p, ring)
Definition: p_polys.h:709
static void p_LmFree(poly p, ring)
Definition: p_polys.h:681
static poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq, const ring r)
Definition: p_polys.h:1181
BOOLEAN pHaveCommonMonoms(poly p, poly q)
Definition: pDebug.cc:178
#define p_Test(p, r)
Definition: p_polys.h:159
static BOOLEAN rField_is_Z(const ring r)
Definition: ring.h:509
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:500
static BOOLEAN rField_is_Domain(const ring r)
Definition: ring.h:487
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:506
#define rField_is_Ring(R)
Definition: ring.h:485