source: git/kernel/kbuckets.cc @ 009d80

spielwiese
Last change on this file since 009d80 was 009d80, checked in by Oliver Wienand <wienand@…>, 17 years ago
new compile switch: HAVE_RINGS at total there are now: HAVE_RINGS, HAVE_RING2TOM, HAVE_RINGMODN kspoly.cc, kstd1.cc, kstd2.cc, kutil.cc, kspoly.cc, polys.cc, ring.*: adaption of new switches and rField_is_Ring method numbers.*: introduced new operation: nDivBy(a, b) = Is a divisible by b? pDebig.cc, pInline1.h: use of new nDivBy method p_Mult_q: new routine for rings p_*, pp_*: template switch HAVE_ZERODIVISOR polys1.cc: poly exponentation, switches rmodulo*: nDivBy implementation, DBTest structs.h: nDivBy, ringtype, ringflaga, ringflagb git-svn-id: file:///usr/local/Singular/svn/trunk@10029 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 31.3 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: kbuckets.cc,v 1.32 2007-05-10 08:12:40 wienand Exp $ */
5
6#include "mod2.h"
7#include "structs.h"
8#include "omalloc.h"
9#include "p_polys.h"
10#include "febase.h"
11#include "pShallowCopyDelete.h"
12#include "numbers.h"
13#include "ring.h"
14#include "p_Procs.h"
15#include "kutil.h" // debug: KDEBUG
16#include "kbuckets.h"
17
18#ifdef HAVE_COEF_BUCKETS
19#define USE_COEF_BUCKETS
20#endif
21
22#ifdef USE_COEF_BUCKETS
23#ifdef HAVE_RINGS_OLD
24#define MULTIPLY_BUCKET(B,I) do                                        \
25  { if (B->coef[I]!=NULL)                                              \
26    {                                                                  \
27      assume(p_IsConstant(b->Coef[i],bucket->bucket->ring));           \
28      B->buckets[I]=p_Mult_q(B->buckets[I],B->coef[I],B->bucket_ring); \
29      B->coef[I]=NULL;                                                 \
30    }                                                                  \
31  } while(0)                                                           \
32    if (rField_is_Ring(currRing)) bucket->buckets_length[i] = pLength(bucket->buckets[i]);
33#else
34#define MULTIPLY_BUCKET(B,I) do                                        \
35  { if (B->coef[I]!=NULL)                                              \
36    {                                                                  \
37      B->buckets[I]=p_Mult_q(B->buckets[I],B->coef[I],B->bucket_ring); \
38      B->coef[I]=NULL;                                                 \
39    }                                                                  \
40  } while(0)
41#endif
42#else
43#define MULTIPLY_BUCKET(B,I)
44#endif
45static omBin kBucket_bin = omGetSpecBin(sizeof(kBucket));
46static int coef_start=1;
47//////////////////////////////////////////////////////////////////////////
48///
49/// Some internal stuff
50///
51
52// returns ceil(log_4(l))
53inline unsigned int pLogLength(unsigned int l)
54{
55  unsigned int i = 0;
56
57  if (l == 0) return 0;
58  l--;
59#ifdef BUCKET_TWO_BASE
60  while ((l = (l >> 1))) i++;
61#else
62  while ((l = (l >> 2))) i++;
63#endif
64  return i+1;
65}
66
67// returns ceil(log_4(pLength(p)))
68inline unsigned int pLogLength(poly p)
69{
70  return pLogLength((unsigned int) pLength(p));
71}
72
73#ifdef KDEBUG
74
75#ifndef HAVE_PSEUDO_BUCKETS
76BOOLEAN kbTest_i(kBucket_pt bucket, int i)
77{//sBucketSortMerge
78  #ifdef USE_COEF_BUCKETS
79  assume(bucket->coef[0]==NULL);
80  if ((bucket->coef[i]!=NULL) && (bucket->buckets[i]==NULL))
81  {
82    dReportError("Bucket %d coef not NULL", i);
83  }
84  if (bucket->coef[i]!=NULL)
85  {
86    assume(bucket->buckets[i]!=NULL);
87    _p_Test(bucket->coef[i],bucket->bucket_ring,PDEBUG);
88    }
89  #endif
90  pFalseReturn(p_Test(bucket->buckets[i], bucket->bucket_ring));
91  if (bucket->buckets_length[i] != pLength(bucket->buckets[i]))
92  {
93    dReportError("Bucket %d lengths difference should:%d has:%d",
94                 i, bucket->buckets_length[i], pLength(bucket->buckets[i]));
95  }
96  else if (i > 0 && (int) pLogLength(bucket->buckets_length[i]) > i)
97  {
98    dReportError("Bucket %d too long %d",
99                 i, bucket->buckets_length[i]);
100  }
101  if (i==0 && bucket->buckets_length[0] > 1)
102  {
103    dReportError("Bucket 0 too long");
104  }
105  return TRUE;
106}
107
108
109BOOLEAN kbTest(kBucket_pt bucket)
110{
111  #ifdef HAVE_COEF_BUCKETS
112  assume(bucket->coef[0]==NULL);
113  #endif
114  int i;
115  poly lm = bucket->buckets[0];
116
117  omCheckAddrBin(bucket, kBucket_bin);
118  if (bucket->bucket_ring!=currRing)
119  {
120     rTest(bucket->bucket_ring);
121  }
122  assume(bucket->buckets_used <= MAX_BUCKET);
123  if (! kbTest_i(bucket, 0)) return FALSE;
124  for (i=1; i<= (int) bucket->buckets_used; i++)
125  {
126    if (!kbTest_i(bucket, i)) return FALSE;
127    if (lm != NULL &&  bucket->buckets[i] != NULL
128        && p_LmCmp(lm, bucket->buckets[i], bucket->bucket_ring) != 1)
129    {
130      dReportError("Bucket %d larger or equal than lm", i);
131      if (p_LmCmp(lm, bucket->buckets[i], bucket->bucket_ring) ==0)
132        dReportError("Bucket %d equal to lm", i);
133      return FALSE;
134    }
135    if (!p_Test(bucket->buckets[i],bucket->bucket_ring))
136    {
137      dReportError("Bucket %d is not =0(4)", i);
138      return FALSE;
139    }
140  }
141
142  for (; i<=MAX_BUCKET; i++)
143  {
144    if (bucket->buckets[i] != NULL || bucket->buckets_length[i] != 0)
145    {
146      dReportError("Bucket %d not zero", i);
147      return FALSE;
148    }
149  }
150  for(i=0;i<=MAX_BUCKET;i++)
151  {
152    if (bucket->buckets[i]!=NULL)
153    {
154      int j;
155      for(j=i+1;j<=MAX_BUCKET;j++)
156      {
157        if (bucket->buckets[j]==bucket->buckets[i])
158        {
159          dReportError("Bucket %d %d equal", i,j);
160          return FALSE;
161        }
162      }
163    }
164    #ifdef HAVE_COEF_BUCKETS
165    if (bucket->coef[i]!=NULL)
166    {
167      int j;
168      for(j=i+1;j<=MAX_BUCKET;j++)
169      {
170        if (bucket->coef[j]==bucket->coef[i])
171        {
172          dReportError("internal coef %d %d equal", i,j);
173          return FALSE;
174        }
175      }
176    }
177    #endif
178  }
179  return TRUE;
180}
181
182#else // HAVE_PSEUDO_BUCKETS
183BOOLEAN kbTest(kBucket_pt bucket)
184{
185  return TRUE;
186}
187#endif // ! HAVE_PSEUDO_BUCKETS
188#endif // KDEBUG
189
190//////////////////////////////////////////////////////////////////////////
191///
192/// Creation/Destruction of buckets
193///
194
195kBucket_pt kBucketCreate(ring bucket_ring)
196{
197  assume(bucket_ring != NULL);
198  kBucket_pt bucket = (kBucket_pt) omAlloc0Bin(kBucket_bin);
199  bucket->bucket_ring = bucket_ring;
200  return bucket;
201}
202void kBucketDestroy(kBucket_pt *bucket_pt)
203{
204  omFreeBin(*bucket_pt, kBucket_bin);
205  *bucket_pt = NULL;
206}
207
208
209void kBucketDeleteAndDestroy(kBucket_pt *bucket_pt)
210{
211  kBucket_pt bucket = *bucket_pt;
212  kbTest(bucket);
213  int i;
214  for (i=0; i<= bucket->buckets_used; i++)
215  {
216
217    if (bucket->buckets[i] != NULL)
218    {
219      p_Delete(&(bucket->buckets[i]), bucket->bucket_ring);
220#ifdef USE_COEF_BUCKETS
221      if (bucket->coef[i]!=NULL)
222        p_Delete(&(bucket->coef[i]), bucket->bucket_ring);
223#endif
224    }
225  }
226  omFreeBin(bucket, kBucket_bin);
227  *bucket_pt = NULL;
228}
229
230/////////////////////////////////////////////////////////////////////////////
231// Convertion from/to Bpolys
232//
233#ifndef HAVE_PSEUDO_BUCKETS
234
235inline void kBucketMergeLm(kBucket_pt bucket)
236{
237  kbTest(bucket);
238  if (bucket->buckets[0] != NULL)
239  {
240    poly lm = bucket->buckets[0];
241    int i = 1;
242#ifdef BUCKET_TWO_BASE
243    int l = 2;
244    while ( bucket->buckets_length[i] >= l)
245    {
246      i++;
247      l = l << 1;
248    }
249#else
250    int l = 4;
251    while ( bucket->buckets_length[i] >= l)
252    {
253      i++;
254      l = l << 2;
255    }
256#endif
257#ifndef USE_COEF_BUCKETS
258    MULTIPLY_BUCKET(bucket,i);
259    pNext(lm) = bucket->buckets[i];
260    bucket->buckets[i] = lm;
261    bucket->buckets_length[i]++;
262    assume(i <= bucket->buckets_used+1);
263    if (i > bucket->buckets_used)  bucket->buckets_used = i;
264    bucket->buckets[0] = NULL;
265    bucket->buckets_length[0] = 0;
266    kbTest(bucket);
267#else
268    if (i > bucket->buckets_used)  bucket->buckets_used = i;
269    assume(i!=0);
270    if (bucket->buckets[i]!=NULL)
271    {
272       MULTIPLY_BUCKET(bucket,i);
273       pNext(lm) = bucket->buckets[i];
274       bucket->buckets[i] = lm;
275       bucket->buckets_length[i]++;
276       assume(i <= bucket->buckets_used+1);
277    }
278    else
279    {
280      #if 1
281       assume(bucket->buckets[i]==NULL);
282       assume(bucket->coef[0]==NULL);
283       assume(pLength(lm)==1);
284       assume(pNext(lm)==NULL);
285       number coef=p_GetCoeff(lm,bucket->bucket_ring);
286       //WARNING: not thread_safe
287       p_SetCoeff0(lm, n_Init(1,bucket->bucket_ring), bucket->bucket_ring);
288       bucket->buckets[i]=lm;
289       bucket->buckets_length[i]=1;
290       bucket->coef[i]=p_NSet(n_Copy(coef,bucket->bucket_ring),bucket->bucket_ring);
291
292       bucket->buckets[i]=lm;
293       bucket->buckets_length[i]=1;
294      #else
295       MULTIPLY_BUCKET(bucket,i);
296       pNext(lm) = bucket->buckets[i];
297       bucket->buckets[i] = lm;
298       bucket->buckets_length[i]++;
299       assume(i <= bucket->buckets_used+1);
300      #endif
301    }
302    bucket->buckets[0]=NULL;
303    bucket->buckets_length[0] = 0;
304    bucket->coef[0]=NULL;
305    kbTest(bucket);
306    #endif
307  }
308
309}
310
311BOOLEAN kBucketIsCleared(kBucket_pt bucket)
312{
313  int i;
314
315  for (i = 0;i<=MAX_BUCKET;i++)
316  {
317    if (bucket->buckets[i] != NULL) return FALSE;
318    #ifdef HAVE_COEF_BUCKETS
319    if (bucket->coef[i] != NULL) return FALSE;
320    #endif
321    if (bucket->buckets_length[i] != 0) return FALSE;
322  }
323  return TRUE;
324}
325
326void kBucketInit(kBucket_pt bucket, poly lm, int length)
327{
328  //assume(false);
329  assume(bucket != NULL);
330  assume(length <= 0 || length == pLength(lm));
331  assume(kBucketIsCleared(bucket));
332
333  if (lm == NULL) return;
334
335  if (length <= 0)
336    length = pLength(lm);
337
338  bucket->buckets[0] = lm;
339  #ifdef HAVE_COEF_BUCKETS
340  assume(bucket->coef[0]==NULL);
341  #endif
342  #ifdef USE_COEF_BUCKETS
343  bucket->coef[0]=NULL;
344  #endif
345  if (lm!=NULL)
346    bucket->buckets_length[0] = 1;
347  else
348    bucket->buckets_length[0]= 0;
349  if (length > 1)
350  {
351    unsigned int i = pLogLength(length-1);
352    bucket->buckets[i] = pNext(lm);
353    pNext(lm) = NULL;
354    bucket->buckets_length[i] = length-1;
355    bucket->buckets_used = i;
356  }
357  else
358  {
359    bucket->buckets_used = 0;
360  }
361}
362
363int kBucketCanonicalize(kBucket_pt bucket)
364{
365  assume(bucket->buckets_used<=MAX_BUCKET);
366  MULTIPLY_BUCKET(bucket,1);
367  kbTest(bucket);
368  poly p = bucket->buckets[1];
369  poly lm;
370  int pl = bucket->buckets_length[1];//, i;
371  int i;
372  bucket->buckets[1] = NULL;
373  bucket->buckets_length[1] = 0;
374  #ifdef USE_COEF_BUCKETS
375    assume(bucket->coef[1]==NULL);
376  #endif
377  ring r=bucket->bucket_ring;
378
379
380  for (i=1; i<=bucket->buckets_used; i++)
381  {
382  #ifdef USE_COEF_BUCKETS
383    if (bucket->coef[i]!=NULL)
384    {
385      assume(bucket->buckets[i]!=NULL);
386      p = p_Plus_mm_Mult_qq(p, bucket->coef[i], bucket->buckets[i],
387                 pl, bucket->buckets_length[i], r);
388      p_Delete(&bucket->coef[i],r);
389      p_Delete(&bucket->buckets[i],r);
390    }
391    else
392    p = p_Add_q(p, bucket->buckets[i],
393                 pl, bucket->buckets_length[i], r);
394  #else
395    p = p_Add_q(p, bucket->buckets[i],
396                 pl, bucket->buckets_length[i], r);
397  #endif
398    if (i==1) continue;
399    bucket->buckets[i] = NULL;
400    bucket->buckets_length[i] = 0;
401  }
402  #ifdef HAVE_COEF_BUCKETS
403  assume(bucket->coef[0]==NULL);
404  #endif
405  lm = bucket->buckets[0];
406  if (lm != NULL)
407  {
408    pNext(lm) = p;
409    p = lm;
410    pl++;
411    bucket->buckets[0] = NULL;
412    bucket->buckets_length[0] = 0;
413  }
414  if (pl > 0)
415  {
416    i = pLogLength(pl);
417    bucket->buckets[i] = p;
418    bucket->buckets_length[i] = pl;
419  }
420  else
421  {
422    i = 0;
423  }
424  bucket->buckets_used = i;
425  assume(bucket->buckets_used <= MAX_BUCKET);
426  #ifdef USE_COEF_BUCKETS
427    assume(bucket->coef[0]==NULL);
428    assume(bucket->coef[i]==NULL);
429  #endif
430  assume(pLength(p) == (int) pl);
431  kbTest(bucket);
432  return i;
433}
434
435void kBucketClear(kBucket_pt bucket, poly *p, int *length)
436{
437  int i = kBucketCanonicalize(bucket);
438  if (i > 0)
439  {
440  #ifdef USE_COEF_BUCKETS
441    MULTIPLY_BUCKET(bucket,i);
442    //bucket->coef[i]=NULL;
443#endif
444    *p = bucket->buckets[i];
445    *length = bucket->buckets_length[i];
446    bucket->buckets[i] = NULL;
447    bucket->buckets_length[i] = 0;
448    bucket->buckets_used = 0;
449
450  }
451  else
452  {
453    *p = NULL;
454    *length = 0;
455  }
456}
457
458void kBucketSetLm(kBucket_pt bucket, poly lm)
459{
460  kBucketMergeLm(bucket);
461  pNext(lm) = NULL;
462  bucket->buckets[0] = lm;
463  bucket->buckets_length[0] = 1;
464}
465
466#else // HAVE_PSEUDO_BUCKETS
467
468void kBucketInit(kBucket_pt bucket, poly lm, int length)
469{
470  int i;
471
472  assume(bucket != NULL);
473  assume(length <= 0 || length == pLength(lm));
474
475  bucket->p = lm;
476  if (length <= 0) bucket->l = pLength(lm);
477  else bucket->l = length;
478
479}
480
481const poly kBucketGetLm(kBucket_pt bucket)
482{
483  return bucket->p;
484}
485
486poly kBucketExtractLm(kBucket_pt bucket)
487{
488  poly lm = bucket->p;
489  assume(pLength(bucket->p) == bucket->l);
490  pIter(bucket->p);
491  (bucket->l)--;
492  pNext(lm) = NULL;
493  return lm;
494}
495
496void kBucketClear(kBucket_pt bucket, poly *p, int *length)
497{
498  assume(pLength(bucket->p) == bucket->l);
499  *p = bucket->p;
500  *length = bucket->l;
501  bucket->p = NULL;
502  bucket->l = 0;
503}
504
505#endif // ! HAVE_PSEUDO_BUCKETS
506//////////////////////////////////////////////////////////////////////////
507///
508/// For changing the ring of the Bpoly to new_tailBin
509///
510void kBucketShallowCopyDelete(kBucket_pt bucket,
511                              ring new_tailRing, omBin new_tailBin,
512                              pShallowCopyDeleteProc p_shallow_copy_delete)
513{
514#ifndef HAVE_PSEUDO_BUCKETS
515  int i;
516
517  kBucketCanonicalize(bucket);
518  for (i=0; i<= bucket->buckets_used; i++)
519    if (bucket->buckets[i] != NULL)
520    {
521      MULTIPLY_BUCKET(bucket,i);
522      bucket->buckets[i] = p_shallow_copy_delete(bucket->buckets[i],
523                                                 bucket->bucket_ring,
524                                                 new_tailRing,
525                                                 new_tailBin);
526    }
527#else
528  bucket->p = p_shallow_copy_delete(p,
529                                    bucket_ring,
530                                    new_tailRing,
531                                    new_tailBin);
532#endif
533  bucket->bucket_ring = new_tailRing;
534}
535
536//////////////////////////////////////////////////////////////////////////
537///
538/// Bucket number i from bucket is out of length sync, resync
539///
540void kBucketAdjust(kBucket_pt bucket, int i) {
541
542  MULTIPLY_BUCKET(bucket,i);
543
544  int l1 = bucket->buckets_length[i];
545  poly p1 = bucket->buckets[i];
546  bucket->buckets[i] = NULL;
547  bucket->buckets_length[i] = 0;
548  i = pLogLength(l1);
549
550  while (bucket->buckets[i] != NULL)
551  {
552    //kbTest(bucket);
553    MULTIPLY_BUCKET(bucket,i);
554    p1 = p_Add_q(p1, bucket->buckets[i],
555                 l1, bucket->buckets_length[i], bucket->bucket_ring);
556    bucket->buckets[i] = NULL;
557    bucket->buckets_length[i] = 0;
558    i = pLogLength(l1);
559  }
560
561  bucket->buckets[i] = p1;
562  bucket->buckets_length[i]=l1;
563  if (i >= bucket->buckets_used)
564    bucket->buckets_used = i;
565  else
566    kBucketAdjustBucketsUsed(bucket);
567}
568
569//////////////////////////////////////////////////////////////////////////
570///
571/// Multiply Bucket by number ,i.e. Bpoly == n*Bpoly
572///
573void kBucket_Mult_n(kBucket_pt bucket, number n)
574{
575#ifndef HAVE_PSEUDO_BUCKETS
576  kbTest(bucket);
577  ring r=bucket->bucket_ring;
578  int i;
579
580  for (i=0; i<= bucket->buckets_used; i++)
581  {
582    if (bucket->buckets[i] != NULL)
583    {
584#ifdef USE_COEF_BUCKETS
585      if (i<coef_start)
586        bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
587#ifdef HAVE_RING2TOM
588        if (r->cring == 1) {
589          bucket->buckets_length[i] = pLength(bucket->buckets[i]);
590          kBucketAdjust(bucket, i);
591        }
592#endif
593      else
594      if (bucket->coef[i]!=NULL)
595      {
596        bucket->coef[i] = p_Mult_nn(bucket->coef[i],n,r);
597#ifdef HAVE_RING2TOM_OLD
598        if (r->cring == 1 && (long) bucket->coef[i] == 0) {
599          bucket->coef[i] = NULL;
600          bucket->buckets[i] = NULL;
601          bucket->buckets_length[i] = 0;
602        }
603#endif
604      }
605      else
606      {
607        bucket->coef[i] = p_NSet(n_Copy(n,r),r);
608      }
609#else
610      bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
611#ifdef HAVE_RINGS
612      if (rField_is_Ring(currRing)) {
613        bucket->buckets_length[i] = pLength(bucket->buckets[i]);
614        kBucketAdjust(bucket, i);
615      }
616#endif
617#endif
618    }
619  }
620  kbTest(bucket);
621#else
622  bucket->p = p_Mult_nn(bucket->p, n, bucket->bucket_ring);
623#ifdef HAVE_RING2TOM_OLD
624  if (r->cring == 1) bucket->buckets_length[i] = pLength(bucket->buckets[i]);
625#endif
626#endif
627}
628
629
630//////////////////////////////////////////////////////////////////////////
631///
632/// Add to Bucket a poly ,i.e. Bpoly == q+Bpoly
633///
634void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
635{
636  if (q == NULL) return;
637  assume(*l <= 0 || pLength(q) == *l);
638
639  int i, l1;
640  ring r = bucket->bucket_ring;
641
642  if (*l <= 0)
643  {
644    l1 = pLength(q);
645    *l = l1;
646  }
647  else
648    l1 = *l;
649
650  kBucketMergeLm(bucket);
651  kbTest(bucket);
652  i = pLogLength(l1);
653
654  while (bucket->buckets[i] != NULL)
655  {
656    //MULTIPLY_BUCKET(bucket,i);
657  #ifdef USE_COEF_BUCKETS
658    if (bucket->coef[i]!=NULL)
659    {
660      q = p_Plus_mm_Mult_qq(q, bucket->coef[i], bucket->buckets[i],
661                 l1, bucket->buckets_length[i], r);
662      p_Delete(&bucket->coef[i],r);
663      p_Delete(&bucket->buckets[i],r);
664    }
665    else
666    q = p_Add_q(q, bucket->buckets[i],
667                 l1, bucket->buckets_length[i], r);
668  #else
669    q = p_Add_q(q, bucket->buckets[i],
670                 l1, bucket->buckets_length[i], r);
671  #endif
672    bucket->buckets[i] = NULL;
673    bucket->buckets_length[i] = 0;
674    i = pLogLength(l1);
675    assume(i<= MAX_BUCKET);
676    assume(bucket->buckets_used<= MAX_BUCKET);
677  }
678
679  kbTest(bucket);
680  bucket->buckets[i] = q;
681  bucket->buckets_length[i]=l1;
682  if (i >= bucket->buckets_used)
683    bucket->buckets_used = i;
684  else
685    kBucketAdjustBucketsUsed(bucket);
686  kbTest(bucket);
687}
688
689
690
691//////////////////////////////////////////////////////////////////////////
692///
693/// Bpoly == Bpoly - m*p; where m is a monom
694/// Does not destroy p and m
695/// assume (*l <= 0 || pLength(p) == *l)
696void kBucket_Minus_m_Mult_p(kBucket_pt bucket, poly m, poly p, int *l,
697                            poly spNoether)
698{
699  assume(*l <= 0 || pLength(p) == *l);
700  int i, l1;
701  poly p1 = p;
702  poly last;
703  ring r = bucket->bucket_ring;
704
705  if (*l <= 0)
706  {
707    l1 = pLength(p1);
708    *l = l1;
709  }
710  else
711    l1 = *l;
712
713  if (m == NULL || p == NULL) return;
714
715#ifndef HAVE_PSEUDO_BUCKETS
716  kBucketMergeLm(bucket);
717  kbTest(bucket);
718  i = pLogLength(l1);
719
720  if ((i <= bucket->buckets_used) && (bucket->buckets[i] != NULL))
721  {
722    assume(pLength(bucket->buckets[i])==bucket->buckets_length[i]);
723//#ifdef USE_COEF_BUCKETS
724//     if(bucket->coef[i]!=NULL)
725//     {
726//       poly mult=p_Mult_mm(bucket->coef[i],m,r);
727//       bucket->coef[i]=NULL;
728//       p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], mult, p1,
729//                               bucket->buckets_length[i], l1,
730//                             spNoether, r);
731//     }
732//     else
733//#endif
734    MULTIPLY_BUCKET(bucket,i);
735    p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], m, p1,
736                            bucket->buckets_length[i], l1,
737                            spNoether, r);
738    l1 = bucket->buckets_length[i];
739    bucket->buckets[i] = NULL;
740    bucket->buckets_length[i] = 0;
741#ifdef HAVE_RINGS
742    if (rField_is_Ring(currRing))
743    {
744      l1 = pLength(p1);
745      assume(pLength(p1) == l1);
746    }
747#endif
748    i = pLogLength(l1);
749  }
750  else
751  {
752    pSetCoeff0(m, nNeg(pGetCoeff(m)));
753    if (spNoether != NULL)
754    {
755      l1 = -1;
756      p1 = r->p_Procs->pp_Mult_mm_Noether(p1, m, spNoether, l1, r, last);
757      i = pLogLength(l1);
758    }
759    else {
760      p1 = r->p_Procs->pp_Mult_mm(p1, m, r, last);
761#ifdef HAVE_RINGS
762      if (rField_is_Ring(currRing))
763      {
764        l1 = pLength(p1);
765        i = pLogLength(l1);
766      }
767#endif
768    }
769    pSetCoeff0(m, nNeg(pGetCoeff(m)));
770  }
771
772  while (bucket->buckets[i] != NULL)
773  {
774    //kbTest(bucket);
775    MULTIPLY_BUCKET(bucket,i);
776    p1 = p_Add_q(p1, bucket->buckets[i],
777                 l1, bucket->buckets_length[i], r);
778    bucket->buckets[i] = NULL;
779    bucket->buckets_length[i] = 0;
780    i = pLogLength(l1);
781  }
782
783  bucket->buckets[i] = p1;
784  bucket->buckets_length[i]=l1;
785  if (i >= bucket->buckets_used)
786    bucket->buckets_used = i;
787  else
788    kBucketAdjustBucketsUsed(bucket);
789#else // HAVE_PSEUDO_BUCKETS
790  bucket->p = p_Minus_mm_Mult_qq(bucket->p, m,  p,
791                               bucket->l, l1,
792                               spNoether, r);
793#endif
794}
795
796//////////////////////////////////////////////////////////////////////////
797///
798/// Bpoly == Bpoly + m*p; where m is a monom
799/// Does not destroy p and m
800/// assume (l <= 0 || pLength(p) == l)
801void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
802{
803    assume((!rIsPluralRing(bucket->bucket_ring))||p_IsConstant(m, bucket->bucket_ring));
804  assume(l <= 0 || pLength(p) == l);
805  int i, l1;
806  poly p1 = p;
807  poly last;
808  ring r = bucket->bucket_ring;
809
810  if (m == NULL || p == NULL) return;
811
812  if (l <= 0)
813  {
814    l1 = pLength(p1);
815    l = l1;
816  }
817  else
818    l1 = l;
819
820  kBucketMergeLm(bucket);
821  kbTest(bucket);
822  i = pLogLength(l1);
823  number n=n_Init(1,r);
824  if (i <= bucket->buckets_used && bucket->buckets[i] != NULL)
825  {
826    //if (FALSE){
827    #ifdef USE_COEF_BUCKETS
828    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
829    {
830      number orig_coef=p_GetCoeff(bucket->coef[i],r);
831      //we take ownership:
832      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
833      number add_coef=n_Copy(p_GetCoeff(m,r),r);
834      number gcd=n_Gcd(add_coef, orig_coef,r);
835
836      if (!(n_IsOne(gcd,r)))
837      {
838        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
839        number add_coef2=n_IntDiv(add_coef, gcd,r);
840        n_Delete(&orig_coef,r);
841        n_Delete(&add_coef,r);
842        orig_coef=orig_coef2;
843        add_coef=add_coef2;
844
845        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
846        n_Delete(&n,r);
847        n=gcd;
848      }
849
850      //assume(n_IsOne(n,r));
851      number backup=p_GetCoeff(m,r);
852
853      p_SetCoeff0(m,add_coef,r);
854      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
855
856      n_Delete(&orig_coef,r);
857      p_Delete(&bucket->coef[i],r);
858
859      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
860                           bucket->buckets_length[i], l1, r);
861      l1=bucket->buckets_length[i];
862      bucket->buckets[i]=NULL;
863      bucket->buckets_length[i] = 0;
864      i = pLogLength(l1);
865      assume(l1==pLength(p1));
866
867      p_SetCoeff(m,backup,r); //deletes add_coef
868    }
869    else
870    #endif
871    {
872      MULTIPLY_BUCKET(bucket,i);
873      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
874                           bucket->buckets_length[i], l1, r);
875      l1 = bucket->buckets_length[i];
876      bucket->buckets[i] = NULL;
877      bucket->buckets_length[i] = 0;
878      i = pLogLength(l1);
879    }
880  }
881
882  else
883  {
884    #ifdef USE_COEF_BUCKETS
885    number swap_n=p_GetCoeff(m,r);
886
887    assume(n_IsOne(n,r));
888    p_SetCoeff0(m,n,r);
889    n=swap_n;
890    //p_SetCoeff0(n, swap_n, r);
891    //p_GetCoeff0(n, swap_n,r);
892    #endif
893    p1 = r->p_Procs->pp_Mult_mm(p1, m, r, last);
894    #ifdef USE_COEF_BUCKETS
895    //m may not be changed
896    p_SetCoeff(m,n_Copy(n,r),r);
897    #endif
898  }
899
900
901  while ((bucket->buckets[i] != NULL) && (p1!=NULL))
902  {
903    assume(i!=0);
904    #ifdef USE_COEF_BUCKETS
905    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
906    {
907      number orig_coef=p_GetCoeff(bucket->coef[i],r);
908      //we take ownership:
909      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
910      number add_coef=n_Copy(n,r);
911      number gcd=n_Gcd(add_coef, orig_coef,r);
912
913      if (!(n_IsOne(gcd,r)))
914      {
915        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
916        number add_coef2=n_IntDiv(add_coef, gcd,r);
917        n_Delete(&orig_coef,r);
918        n_Delete(&n,r);
919        n_Delete(&add_coef,r);
920        orig_coef=orig_coef2;
921        add_coef=add_coef2;
922        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
923        n=gcd;
924      }
925      //assume(n_IsOne(n,r));
926      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
927      p1=p_Mult_nn(p1,add_coef,r);
928
929      p1 = p_Add_q(p1, bucket->buckets[i],r);
930      l1=pLength(p1);
931
932      bucket->buckets[i]=NULL;
933      n_Delete(&orig_coef,r);
934      p_Delete(&bucket->coef[i],r);
935      //l1=bucket->buckets_length[i];
936      assume(l1==pLength(p1));
937    }
938    else
939    #endif
940    {
941      //don't do that, pull out gcd
942      #ifdef USE_COEF_BUCKETS
943      if(!(n_IsOne(n,r)))
944      {
945        p1=p_Mult_nn(p1, n, r);
946        n_Delete(&n,r);
947        n=n_Init(1,r);
948      }
949      #endif
950      MULTIPLY_BUCKET(bucket,i);
951      p1 = p_Add_q(p1, bucket->buckets[i],
952                 l1, bucket->buckets_length[i], r);
953      bucket->buckets[i] = NULL;
954      bucket->buckets_length[i] = 0;
955    }
956    i = pLogLength(l1);
957  }
958
959  bucket->buckets[i] = p1;
960#ifdef USE_COEF_BUCKETS
961  assume(bucket->coef[i]==NULL);
962
963  if (!(n_IsOne(n,r)))
964  {
965    bucket->coef[i]=p_NSet(n,r);
966  }
967  else
968  {
969    bucket->coef[i]=NULL;
970    n_Delete(&n,r);
971  }
972
973  if ((p1==NULL) && (bucket->coef[i]!=NULL))
974    p_Delete(&bucket->coef[i],r);
975#endif
976  bucket->buckets_length[i]=l1;
977  if (i >= bucket->buckets_used)
978    bucket->buckets_used = i;
979  else
980    kBucketAdjustBucketsUsed(bucket);
981
982  kbTest(bucket);
983}
984
985poly kBucket_ExtractLarger(kBucket_pt bucket, poly q, poly append)
986{
987  if (q == NULL) return append;
988  poly lm;
989  loop
990  {
991    lm = kBucketGetLm(bucket);
992    if (lm == NULL) return append;
993    if (p_LmCmp(lm, q, bucket->bucket_ring) == 1)
994    {
995      lm = kBucketExtractLm(bucket);
996      pNext(append) = lm;
997      pIter(append);
998    }
999    else
1000    {
1001      return append;
1002    }
1003  }
1004}
1005
1006/////////////////////////////////////////////////////////////////////////////
1007//
1008// Extract all monomials from bucket with component comp
1009// Return as a polynomial *p with length *l
1010// In other words, afterwards
1011// Bpoly = Bpoly - (poly consisting of all monomials with component comp)
1012// and components of monomials of *p are all 0
1013//
1014
1015// Hmm... for now I'm too lazy to implement those independent of currRing
1016// But better declare it extern than including polys.h
1017extern void pTakeOutComp(poly *p, Exponent_t comp, poly *q, int *lq);
1018void pDecrOrdTakeOutComp(poly *p, Exponent_t comp, Order_t order,
1019                         poly *q, int *lq);
1020
1021void kBucketTakeOutComp(kBucket_pt bucket,
1022                        Exponent_t comp,
1023                        poly *r_p, int *l)
1024{
1025  poly p = NULL, q;
1026  int i, lp = 0, lq;
1027
1028#ifndef HAVE_PSEUDO_BUCKETS
1029  kBucketMergeLm(bucket);
1030  for (i=1; i<=bucket->buckets_used; i++)
1031  {
1032    if (bucket->buckets[i] != NULL)
1033    {
1034      MULTIPLY_BUCKET(bucket,i);
1035      pTakeOutComp(&(bucket->buckets[i]), comp, &q, &lq);
1036      if (q != NULL)
1037      {
1038        assume(pLength(q) == lq);
1039        bucket->buckets_length[i] -= lq;
1040        assume(pLength(bucket->buckets[i]) == bucket->buckets_length[i]);
1041        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1042      }
1043    }
1044  }
1045  kBucketAdjustBucketsUsed(bucket);
1046#else
1047  pTakeOutComp(&(bucket->p), comp, &p, &lp);
1048  (bucket->l) -= lp;
1049#endif
1050  *r_p = p;
1051  *l = lp;
1052
1053  kbTest(bucket);
1054}
1055
1056void kBucketDecrOrdTakeOutComp(kBucket_pt bucket,
1057                               Exponent_t comp, Order_t order,
1058                               poly *r_p, int *l)
1059{
1060  poly p = NULL, q;
1061  int i, lp = 0, lq;
1062
1063#ifndef HAVE_PSEUDO_BUCKETS
1064  kBucketMergeLm(bucket);
1065  for (i=1; i<=bucket->buckets_used; i++)
1066  {
1067    if (bucket->buckets[i] != NULL)
1068    {
1069      MULTIPLY_BUCKET(bucket,i);
1070      pDecrOrdTakeOutComp(&(bucket->buckets[i]), comp, order, &q, &lq);
1071      if (q != NULL)
1072      {
1073        bucket->buckets_length[i] -= lq;
1074        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1075      }
1076    }
1077  }
1078  kBucketAdjustBucketsUsed(bucket);
1079#else
1080  pDecrOrdTakeOutComp(&(bucket->p), comp, order, &p, &lp);
1081  (bucket->l) -= lp;
1082#endif
1083
1084  *r_p = p;
1085  *l = lp;
1086}
1087
1088/////////////////////////////////////////////////////////////////////////////
1089// Reduction of Bpoly with a given poly
1090//
1091
1092extern int ksCheckCoeff(number *a, number *b);
1093
1094number kBucketPolyRed(kBucket_pt bucket,
1095                      poly p1, int l1,
1096                      poly spNoether)
1097{
1098  assume((!rIsPluralRing(bucket->bucket_ring))||p_LmEqual(p1,kBucketGetLm(bucket), bucket->bucket_ring));
1099  assume(p1 != NULL &&
1100         p_DivisibleBy(p1,  kBucketGetLm(bucket), bucket->bucket_ring));
1101  assume(pLength(p1) == (int) l1);
1102
1103  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
1104  BOOLEAN reset_vec=FALSE;
1105  number rn;
1106
1107  if(a1==NULL)
1108  {
1109    p_DeleteLm(&lm, bucket->bucket_ring);
1110    return nInit(1);
1111  }
1112
1113  if (! nIsOne(pGetCoeff(p1)))
1114  {
1115    number an = pGetCoeff(p1), bn = pGetCoeff(lm);
1116    int ct = ksCheckCoeff(&an, &bn);
1117    p_SetCoeff(lm, bn, bucket->bucket_ring);
1118    if ((ct == 0) || (ct == 2)) kBucket_Mult_n(bucket, an);
1119    rn = an;
1120  }
1121  else
1122  {
1123    rn = nInit(1);
1124  }
1125
1126  if (p_GetComp(p1, bucket->bucket_ring) != p_GetComp(lm, bucket->bucket_ring))
1127  {
1128    p_SetCompP(a1, p_GetComp(lm, bucket->bucket_ring), bucket->bucket_ring);
1129    reset_vec = TRUE;
1130    p_SetComp(lm, p_GetComp(p1, bucket->bucket_ring), bucket->bucket_ring);
1131    p_Setm(lm, bucket->bucket_ring);
1132  }
1133
1134  p_ExpVectorSub(lm, p1, bucket->bucket_ring);
1135  l1--;
1136
1137  assume(l1==pLength(a1));
1138  BOOLEAN backuped=FALSE;
1139  number coef;
1140  #if 0
1141  //@Viktor, don't ignore coefficients on monomials
1142  if(l1==1) {
1143
1144    //if (rField_is_Q(bucket->bucket_ring)) {
1145      //avoid this for function fields, as gcds are expensive at the moment
1146
1147
1148      coef=p_GetCoeff(a1,bucket->bucket_ring);
1149      lm=p_Mult_nn(lm, coef, bucket->bucket_ring);
1150      p_SetCoeff0(a1, n_Init(1,bucket->bucket_ring), bucket->bucket_ring);
1151      backuped=TRUE;
1152      //WARNING: not thread_safe
1153    //deletes coef as side effect
1154    //}
1155  }
1156  #endif
1157
1158  kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1159
1160  if (backuped)
1161    p_SetCoeff0(a1,coef,bucket->bucket_ring);
1162  p_DeleteLm(&lm, bucket->bucket_ring);
1163  if (reset_vec) p_SetCompP(a1, 0, bucket->bucket_ring);
1164  kbTest(bucket);
1165  return rn;
1166}
1167static BOOLEAN nIsPseudoUnit(number n, ring r)
1168{
1169  if (rField_is_Zp(r))
1170    return TRUE;
1171  if (r->parameter==NULL)
1172  {
1173    if (r->cf->nSize(n)==1)
1174      return TRUE;
1175    else
1176      return FALSE;
1177  }
1178  //if (r->parameter!=NULL)
1179  number one=n_Init(1,r);
1180  if (n_Equal(n,one,r))
1181  {
1182    n_Delete(&one,r);
1183    return TRUE;
1184  }
1185  n_Delete(&one,r);
1186  number minus_one=n_Init(-1,r);
1187  if (n_Equal(n,minus_one,r))
1188  {
1189    n_Delete(&minus_one,r);
1190    return TRUE;
1191  }
1192  return FALSE;
1193}
1194
1195void kBucketSimpleContent(kBucket_pt bucket)
1196{
1197  #ifdef USE_COEF_BUCKETS
1198  ring r=bucket->bucket_ring;
1199  int i;
1200  //PrintS("HHHHHHHHHHHHH");
1201  for (i=0;i<=MAX_BUCKET;i++)
1202  {
1203    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1204    //    PrintS("H2H2H2");
1205    if (i==0)
1206    {
1207      assume(bucket->buckets[i]==NULL);
1208    }
1209    if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]==NULL))
1210      return;
1211  }
1212  for (i=0;i<=MAX_BUCKET;i++)
1213  {
1214    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1215    //    PrintS("H2H2H2");
1216    if (i==0)
1217    {
1218      assume(bucket->buckets[i]==NULL);
1219    }
1220    if ((bucket->buckets[i]!=NULL)
1221    && (nIsPseudoUnit(p_GetCoeff(bucket->coef[i],r),r)))
1222      return;
1223  }
1224  //return;
1225
1226  number coef=n_Init(0,r);
1227  //ATTENTION: will not work correct for GB over ring
1228  //if (TEST_OPT_PROT)
1229  //    PrintS("CCCCCCCCCCCCC");
1230  for (i=MAX_BUCKET;i>=0;i--)
1231  {
1232    if (i==0)
1233    {
1234      assume(bucket->buckets[i]==NULL);
1235    }
1236    if (bucket->buckets[i]!=NULL)
1237    {
1238      assume(bucket->coef[i]!=NULL);
1239      assume(!(n_IsZero(pGetCoeff(bucket->coef[i]),r)));
1240
1241      //in this way it should crash on programming errors, yeah
1242      number temp=nGcd(coef, pGetCoeff(bucket->coef[i]),r);
1243      n_Delete(&coef,r );
1244      coef=temp;
1245      if (nIsPseudoUnit(coef,r))
1246      {
1247        n_Delete(&coef,r);
1248        return;
1249      }
1250      assume(!(n_IsZero(coef,r)));
1251    }
1252  }
1253  if (n_IsZero(coef,r)){
1254    n_Delete(&coef,r);
1255    return;
1256    }
1257  if (TEST_OPT_PROT)
1258    PrintS("S");
1259  for(i=0;i<=MAX_BUCKET;i++)
1260  {
1261    if (bucket->buckets[i]!=NULL)
1262    {
1263      assume(!(n_IsZero(coef,r)));
1264      assume(bucket->coef[i]!=NULL);
1265      number lc=p_GetCoeff(bucket->coef[i],r);
1266      p_SetCoeff(bucket->coef[i], n_IntDiv(lc,coef,r),r);
1267      assume(!(n_IsZero(p_GetCoeff(bucket->coef[i],r),r)));
1268    }
1269  }
1270  n_Delete(&coef,r);
1271  #endif
1272}
1273
1274
1275poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i)
1276{
1277  assume(bucket->buckets[i]!=NULL);
1278
1279  ring r=bucket->bucket_ring;
1280  poly p=bucket->buckets[i];
1281  bucket->buckets_length[i]--;
1282#ifdef USE_COEF_BUCKETS
1283  if (bucket->coef[i]!=NULL)
1284  {
1285    poly next=pNext(p);
1286    if (next==NULL)
1287    {
1288      MULTIPLY_BUCKET(bucket,i);
1289      p=bucket->buckets[i];
1290      bucket->buckets[i]=NULL;
1291      return p;
1292    }
1293    else
1294    {
1295      bucket->buckets[i]=next;
1296      number c=p_GetCoeff(bucket->coef[i],r);
1297      pNext(p)=NULL;
1298      p=p_Mult_nn(p,c,r);
1299      assume(p!=NULL);
1300      return p;
1301    }
1302  }
1303  else
1304#endif
1305  {
1306    bucket->buckets[i]=pNext(bucket->buckets[i]);
1307    pNext(p)=NULL;
1308    assume(p!=NULL);
1309    return p;
1310  }
1311}
Note: See TracBrowser for help on using the repository browser.