source: git/kernel/kbuckets.cc @ 8227a9

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