source: git/kernel/kbuckets.cc @ be5fcb

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