source: git/kernel/kbuckets.cc @ 3a0e1a

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