source: git/kernel/kbuckets.cc @ 1b38b23

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