source: git/kernel/kbuckets.cc @ cfbee28

spielwiese
Last change on this file since cfbee28 was cfbee28, checked in by Frank Seelisch <seelisch@…>, 14 years ago
fixes trac ticket #200 git-svn-id: file:///usr/local/Singular/svn/trunk@12627 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$ */
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  //if (TEST_OPT_PROT) { Print("C(%d)",pl); }
432  kbTest(bucket);
433  return i;
434}
435
436void kBucketClear(kBucket_pt bucket, poly *p, int *length)
437{
438  int i = kBucketCanonicalize(bucket);
439  if (i > 0)
440  {
441  #ifdef USE_COEF_BUCKETS
442    MULTIPLY_BUCKET(bucket,i);
443    //bucket->coef[i]=NULL;
444  #endif
445    *p = bucket->buckets[i];
446    *length = bucket->buckets_length[i];
447    bucket->buckets[i] = NULL;
448    bucket->buckets_length[i] = 0;
449    bucket->buckets_used = 0;
450
451  }
452  else
453  {
454    *p = NULL;
455    *length = 0;
456  }
457}
458
459void kBucketSetLm(kBucket_pt bucket, poly lm)
460{
461  kBucketMergeLm(bucket);
462  pNext(lm) = NULL;
463  bucket->buckets[0] = lm;
464  bucket->buckets_length[0] = 1;
465}
466
467#else // HAVE_PSEUDO_BUCKETS
468
469void kBucketInit(kBucket_pt bucket, poly lm, int length)
470{
471  int i;
472
473  assume(bucket != NULL);
474  assume(length <= 0 || length == pLength(lm));
475
476  bucket->p = lm;
477  if (length <= 0) bucket->l = pLength(lm);
478  else bucket->l = length;
479
480}
481
482const poly kBucketGetLm(kBucket_pt bucket)
483{
484  return bucket->p;
485}
486
487poly kBucketExtractLm(kBucket_pt bucket)
488{
489  poly lm = bucket->p;
490  assume(pLength(bucket->p) == bucket->l);
491  pIter(bucket->p);
492  (bucket->l)--;
493  pNext(lm) = NULL;
494  return lm;
495}
496
497void kBucketClear(kBucket_pt bucket, poly *p, int *length)
498{
499  assume(pLength(bucket->p) == bucket->l);
500  *p = bucket->p;
501  *length = bucket->l;
502  bucket->p = NULL;
503  bucket->l = 0;
504}
505
506#endif // ! HAVE_PSEUDO_BUCKETS
507//////////////////////////////////////////////////////////////////////////
508///
509/// For changing the ring of the Bpoly to new_tailBin
510///
511void kBucketShallowCopyDelete(kBucket_pt bucket,
512                              ring new_tailRing, omBin new_tailBin,
513                              pShallowCopyDeleteProc p_shallow_copy_delete)
514{
515#ifndef HAVE_PSEUDO_BUCKETS
516  int i;
517
518  kBucketCanonicalize(bucket);
519  for (i=0; i<= bucket->buckets_used; i++)
520    if (bucket->buckets[i] != NULL)
521    {
522      MULTIPLY_BUCKET(bucket,i);
523      bucket->buckets[i] = p_shallow_copy_delete(bucket->buckets[i],
524                                                 bucket->bucket_ring,
525                                                 new_tailRing,
526                                                 new_tailBin);
527    }
528#else
529  bucket->p = p_shallow_copy_delete(p,
530                                    bucket_ring,
531                                    new_tailRing,
532                                    new_tailBin);
533#endif
534  bucket->bucket_ring = new_tailRing;
535}
536
537//////////////////////////////////////////////////////////////////////////
538///
539/// Bucket number i from bucket is out of length sync, resync
540///
541void kBucketAdjust(kBucket_pt bucket, int i) {
542
543  MULTIPLY_BUCKET(bucket,i);
544
545  int l1 = bucket->buckets_length[i];
546  poly p1 = bucket->buckets[i];
547  bucket->buckets[i] = NULL;
548  bucket->buckets_length[i] = 0;
549  i = pLogLength(l1);
550
551  while (bucket->buckets[i] != NULL)
552  {
553    //kbTest(bucket);
554    MULTIPLY_BUCKET(bucket,i);
555    p1 = p_Add_q(p1, bucket->buckets[i],
556                 l1, bucket->buckets_length[i], bucket->bucket_ring);
557    bucket->buckets[i] = NULL;
558    bucket->buckets_length[i] = 0;
559    i = pLogLength(l1);
560  }
561
562  bucket->buckets[i] = p1;
563  bucket->buckets_length[i]=l1;
564  if (i >= bucket->buckets_used)
565    bucket->buckets_used = i;
566  else
567    kBucketAdjustBucketsUsed(bucket);
568}
569
570//////////////////////////////////////////////////////////////////////////
571///
572/// Multiply Bucket by number ,i.e. Bpoly == n*Bpoly
573///
574void kBucket_Mult_n(kBucket_pt bucket, number n)
575{
576#ifndef HAVE_PSEUDO_BUCKETS
577  kbTest(bucket);
578  ring r=bucket->bucket_ring;
579  int i;
580
581  for (i=0; i<= bucket->buckets_used; i++)
582  {
583    if (bucket->buckets[i] != NULL)
584    {
585#ifdef USE_COEF_BUCKETS
586      if (i<coef_start)
587      {
588        bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
589#ifdef HAVE_RINGS
590        if (rField_is_Ring(r) && !(rField_is_Domain(r))) {
591          bucket->buckets_length[i] = pLength(bucket->buckets[i]);
592          kBucketAdjust(bucket, i);
593        }
594#endif
595      }
596      else
597      if (bucket->coef[i]!=NULL)
598      {
599        bucket->coef[i] = p_Mult_nn(bucket->coef[i],n,r);
600      }
601      else
602      {
603        bucket->coef[i] = p_NSet(n_Copy(n,r),r);
604      }
605#else
606      bucket->buckets[i] = p_Mult_nn(bucket->buckets[i], n, r);
607#ifdef HAVE_RINGS
608      if (rField_is_Ring(currRing) && !(rField_is_Domain(currRing))) {
609        bucket->buckets_length[i] = pLength(bucket->buckets[i]);
610        kBucketAdjust(bucket, i);
611      }
612#endif
613#endif
614    }
615  }
616  kbTest(bucket);
617#else
618  bucket->p = p_Mult_nn(bucket->p, n, bucket->bucket_ring);
619#endif
620}
621
622
623//////////////////////////////////////////////////////////////////////////
624///
625/// Add to Bucket a poly ,i.e. Bpoly == q+Bpoly
626///
627void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
628{
629  if (q == NULL) return;
630  assume(*l <= 0 || pLength(q) == *l);
631
632  int i, l1;
633  ring r = bucket->bucket_ring;
634
635  if (*l <= 0)
636  {
637    l1 = pLength(q);
638    *l = l1;
639  }
640  else
641    l1 = *l;
642
643  kBucketMergeLm(bucket);
644  kbTest(bucket);
645  i = pLogLength(l1);
646
647  while (bucket->buckets[i] != NULL)
648  {
649    //MULTIPLY_BUCKET(bucket,i);
650  #ifdef USE_COEF_BUCKETS
651    if (bucket->coef[i]!=NULL)
652    {
653      q = p_Plus_mm_Mult_qq(q, bucket->coef[i], bucket->buckets[i],
654                 l1, bucket->buckets_length[i], r);
655      p_Delete(&bucket->coef[i],r);
656      p_Delete(&bucket->buckets[i],r);
657    }
658    else
659    q = p_Add_q(q, bucket->buckets[i],
660                 l1, bucket->buckets_length[i], r);
661  #else
662    q = p_Add_q(q, bucket->buckets[i],
663                 l1, bucket->buckets_length[i], r);
664  #endif
665    bucket->buckets[i] = NULL;
666    bucket->buckets_length[i] = 0;
667    i = pLogLength(l1);
668    assume(i<= MAX_BUCKET);
669    assume(bucket->buckets_used<= MAX_BUCKET);
670  }
671
672  kbTest(bucket);
673  bucket->buckets[i] = q;
674  bucket->buckets_length[i]=l1;
675  if (i >= bucket->buckets_used)
676    bucket->buckets_used = i;
677  else
678    kBucketAdjustBucketsUsed(bucket);
679  kbTest(bucket);
680}
681
682
683
684//////////////////////////////////////////////////////////////////////////
685///
686/// Bpoly == Bpoly - m*p; where m is a monom
687/// Does not destroy p and m
688/// assume (*l <= 0 || pLength(p) == *l)
689void kBucket_Minus_m_Mult_p(kBucket_pt bucket, poly m, poly p, int *l,
690                            poly spNoether)
691{
692  assume(*l <= 0 || pLength(p) == *l);
693  int i, l1;
694  poly p1 = p;
695  poly last;
696  ring r = bucket->bucket_ring;
697
698  if (*l <= 0)
699  {
700    l1 = pLength(p1);
701    *l = l1;
702  }
703  else
704    l1 = *l;
705
706  if (m == NULL || p == NULL) return;
707
708#ifndef HAVE_PSEUDO_BUCKETS
709  kBucketMergeLm(bucket);
710  kbTest(bucket);
711  i = pLogLength(l1);
712
713  if ((i <= bucket->buckets_used) && (bucket->buckets[i] != NULL))
714  {
715    assume(pLength(bucket->buckets[i])==bucket->buckets_length[i]);
716//#ifdef USE_COEF_BUCKETS
717//     if(bucket->coef[i]!=NULL)
718//     {
719//       poly mult=p_Mult_mm(bucket->coef[i],m,r);
720//       bucket->coef[i]=NULL;
721//       p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], mult, p1,
722//                               bucket->buckets_length[i], l1,
723//                             spNoether, r);
724//     }
725//     else
726//#endif
727    MULTIPLY_BUCKET(bucket,i);
728    p1 = p_Minus_mm_Mult_qq(bucket->buckets[i], m, p1,
729                            bucket->buckets_length[i], l1,
730                            spNoether, r);
731    l1 = bucket->buckets_length[i];
732    bucket->buckets[i] = NULL;
733    bucket->buckets_length[i] = 0;
734#ifdef HAVE_RINGS
735    if (rField_is_Ring(currRing) && !(rField_is_Domain(currRing)))
736    {
737      l1 = pLength(p1);
738      assume(pLength(p1) == l1);
739    }
740#endif
741    i = pLogLength(l1);
742  }
743  else
744  {
745    pSetCoeff0(m, nNeg(pGetCoeff(m)));
746    if (spNoether != NULL)
747    {
748      l1 = -1;
749      p1 = r->p_Procs->pp_Mult_mm_Noether(p1, m, spNoether, l1, r, last);
750      i = pLogLength(l1);
751    }
752    else {
753      p1 = r->p_Procs->pp_Mult_mm(p1, m, r, last);
754#ifdef HAVE_RINGS
755      if (rField_is_Ring(currRing) && !(rField_is_Domain(currRing)))
756      {
757        l1 = pLength(p1);
758        i = pLogLength(l1);
759      }
760#endif
761    }
762    pSetCoeff0(m, nNeg(pGetCoeff(m)));
763  }
764
765  while (bucket->buckets[i] != NULL)
766  {
767    //kbTest(bucket);
768    MULTIPLY_BUCKET(bucket,i);
769    p1 = p_Add_q(p1, bucket->buckets[i],
770                 l1, bucket->buckets_length[i], r);
771    bucket->buckets[i] = NULL;
772    bucket->buckets_length[i] = 0;
773    i = pLogLength(l1);
774  }
775
776  bucket->buckets[i] = p1;
777  bucket->buckets_length[i]=l1;
778  if (i >= bucket->buckets_used)
779    bucket->buckets_used = i;
780  else
781    kBucketAdjustBucketsUsed(bucket);
782#else // HAVE_PSEUDO_BUCKETS
783  bucket->p = p_Minus_mm_Mult_qq(bucket->p, m,  p,
784                               bucket->l, l1,
785                               spNoether, r);
786#endif
787}
788
789//////////////////////////////////////////////////////////////////////////
790///
791/// Bpoly == Bpoly + m*p; where m is a monom
792/// Does not destroy p and m
793/// assume (l <= 0 || pLength(p) == l)
794void kBucket_Plus_mm_Mult_pp(kBucket_pt bucket, poly m, poly p, int l)
795{
796    assume((!rIsPluralRing(bucket->bucket_ring))||p_IsConstant(m, bucket->bucket_ring));
797  assume(l <= 0 || pLength(p) == l);
798  int i, l1;
799  poly p1 = p;
800  poly last;
801  ring r = bucket->bucket_ring;
802
803  if (m == NULL || p == NULL) return;
804
805  if (l <= 0)
806  {
807    l1 = pLength(p1);
808    l = l1;
809  }
810  else
811    l1 = l;
812
813  kBucketMergeLm(bucket);
814  kbTest(bucket);
815  i = pLogLength(l1);
816  number n=n_Init(1,r);
817  if (i <= bucket->buckets_used && bucket->buckets[i] != NULL)
818  {
819    //if (FALSE){
820    #ifdef USE_COEF_BUCKETS
821    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
822    {
823      number orig_coef=p_GetCoeff(bucket->coef[i],r);
824      //we take ownership:
825      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
826      number add_coef=n_Copy(p_GetCoeff(m,r),r);
827      number gcd=n_Gcd(add_coef, orig_coef,r);
828
829      if (!(n_IsOne(gcd,r)))
830      {
831        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
832        number add_coef2=n_IntDiv(add_coef, gcd,r);
833        n_Delete(&orig_coef,r);
834        n_Delete(&add_coef,r);
835        orig_coef=orig_coef2;
836        add_coef=add_coef2;
837
838        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
839        n_Delete(&n,r);
840        n=gcd;
841      }
842
843      //assume(n_IsOne(n,r));
844      number backup=p_GetCoeff(m,r);
845
846      p_SetCoeff0(m,add_coef,r);
847      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
848
849      n_Delete(&orig_coef,r);
850      p_Delete(&bucket->coef[i],r);
851
852      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
853                           bucket->buckets_length[i], l1, r);
854      l1=bucket->buckets_length[i];
855      bucket->buckets[i]=NULL;
856      bucket->buckets_length[i] = 0;
857      i = pLogLength(l1);
858      assume(l1==pLength(p1));
859
860      p_SetCoeff(m,backup,r); //deletes add_coef
861    }
862    else
863    #endif
864    {
865      MULTIPLY_BUCKET(bucket,i);
866      p1 = p_Plus_mm_Mult_qq(bucket->buckets[i], m, p1,
867                           bucket->buckets_length[i], l1, r);
868      l1 = bucket->buckets_length[i];
869      bucket->buckets[i] = NULL;
870      bucket->buckets_length[i] = 0;
871      i = pLogLength(l1);
872    }
873  }
874
875  else
876  {
877    #ifdef USE_COEF_BUCKETS
878    number swap_n=p_GetCoeff(m,r);
879
880    assume(n_IsOne(n,r));
881    p_SetCoeff0(m,n,r);
882    n=swap_n;
883    //p_SetCoeff0(n, swap_n, r);
884    //p_GetCoeff0(n, swap_n,r);
885    #endif
886    p1 = r->p_Procs->pp_Mult_mm(p1, m, r, last);
887    #ifdef USE_COEF_BUCKETS
888    //m may not be changed
889    p_SetCoeff(m,n_Copy(n,r),r);
890    #endif
891  }
892
893
894  while ((bucket->buckets[i] != NULL) && (p1!=NULL))
895  {
896    assume(i!=0);
897    #ifdef USE_COEF_BUCKETS
898    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
899    {
900      number orig_coef=p_GetCoeff(bucket->coef[i],r);
901      //we take ownership:
902      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
903      number add_coef=n_Copy(n,r);
904      number gcd=n_Gcd(add_coef, orig_coef,r);
905
906      if (!(n_IsOne(gcd,r)))
907      {
908        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
909        number add_coef2=n_IntDiv(add_coef, gcd,r);
910        n_Delete(&orig_coef,r);
911        n_Delete(&n,r);
912        n_Delete(&add_coef,r);
913        orig_coef=orig_coef2;
914        add_coef=add_coef2;
915        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
916        n=gcd;
917      }
918      //assume(n_IsOne(n,r));
919      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
920      p1=p_Mult_nn(p1,add_coef,r);
921
922      p1 = p_Add_q(p1, bucket->buckets[i],r);
923      l1=pLength(p1);
924
925      bucket->buckets[i]=NULL;
926      n_Delete(&orig_coef,r);
927      p_Delete(&bucket->coef[i],r);
928      //l1=bucket->buckets_length[i];
929      assume(l1==pLength(p1));
930    }
931    else
932    #endif
933    {
934      //don't do that, pull out gcd
935      #ifdef USE_COEF_BUCKETS
936      if(!(n_IsOne(n,r)))
937      {
938        p1=p_Mult_nn(p1, n, r);
939        n_Delete(&n,r);
940        n=n_Init(1,r);
941      }
942      #endif
943      MULTIPLY_BUCKET(bucket,i);
944      p1 = p_Add_q(p1, bucket->buckets[i],
945                 l1, bucket->buckets_length[i], r);
946      bucket->buckets[i] = NULL;
947      bucket->buckets_length[i] = 0;
948    }
949    i = pLogLength(l1);
950  }
951
952  bucket->buckets[i] = p1;
953#ifdef USE_COEF_BUCKETS
954  assume(bucket->coef[i]==NULL);
955
956  if (!(n_IsOne(n,r)))
957  {
958    bucket->coef[i]=p_NSet(n,r);
959  }
960  else
961  {
962    bucket->coef[i]=NULL;
963    n_Delete(&n,r);
964  }
965
966  if ((p1==NULL) && (bucket->coef[i]!=NULL))
967    p_Delete(&bucket->coef[i],r);
968#endif
969  bucket->buckets_length[i]=l1;
970  if (i >= bucket->buckets_used)
971    bucket->buckets_used = i;
972  else
973    kBucketAdjustBucketsUsed(bucket);
974
975  kbTest(bucket);
976}
977
978poly kBucket_ExtractLarger(kBucket_pt bucket, poly q, poly append)
979{
980  if (q == NULL) return append;
981  poly lm;
982  loop
983  {
984    lm = kBucketGetLm(bucket);
985    if (lm == NULL) return append;
986    if (p_LmCmp(lm, q, bucket->bucket_ring) == 1)
987    {
988      lm = kBucketExtractLm(bucket);
989      pNext(append) = lm;
990      pIter(append);
991    }
992    else
993    {
994      return append;
995    }
996  }
997}
998
999/////////////////////////////////////////////////////////////////////////////
1000//
1001// Extract all monomials from bucket with component comp
1002// Return as a polynomial *p with length *l
1003// In other words, afterwards
1004// Bpoly = Bpoly - (poly consisting of all monomials with component comp)
1005// and components of monomials of *p are all 0
1006//
1007
1008// Hmm... for now I'm too lazy to implement those independent of currRing
1009// But better declare it extern than including polys.h
1010extern void pTakeOutComp(poly *p, Exponent_t comp, poly *q, int *lq);
1011void pDecrOrdTakeOutComp(poly *p, Exponent_t comp, long order,
1012                         poly *q, int *lq);
1013
1014void kBucketTakeOutComp(kBucket_pt bucket,
1015                        Exponent_t comp,
1016                        poly *r_p, int *l)
1017{
1018  poly p = NULL, q;
1019  int i, lp = 0, lq;
1020
1021#ifndef HAVE_PSEUDO_BUCKETS
1022  kBucketMergeLm(bucket);
1023  for (i=1; i<=bucket->buckets_used; i++)
1024  {
1025    if (bucket->buckets[i] != NULL)
1026    {
1027      MULTIPLY_BUCKET(bucket,i);
1028      pTakeOutComp(&(bucket->buckets[i]), comp, &q, &lq);
1029      if (q != NULL)
1030      {
1031        assume(pLength(q) == lq);
1032        bucket->buckets_length[i] -= lq;
1033        assume(pLength(bucket->buckets[i]) == bucket->buckets_length[i]);
1034        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1035      }
1036    }
1037  }
1038  kBucketAdjustBucketsUsed(bucket);
1039#else
1040  pTakeOutComp(&(bucket->p), comp, &p, &lp);
1041  (bucket->l) -= lp;
1042#endif
1043  *r_p = p;
1044  *l = lp;
1045
1046  kbTest(bucket);
1047}
1048
1049void kBucketDecrOrdTakeOutComp(kBucket_pt bucket,
1050                               Exponent_t comp, long order,
1051                               poly *r_p, int *l)
1052{
1053  poly p = NULL, q;
1054  int i, lp = 0, lq;
1055
1056#ifndef HAVE_PSEUDO_BUCKETS
1057  kBucketMergeLm(bucket);
1058  for (i=1; i<=bucket->buckets_used; i++)
1059  {
1060    if (bucket->buckets[i] != NULL)
1061    {
1062      MULTIPLY_BUCKET(bucket,i);
1063      pDecrOrdTakeOutComp(&(bucket->buckets[i]), comp, order, &q, &lq);
1064      if (q != NULL)
1065      {
1066        bucket->buckets_length[i] -= lq;
1067        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1068      }
1069    }
1070  }
1071  kBucketAdjustBucketsUsed(bucket);
1072#else
1073  pDecrOrdTakeOutComp(&(bucket->p), comp, order, &p, &lp);
1074  (bucket->l) -= lp;
1075#endif
1076
1077  *r_p = p;
1078  *l = lp;
1079}
1080
1081/////////////////////////////////////////////////////////////////////////////
1082// Reduction of Bpoly with a given poly
1083//
1084
1085extern int ksCheckCoeff(number *a, number *b);
1086
1087number kBucketPolyRed(kBucket_pt bucket,
1088                      poly p1, int l1,
1089                      poly spNoether)
1090{
1091  assume((!rIsPluralRing(bucket->bucket_ring))||p_LmEqual(p1,kBucketGetLm(bucket), bucket->bucket_ring));
1092  assume(p1 != NULL &&
1093         p_DivisibleBy(p1,  kBucketGetLm(bucket), bucket->bucket_ring));
1094  assume(pLength(p1) == (int) l1);
1095
1096  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
1097  BOOLEAN reset_vec=FALSE;
1098  number rn;
1099
1100  if(a1==NULL)
1101  {
1102    p_DeleteLm(&lm, bucket->bucket_ring);
1103    return nInit(1);
1104  }
1105
1106  if (! nIsOne(pGetCoeff(p1)))
1107  {
1108    number an = pGetCoeff(p1), bn = pGetCoeff(lm);
1109//StringSetS("##### an = "); nWrite(an); PrintS(StringAppend("\n"));
1110//StringSetS("##### bn = "); nWrite(bn); PrintS(StringAppend("\n"));
1111    int ct = ksCheckCoeff(&an, &bn);
1112    p_SetCoeff(lm, bn, bucket->bucket_ring);
1113    if ((ct == 0) || (ct == 2))
1114    {
1115      kBucket_Mult_n(bucket, an);
1116      lm = p_Mult_nn(lm, an, bucket->bucket_ring);
1117    }
1118    rn = an;
1119  }
1120  else
1121  {
1122    rn = nInit(1);
1123  }
1124
1125  if (p_GetComp(p1, bucket->bucket_ring) != p_GetComp(lm, bucket->bucket_ring))
1126  {
1127    p_SetCompP(a1, p_GetComp(lm, bucket->bucket_ring), bucket->bucket_ring);
1128    reset_vec = TRUE;
1129    p_SetComp(lm, p_GetComp(p1, bucket->bucket_ring), bucket->bucket_ring);
1130    p_Setm(lm, bucket->bucket_ring);
1131  }
1132
1133  p_ExpVectorSub(lm, p1, bucket->bucket_ring);
1134  l1--;
1135
1136  assume(l1==pLength(a1));
1137  BOOLEAN backuped=FALSE;
1138  number coef;
1139  #if 0
1140  //@Viktor, don't ignore coefficients on monomials
1141  if(l1==1) {
1142
1143    //if (rField_is_Q(bucket->bucket_ring)) {
1144      //avoid this for function fields, as gcds are expensive at the moment
1145
1146
1147      coef=p_GetCoeff(a1,bucket->bucket_ring);
1148      lm=p_Mult_nn(lm, coef, bucket->bucket_ring);
1149      p_SetCoeff0(a1, n_Init(1,bucket->bucket_ring), bucket->bucket_ring);
1150      backuped=TRUE;
1151      //WARNING: not thread_safe
1152    //deletes coef as side effect
1153    //}
1154  }
1155  #endif
1156  kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1157
1158  if (backuped)
1159    p_SetCoeff0(a1,coef,bucket->bucket_ring);
1160  p_DeleteLm(&lm, bucket->bucket_ring);
1161  if (reset_vec) p_SetCompP(a1, 0, bucket->bucket_ring);
1162  kbTest(bucket);
1163  return rn;
1164}
1165static BOOLEAN nIsPseudoUnit(number n, ring r)
1166{
1167  if (rField_is_Zp(r))
1168    return TRUE;
1169  if (r->parameter==NULL)
1170  {
1171    if (r->cf->nSize(n)==1)
1172      return TRUE;
1173    else
1174      return FALSE;
1175  }
1176  //if (r->parameter!=NULL)
1177  number one=n_Init(1,r);
1178  if (n_Equal(n,one,r))
1179  {
1180    n_Delete(&one,r);
1181    return TRUE;
1182  }
1183  n_Delete(&one,r);
1184  number minus_one=n_Init(-1,r);
1185  if (n_Equal(n,minus_one,r))
1186  {
1187    n_Delete(&minus_one,r);
1188    return TRUE;
1189  }
1190  return FALSE;
1191}
1192
1193void kBucketSimpleContent(kBucket_pt bucket)
1194{
1195  #ifdef USE_COEF_BUCKETS
1196  ring r=bucket->bucket_ring;
1197  int i;
1198  //PrintS("HHHHHHHHHHHHH");
1199  for (i=0;i<=MAX_BUCKET;i++)
1200  {
1201    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1202    //    PrintS("H2H2H2");
1203    if (i==0)
1204    {
1205      assume(bucket->buckets[i]==NULL);
1206    }
1207    if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]==NULL))
1208      return;
1209  }
1210  for (i=0;i<=MAX_BUCKET;i++)
1211  {
1212    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1213    //    PrintS("H2H2H2");
1214    if (i==0)
1215    {
1216      assume(bucket->buckets[i]==NULL);
1217    }
1218    if ((bucket->buckets[i]!=NULL)
1219    && (nIsPseudoUnit(p_GetCoeff(bucket->coef[i],r),r)))
1220      return;
1221  }
1222  //return;
1223
1224  number coef=n_Init(0,r);
1225  //ATTENTION: will not work correct for GB over ring
1226  //if (TEST_OPT_PROT)
1227  //    PrintS("CCCCCCCCCCCCC");
1228  for (i=MAX_BUCKET;i>=0;i--)
1229  {
1230    if (i==0)
1231    {
1232      assume(bucket->buckets[i]==NULL);
1233    }
1234    if (bucket->buckets[i]!=NULL)
1235    {
1236      assume(bucket->coef[i]!=NULL);
1237      assume(!(n_IsZero(pGetCoeff(bucket->coef[i]),r)));
1238
1239      //in this way it should crash on programming errors, yeah
1240      number temp=nGcd(coef, pGetCoeff(bucket->coef[i]),r);
1241      n_Delete(&coef,r );
1242      coef=temp;
1243      if (nIsPseudoUnit(coef,r))
1244      {
1245        n_Delete(&coef,r);
1246        return;
1247      }
1248      assume(!(n_IsZero(coef,r)));
1249    }
1250  }
1251  if (n_IsZero(coef,r)){
1252    n_Delete(&coef,r);
1253    return;
1254    }
1255  if (TEST_OPT_PROT)
1256    PrintS("S");
1257  for(i=0;i<=MAX_BUCKET;i++)
1258  {
1259    if (bucket->buckets[i]!=NULL)
1260    {
1261      assume(!(n_IsZero(coef,r)));
1262      assume(bucket->coef[i]!=NULL);
1263      number lc=p_GetCoeff(bucket->coef[i],r);
1264      p_SetCoeff(bucket->coef[i], n_IntDiv(lc,coef,r),r);
1265      assume(!(n_IsZero(p_GetCoeff(bucket->coef[i],r),r)));
1266    }
1267  }
1268  n_Delete(&coef,r);
1269  #endif
1270}
1271
1272
1273poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i)
1274{
1275  assume(bucket->buckets[i]!=NULL);
1276
1277  ring r=bucket->bucket_ring;
1278  poly p=bucket->buckets[i];
1279  bucket->buckets_length[i]--;
1280#ifdef USE_COEF_BUCKETS
1281  if (bucket->coef[i]!=NULL)
1282  {
1283    poly next=pNext(p);
1284    if (next==NULL)
1285    {
1286      MULTIPLY_BUCKET(bucket,i);
1287      p=bucket->buckets[i];
1288      bucket->buckets[i]=NULL;
1289      return p;
1290    }
1291    else
1292    {
1293      bucket->buckets[i]=next;
1294      number c=p_GetCoeff(bucket->coef[i],r);
1295      pNext(p)=NULL;
1296      p=p_Mult_nn(p,c,r);
1297      assume(p!=NULL);
1298      return p;
1299    }
1300  }
1301  else
1302#endif
1303  {
1304    bucket->buckets[i]=pNext(bucket->buckets[i]);
1305    pNext(p)=NULL;
1306    assume(p!=NULL);
1307    return p;
1308  }
1309}
Note: See TracBrowser for help on using the repository browser.