source: git/kernel/kbuckets.cc @ 42ddd0

spielwiese
Last change on this file since 42ddd0 was 42ddd0, checked in by Michael Brickenstein <bricken@…>, 18 years ago
*bricken: possible to compile without HAVE_COEF_BUCKETS git-svn-id: file:///usr/local/Singular/svn/trunk@8939 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.25 2006-02-13 19:52:26 bricken 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  #ifdef USE_COEF_BUCKETS
874  else
875  {
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
884    p1 = r->p_Procs->pp_Mult_mm(p1, m, r, last);
885    //m may not be changed
886    p_SetCoeff(m,n_Copy(n,r),r);
887  }
888  #endif
889
890  while ((bucket->buckets[i] != NULL) && (p1!=NULL))
891  {
892    assume(i!=0);
893    #ifdef USE_COEF_BUCKETS
894    if ((bucket->coef[i]!=NULL) &&(i>=coef_start))
895    {
896      number orig_coef=p_GetCoeff(bucket->coef[i],r);
897      //we take ownership:
898      p_SetCoeff0(bucket->coef[i],n_Init(0,r),r);
899      number add_coef=n_Copy(n,r);
900      number gcd=n_Gcd(add_coef, orig_coef,r);
901
902      if (!(n_IsOne(gcd,r)))
903      {
904        number orig_coef2=n_IntDiv(orig_coef,gcd,r);
905        number add_coef2=n_IntDiv(add_coef, gcd,r);
906        n_Delete(&orig_coef,r);
907        n_Delete(&n,r);
908        n_Delete(&add_coef,r);
909        orig_coef=orig_coef2;
910        add_coef=add_coef2;
911        //p_Mult_nn(bucket->buckets[i], orig_coef,r);
912        n=gcd;
913      }
914      //assume(n_IsOne(n,r));
915      bucket->buckets[i]=p_Mult_nn(bucket->buckets[i],orig_coef,r);
916      p1=p_Mult_nn(p1,add_coef,r);
917
918      p1 = p_Add_q(p1, bucket->buckets[i],r);
919      l1=pLength(p1);
920
921      bucket->buckets[i]=NULL;
922      n_Delete(&orig_coef,r);
923      p_Delete(&bucket->coef[i],r);
924      //l1=bucket->buckets_length[i];
925      assume(l1==pLength(p1));
926    }
927    else
928    #endif
929    {
930      //don't do that, pull out gcd
931      #ifdef USE_COEF_BUCKETS
932      if(!(n_IsOne(n,r)))
933      {
934        p1=p_Mult_nn(p1, n, r);
935        n_Delete(&n,r);
936        n=n_Init(1,r);
937      }
938      #endif
939      MULTIPLY_BUCKET(bucket,i);
940      p1 = p_Add_q(p1, bucket->buckets[i],
941                 l1, bucket->buckets_length[i], r);
942      bucket->buckets[i] = NULL;
943      bucket->buckets_length[i] = 0;
944    }
945    i = pLogLength(l1);
946  }
947
948  bucket->buckets[i] = p1;
949#ifdef USE_COEF_BUCKETS
950  assume(bucket->coef[i]==NULL);
951 
952  if (!(n_IsOne(n,r)))
953  {
954    bucket->coef[i]=p_NSet(n,r);
955  }
956  else
957  {
958    bucket->coef[i]=NULL;
959    n_Delete(&n,r);
960  }
961 
962  if ((p1==NULL) && (bucket->coef[i]!=NULL))
963    p_Delete(&bucket->coef[i],r);
964#endif
965  bucket->buckets_length[i]=l1;
966  if (i >= bucket->buckets_used)
967    bucket->buckets_used = i;
968  else
969    kBucketAdjustBucketsUsed(bucket);
970
971  kbTest(bucket);
972}
973
974poly kBucket_ExtractLarger(kBucket_pt bucket, poly q, poly append)
975{
976  if (q == NULL) return append;
977  poly lm;
978  loop
979  {
980    lm = kBucketGetLm(bucket);
981    if (lm == NULL) return append;
982    if (p_LmCmp(lm, q, bucket->bucket_ring) == 1)
983    {
984      lm = kBucketExtractLm(bucket);
985      pNext(append) = lm;
986      pIter(append);
987    }
988    else
989    {
990      return append;
991    }
992  }
993}
994
995/////////////////////////////////////////////////////////////////////////////
996//
997// Extract all monomials from bucket with component comp
998// Return as a polynomial *p with length *l
999// In other words, afterwards
1000// Bpoly = Bpoly - (poly consisting of all monomials with component comp)
1001// and components of monomials of *p are all 0
1002//
1003
1004// Hmm... for now I'm too lazy to implement those independent of currRing
1005// But better declare it extern than including polys.h
1006extern void pTakeOutComp(poly *p, Exponent_t comp, poly *q, int *lq);
1007void pDecrOrdTakeOutComp(poly *p, Exponent_t comp, Order_t order,
1008                         poly *q, int *lq);
1009
1010void kBucketTakeOutComp(kBucket_pt bucket,
1011                        Exponent_t comp,
1012                        poly *r_p, int *l)
1013{
1014  poly p = NULL, q;
1015  int i, lp = 0, lq;
1016
1017#ifndef HAVE_PSEUDO_BUCKETS
1018  kBucketMergeLm(bucket);
1019  for (i=1; i<=bucket->buckets_used; i++)
1020  {
1021    if (bucket->buckets[i] != NULL)
1022    {
1023      MULTIPLY_BUCKET(bucket,i);
1024      pTakeOutComp(&(bucket->buckets[i]), comp, &q, &lq);
1025      if (q != NULL)
1026      {
1027        assume(pLength(q) == lq);
1028        bucket->buckets_length[i] -= lq;
1029        assume(pLength(bucket->buckets[i]) == bucket->buckets_length[i]);
1030        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1031      }
1032    }
1033  }
1034  kBucketAdjustBucketsUsed(bucket);
1035#else
1036  pTakeOutComp(&(bucket->p), comp, &p, &lp);
1037  (bucket->l) -= lp;
1038#endif
1039  *r_p = p;
1040  *l = lp;
1041
1042  kbTest(bucket);
1043}
1044
1045void kBucketDecrOrdTakeOutComp(kBucket_pt bucket,
1046                               Exponent_t comp, Order_t order,
1047                               poly *r_p, int *l)
1048{
1049  poly p = NULL, q;
1050  int i, lp = 0, lq;
1051
1052#ifndef HAVE_PSEUDO_BUCKETS
1053  kBucketMergeLm(bucket);
1054  for (i=1; i<=bucket->buckets_used; i++)
1055  {
1056    if (bucket->buckets[i] != NULL)
1057    {
1058      MULTIPLY_BUCKET(bucket,i);
1059      pDecrOrdTakeOutComp(&(bucket->buckets[i]), comp, order, &q, &lq);
1060      if (q != NULL)
1061      {
1062        bucket->buckets_length[i] -= lq;
1063        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1064      }
1065    }
1066  }
1067  kBucketAdjustBucketsUsed(bucket);
1068#else
1069  pDecrOrdTakeOutComp(&(bucket->p), comp, order, &p, &lp);
1070  (bucket->l) -= lp;
1071#endif
1072
1073  *r_p = p;
1074  *l = lp;
1075}
1076
1077/////////////////////////////////////////////////////////////////////////////
1078// Reduction of Bpoly with a given poly
1079//
1080
1081extern int ksCheckCoeff(number *a, number *b);
1082
1083number kBucketPolyRed(kBucket_pt bucket,
1084                      poly p1, int l1,
1085                      poly spNoether)
1086{
1087  assume((!rIsPluralRing(bucket->bucket_ring))||p_LmEqual(p1,kBucketGetLm(bucket), bucket->bucket_ring)); 
1088  assume(p1 != NULL &&
1089         p_DivisibleBy(p1,  kBucketGetLm(bucket), bucket->bucket_ring));
1090  assume(pLength(p1) == (int) l1);
1091
1092  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
1093  BOOLEAN reset_vec=FALSE;
1094  number rn;
1095
1096  if(a1==NULL)
1097  {
1098    p_DeleteLm(&lm, bucket->bucket_ring);
1099    return nInit(1);
1100  }
1101
1102  if (! nIsOne(pGetCoeff(p1)))
1103  {
1104    number an = pGetCoeff(p1), bn = pGetCoeff(lm);
1105    int ct = ksCheckCoeff(&an, &bn);
1106    p_SetCoeff(lm, bn, bucket->bucket_ring);
1107    if ((ct == 0) || (ct == 2)) kBucket_Mult_n(bucket, an);
1108    rn = an;
1109  }
1110  else
1111  {
1112    rn = nInit(1);
1113  }
1114
1115  if (p_GetComp(p1, bucket->bucket_ring) != p_GetComp(lm, bucket->bucket_ring))
1116  {
1117    p_SetCompP(a1, p_GetComp(lm, bucket->bucket_ring), bucket->bucket_ring);
1118    reset_vec = TRUE;
1119    p_SetComp(lm, p_GetComp(p1, bucket->bucket_ring), bucket->bucket_ring);
1120    p_Setm(lm, bucket->bucket_ring);
1121  }
1122
1123  p_ExpVectorSub(lm, p1, bucket->bucket_ring);
1124  l1--;
1125
1126  assume(l1==pLength(a1));
1127  BOOLEAN backuped=FALSE;
1128  number coef;
1129  if(l1==1) {
1130   
1131    //if (rField_is_Q(bucket->bucket_ring)) {
1132      //avoid this for function fields, as gcds are expensive at the moment
1133     
1134     
1135      coef=p_GetCoeff(a1,bucket->bucket_ring);
1136      lm=p_Mult_nn(lm, coef, bucket->bucket_ring);
1137      p_SetCoeff0(a1, n_Init(1,bucket->bucket_ring), bucket->bucket_ring);
1138      backuped=TRUE;
1139      //WARNING: not thread_safe
1140    //deletes coef as side effect
1141    //}
1142  }
1143
1144 
1145  kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1146 
1147  if (backuped)
1148    p_SetCoeff0(a1,coef,bucket->bucket_ring);
1149  p_DeleteLm(&lm, bucket->bucket_ring);
1150  if (reset_vec) p_SetCompP(a1, 0, bucket->bucket_ring);
1151  kbTest(bucket);
1152  return rn;
1153}
1154static BOOLEAN nIsPseudoUnit(number n, ring r){
1155    if (rField_is_Zp(r))
1156        return TRUE;
1157    if (r->parameter==NULL)
1158    {
1159        if (r->cf->nSize(n)==1)
1160            return TRUE;
1161        else
1162            return FALSE;
1163    }
1164    //if (r->parameter!=NULL)
1165    number one=n_Init(1,r);
1166    if (n_Equal(n,one,r)) {
1167    n_Delete(&one,r);
1168    return TRUE;
1169    }
1170    n_Delete(&one,r);
1171    number minus_one=n_Init(-1,r);
1172    if (n_Equal(n,minus_one,r)){
1173        n_Delete(&minus_one,r);
1174        return TRUE;
1175    }
1176    return FALSE;
1177}
1178
1179void kBucketSimpleContent(kBucket_pt bucket)
1180{
1181  #ifdef USE_COEF_BUCKETS
1182  ring r=bucket->bucket_ring;
1183  int i;
1184  //PrintS("HHHHHHHHHHHHH");
1185  for (i=0;i<=MAX_BUCKET;i++)
1186  {
1187    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1188    //    PrintS("H2H2H2");
1189    if (i==0)
1190    {
1191      assume(bucket->buckets[i]==NULL);
1192    }
1193    if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]==NULL))
1194      return;
1195  }
1196  for (i=0;i<=MAX_BUCKET;i++)
1197  {
1198    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1199    //    PrintS("H2H2H2");
1200    if (i==0)
1201    {
1202      assume(bucket->buckets[i]==NULL);
1203    }
1204    if ((bucket->buckets[i]!=NULL) &&     (nIsPseudoUnit(p_GetCoeff(bucket->coef[i],r),r)))
1205
1206      return;
1207  }
1208  //return;
1209 
1210  number coef=n_Init(0,r);
1211  //ATTENTION: will not work correct for GB over ring
1212  //if (TEST_OPT_PROT)
1213  //    PrintS("CCCCCCCCCCCCC");
1214  for (i=MAX_BUCKET;i>=0;i--)
1215  {
1216    if (i==0)
1217    {
1218      assume(bucket->buckets[i]==NULL);
1219    }
1220    if (bucket->buckets[i]!=NULL)
1221    {
1222      assume(bucket->coef[i]!=NULL);
1223      assume(!(n_IsZero(pGetCoeff(bucket->coef[i]),r)));
1224     
1225     
1226      //in this way it should crash on programming errors, yeah
1227      number temp=nGcd(coef, pGetCoeff(bucket->coef[i]),r);
1228      n_Delete(&coef,r );
1229      coef=temp;
1230      if (nIsPseudoUnit(coef,r))
1231      {
1232        n_Delete(&coef,r);
1233        return;
1234      }
1235      assume(!(n_IsZero(coef,r)));
1236    }
1237  }
1238  if (n_IsZero(coef,r)){
1239    n_Delete(&coef,r);
1240    return;
1241    }
1242  if (TEST_OPT_PROT)
1243    PrintS("S");
1244  for(i=0;i<=MAX_BUCKET;i++)
1245  {
1246   
1247    if (bucket->buckets[i]!=NULL)
1248    {
1249      assume(!(n_IsZero(coef,r)));
1250      assume(bucket->coef[i]!=NULL);
1251      number lc=p_GetCoeff(bucket->coef[i],r);
1252      p_SetCoeff(bucket->coef[i], n_IntDiv(lc,coef,r),r);
1253      assume(!(n_IsZero(p_GetCoeff(bucket->coef[i],r),r)));
1254    }
1255  }
1256  n_Delete(&coef,r);
1257  #endif
1258}
1259
1260
1261poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i)
1262{
1263  assume(bucket->buckets[i]!=NULL);
1264
1265  ring r=bucket->bucket_ring;
1266  poly p=bucket->buckets[i];
1267  bucket->buckets_length[i]--;
1268#ifdef USE_COEF_BUCKETS
1269  if (bucket->coef[i]!=NULL)
1270  {
1271    poly next=pNext(p);
1272    if (next==NULL)
1273    {
1274      MULTIPLY_BUCKET(bucket,i);
1275      p=bucket->buckets[i];
1276      bucket->buckets[i]=NULL;
1277      return p;
1278    }
1279    else
1280    {
1281      bucket->buckets[i]=next;
1282      number c=p_GetCoeff(bucket->coef[i],r);
1283      pNext(p)=NULL;
1284      p=p_Mult_nn(p,c,r);
1285      assume(p!=NULL);
1286      return p;
1287    }
1288  }
1289  else
1290#endif
1291  {
1292    bucket->buckets[i]=pNext(bucket->buckets[i]);
1293    pNext(p)=NULL;
1294    assume(p!=NULL);
1295    return p;
1296  }
1297}
Note: See TracBrowser for help on using the repository browser.