source: git/kernel/p_Mult_q.cc @ 171950

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