# source:git/libpolys/polys/operations/p_Mult_q.cc@5b53cb

spielwiese
Last change on this file since 5b53cb was 5b53cb, checked in by Hans Schoenemann <hannes@…>, 9 years ago
fix: typo
• Property mode set to `100644`
File size: 6.3 KB
Line
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
12
13#include <misc/auxiliary.h>
14
15#include <factory/factory.h>
16
17#include <misc/options.h>
18
19#include <polys/monomials/p_polys.h>
20#include <polys/templates/p_Procs.h>
21#include <polys/templates/p_MemCmp.h>
23#include <polys/templates/p_MemCopy.h>
24#include <polys/templates/p_Numbers.h>
25#include <polys/kbuckets.h>
26
27#include "p_Mult_q.h"
28
29
30BOOLEAN pqLength(poly p, poly q, int &lp, int &lq, const int min)
31{
32  int l = 0;
33
34  do
35  {
36    if (p == NULL)
37    {
38      lp = l;
39      if (l < min)
40      {
41        if (q != NULL)
42          lq = l+1;
43        else
44          lq = l;
45        return FALSE;
46      }
47      lq = l + pLength(q);
48      return TRUE;
49    }
50    pIter(p);
51    if (q == NULL)
52    {
53      lq = l;
54      if (l < min)
55      {
56        lp = l+1;
57        return FALSE;
58      }
59      lp = l + 1 + pLength(p);
60      return TRUE;
61    }
62    pIter(q);
63    l++;
64  }
65  while (1);
66}
67
68
69static poly _p_Mult_q_Bucket(poly p, const int lp,
70                             poly q, const int lq,
71                             const int copy, const ring r)
72{
73  assume(p != NULL && pNext(p) != NULL && q != NULL && pNext(q) != NULL);
74  pAssume1(! pHaveCommonMonoms(p, q));
75#ifdef HAVE_RINGS
76  assume(!rField_is_Ring(r) || rField_is_Domain(r));
77#endif
78  assume(lp >= 1 && lq >= 1);
79  p_Test(p, r);
80  p_Test(q, r);
81
82  poly res = pp_Mult_mm(p,q,r);     // holds initially q1*p
83  poly qq = pNext(q);               // we iter of this
84  poly qn = pp_Mult_mm(qq, p,r);    // holds p1*qi
85  poly pp = pNext(p);               // used for Lm(qq)*pp
86  poly rr = res;                    // last monom which is surely not NULL
87  poly rn = pNext(res);             // pNext(rr)
88  number n, n1;
89
90  kBucket_pt bucket = kBucketCreate(r);
91
92  // initialize bucket
93  kBucketInit(bucket, pNext(rn), lp - 2);
94  pNext(rn) = NULL;
95
96  // now the main loop
97  Top:
98  if (rn == NULL) goto Smaller;
99  p_LmCmpAction(rn, qn, r, goto Equal, goto Greater, goto Smaller);
100
101  Greater:
102  // rn > qn, so iter
103  rr = rn;
104  pNext(rn) = kBucketExtractLm(bucket);
105  pIter(rn);
106  goto Top;
107
108  // rn < qn, append qn to rr, and compute next Lm(qq)*pp
109  Smaller:
110  pNext(rr) = qn;
111  rr = qn;
112  pIter(qn);
113  Work: // compute res + Lm(qq)*pp
114  if (rn == NULL)
115  {
116    pNext(rr) = pp_Mult_mm(pp, qq, r);
117    kBucketInit(bucket, pNext(pNext(rr)), lp - 2);
118    pNext(pNext(rr)) = NULL;
119  }
120  else
121  {
122    kBucketSetLm(bucket, rn);
123    kBucket_Plus_mm_Mult_pp(bucket, qq, pp, lp - 1);
124    pNext(rr) = kBucketExtractLm(bucket);
125  }
126
127  pIter(qq);
128  if (qq == NULL) goto Finish;
129  rn = pNext(rr);
130  goto Top;
131
132  Equal:
133  n1 = pGetCoeff(rn);
134  n = n_Add(n1, pGetCoeff(qn), r->cf);
135  n_Delete(&n1, r->cf);
136  if (n_IsZero(n, r->cf))
137  {
138    n_Delete(&n, r->cf);
139    p_LmFree(rn, r);
140  }
141  else
142  {
143    pSetCoeff0(rn, n);
144    rr = rn;
145  }
146  rn = kBucketExtractLm(bucket);
147  n_Delete(&pGetCoeff(qn),r->cf);
148  qn = p_LmFreeAndNext(qn, r);
149  goto Work;
150
151  Finish:
152  assume(rr != NULL && pNext(rr) != NULL);
153  pNext(pNext(rr)) = kBucketClear(bucket);
154  kBucketDestroy(&bucket);
155
156  if (!copy)
157  {
158    p_Delete(&p, r);
159    p_Delete(&q, r);
160  }
161  p_Test(res, r);
162  return res;
163}
164
165#ifdef HAVE_RINGS
166static poly _p_Mult_q_Normal_ZeroDiv(poly p, poly q, const int copy, const ring r)
167{
168  assume(p != NULL && pNext(p) != NULL && q != NULL && pNext(q) != NULL);
169  pAssume1(! pHaveCommonMonoms(p, q));
170  p_Test(p, r);
171  p_Test(q, r);
172
173  poly res = pp_Mult_mm(p,q,r);     // holds initially q1*p
174  poly qq = pNext(q);               // we iter of this
175
176  while (qq != NULL)
177  {
178    res = p_Plus_mm_Mult_qq(res, qq, p, r);
179    pIter(qq);
180  }
181
182  if (!copy)
183  {
184    p_Delete(&p, r);
185    p_Delete(&q, r);
186  }
187
188  p_Test(res, r);
189
190  return res;
191}
192#endif
193
194static poly _p_Mult_q_Normal(poly p, poly q, const int copy, const ring r)
195{
196  assume(r != NULL);
197  assume(p != NULL && pNext(p) != NULL && q != NULL && pNext(q) != NULL);
198#ifdef HAVE_RINGS
199  assume(nCoeff_is_Domain(r->cf));
200#endif
201  pAssume1(! p_HaveCommonMonoms(p, q, r));
202  p_Test(p, r);
203  p_Test(q, r);
204
205  poly res = pp_Mult_mm(p,q,r);     // holds initially q1*p
206  poly qq = pNext(q);               // we iter of this
207  poly qn = pp_Mult_mm(qq, p,r);    // holds p1*qi
208  poly pp = pNext(p);               // used for Lm(qq)*pp
209  poly rr = res;                    // last monom which is surely not NULL
210  poly rn = pNext(res);             // pNext(rr)
211  number n, n1;
212
213  // now the main loop
214  Top:
215  if (rn == NULL) goto Smaller;
216  p_LmCmpAction(rn, qn, r, goto Equal, goto Greater, goto Smaller);
217
218  Greater:
219  // rn > qn, so iter
220  rr = rn;
221  pIter(rn);
222  goto Top;
223
224  // rn < qn, append qn to rr, and compute next Lm(qq)*pp
225  Smaller:
226  pNext(rr) = qn;
227  rr = qn;
228  pIter(qn);
229
230  Work: // compute res + Lm(qq)*pp
231  if (rn == NULL)
232    pNext(rr) = pp_Mult_mm(pp, qq, r);
233  else
234  {
235    pNext(rr) = p_Plus_mm_Mult_qq(rn, qq, pp, r);
236  }
237
238  pIter(qq);
239  if (qq == NULL) goto Finish;
240  rn = pNext(rr);
241  goto Top;
242
243  Equal:
244  n1 = pGetCoeff(rn);
245  n = n_Add(n1, pGetCoeff(qn), r->cf);
246  n_Delete(&n1, r->cf);
247  if (n_IsZero(n, r->cf))
248  {
249    n_Delete(&n, r->cf);
250    rn = p_LmFreeAndNext(rn, r);
251  }
252  else
253  {
254    pSetCoeff0(rn, n);
255    rr = rn;
256    pIter(rn);
257  }
258  n_Delete(&pGetCoeff(qn),r->cf);
259  qn = p_LmFreeAndNext(qn, r);
260  goto Work;
261
262  Finish:
263  if (!copy)
264  {
265    p_Delete(&p, r);
266    p_Delete(&q, r);
267  }
268  p_Test(res, r);
269  return res;
270}
271
272
273/// Returns:  p * q,
274/// Destroys: if !copy then p, q
275/// Assumes: pLength(p) >= 2 pLength(q) >=2
276poly _p_Mult_q(poly p, poly q, const int copy, const ring r)
277{
278  assume(r != NULL);
279#ifdef HAVE_RINGS
280  if (!nCoeff_is_Domain(r->cf))
281    return _p_Mult_q_Normal_ZeroDiv(p, q, copy, r);
282#endif
283  int lp, lq, l;
284  poly pt;
285
286  pqLength(p, q, lp, lq, MIN_LENGTH_BUCKET);
287
288  if (lp < lq)
289  {
290    pt = p;
291    p =  q;
292    q = pt;
293    l = lp;
294    lp = lq;
295    lq = l;
296  }
297  if (lq < MIN_LENGTH_BUCKET || TEST_OPT_NOT_BUCKETS)
298    return _p_Mult_q_Normal(p, q, copy, r);
299  else
300  {
301    assume(lp == pLength(p));
302    assume(lq == pLength(q));
303    return _p_Mult_q_Bucket(p, lp, q, lq, copy, r);
304  }
305}
Note: See TracBrowser for help on using the repository browser.