source: git/kernel/kbuckets.cc @ 086a96

fieker-DuValspielwiese
Last change on this file since 086a96 was 086a96, checked in by Hans Schoenemann <hannes@…>, 13 years ago
simplify nIsPseudoUnit git-svn-id: file:///usr/local/Singular/svn/trunk@14034 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 31.3 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5
6#include <kernel/mod2.h>
7#include <kernel/structs.h>
8#include <omalloc/omalloc.h>
9#include <kernel/p_polys.h>
10#include <kernel/febase.h>
11#include <kernel/pShallowCopyDelete.h>
12#include <kernel/numbers.h>
13#include <kernel/ring.h>
14#include <kernel/p_Procs.h>
15#include <kernel/kutil.h>
16#include <kernel/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],B->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)) B->buckets_length[i] = pLength(B->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);
1016
1017void kBucketTakeOutComp(kBucket_pt bucket,
1018                        long comp,
1019                        poly *r_p, int *l)
1020{
1021  poly p = NULL, q;
1022  int i, lp = 0, lq;
1023
1024#ifndef HAVE_PSEUDO_BUCKETS
1025  kBucketMergeLm(bucket);
1026  for (i=1; i<=bucket->buckets_used; i++)
1027  {
1028    if (bucket->buckets[i] != NULL)
1029    {
1030      MULTIPLY_BUCKET(bucket,i);
1031      pTakeOutComp(&(bucket->buckets[i]), comp, &q, &lq);
1032      if (q != NULL)
1033      {
1034        assume(pLength(q) == lq);
1035        bucket->buckets_length[i] -= lq;
1036        assume(pLength(bucket->buckets[i]) == bucket->buckets_length[i]);
1037        p = p_Add_q(p, q, lp, lq, bucket->bucket_ring);
1038      }
1039    }
1040  }
1041  kBucketAdjustBucketsUsed(bucket);
1042#else
1043  pTakeOutComp(&(bucket->p), comp, &p, &lp);
1044  (bucket->l) -= lp;
1045#endif
1046  *r_p = p;
1047  *l = lp;
1048
1049  kbTest(bucket);
1050}
1051
1052/////////////////////////////////////////////////////////////////////////////
1053// Reduction of Bpoly with a given poly
1054//
1055
1056extern int ksCheckCoeff(number *a, number *b);
1057
1058number kBucketPolyRed(kBucket_pt bucket,
1059                      poly p1, int l1,
1060                      poly spNoether)
1061{
1062  assume((!rIsPluralRing(bucket->bucket_ring))||p_LmEqual(p1,kBucketGetLm(bucket), bucket->bucket_ring));
1063  assume(p1 != NULL &&
1064         p_DivisibleBy(p1,  kBucketGetLm(bucket), bucket->bucket_ring));
1065  assume(pLength(p1) == (int) l1);
1066
1067  poly a1 = pNext(p1), lm = kBucketExtractLm(bucket);
1068  BOOLEAN reset_vec=FALSE;
1069  number rn;
1070
1071  /* we shall reduce bucket=bn*lm+... by p1=an*t+a1 where t=lm(p1)
1072     and an,bn shall be defined further down only if lc(p1)!=1
1073     we already know: an|bn and t|lm */
1074  if(a1==NULL)
1075  {
1076    p_LmDelete(&lm, bucket->bucket_ring);
1077    return nInit(1);
1078  }
1079
1080  if (! nIsOne(pGetCoeff(p1)))
1081  {
1082    number an = pGetCoeff(p1), bn = pGetCoeff(lm);
1083//StringSetS("##### an = "); nWrite(an); PrintS(StringAppend("\n"));
1084//StringSetS("##### bn = "); nWrite(bn); PrintS(StringAppend("\n"));
1085    /* ksCheckCoeff: divide out gcd from an and bn: */
1086    int ct = ksCheckCoeff(&an, &bn);
1087    /* the previous command returns ct=0 or ct=2 iff an!=1
1088       note: an is now 1 or -1 */
1089
1090    /* setup factor for p1 which cancels leading terms */
1091    p_SetCoeff(lm, bn, bucket->bucket_ring);
1092    if ((ct == 0) || (ct == 2))
1093    {
1094      /* next line used to be here before but is WRONG:
1095      kBucket_Mult_n(bucket, an);
1096        its use would result in a wrong sign for the tail of bucket
1097        in the reduction */
1098
1099      /* correct factor for cancelation by changing sign if an=-1 */
1100      if (rField_is_Ring())
1101        lm = p_Mult_nn(lm, an, bucket->bucket_ring);
1102      else
1103        kBucket_Mult_n(bucket, an);
1104    }
1105    rn = an;
1106  }
1107  else
1108  {
1109    rn = nInit(1);
1110  }
1111
1112  if (p_GetComp(p1, bucket->bucket_ring) != p_GetComp(lm, bucket->bucket_ring))
1113  {
1114    p_SetCompP(a1, p_GetComp(lm, bucket->bucket_ring), bucket->bucket_ring);
1115    reset_vec = TRUE;
1116    p_SetComp(lm, p_GetComp(p1, bucket->bucket_ring), bucket->bucket_ring);
1117    p_Setm(lm, bucket->bucket_ring);
1118  }
1119
1120  p_ExpVectorSub(lm, p1, bucket->bucket_ring);
1121  l1--;
1122
1123  assume(l1==pLength(a1));
1124  BOOLEAN backuped=FALSE;
1125  number coef;
1126  #if 0
1127  //@Viktor, don't ignore coefficients on monomials
1128  if(l1==1) {
1129
1130    //if (rField_is_Q(bucket->bucket_ring)) {
1131      //avoid this for function fields, as gcds are expensive at the moment
1132
1133
1134      coef=p_GetCoeff(a1,bucket->bucket_ring);
1135      lm=p_Mult_nn(lm, coef, bucket->bucket_ring);
1136      p_SetCoeff0(a1, n_Init(1,bucket->bucket_ring), bucket->bucket_ring);
1137      backuped=TRUE;
1138      //WARNING: not thread_safe
1139    //deletes coef as side effect
1140    //}
1141  }
1142  #endif
1143
1144  kBucket_Minus_m_Mult_p(bucket, lm, a1, &l1, spNoether);
1145
1146  if (backuped)
1147    p_SetCoeff0(a1,coef,bucket->bucket_ring);
1148  p_LmDelete(&lm, bucket->bucket_ring);
1149  if (reset_vec) p_SetCompP(a1, 0, bucket->bucket_ring);
1150  kbTest(bucket);
1151  return rn;
1152}
1153static BOOLEAN nIsPseudoUnit(number n, ring r)
1154{
1155  if (rField_is_Zp(r))
1156    return TRUE;
1157  if (r->parameter==NULL)
1158  {
1159    return (r->cf->nSize(n)==1);
1160  }
1161  //if (r->parameter!=NULL)
1162  return ((n_IsOne(n,r)) || (n_IsMOne(n,r)));
1163}
1164
1165void kBucketSimpleContent(kBucket_pt bucket)
1166{
1167  #ifdef USE_COEF_BUCKETS
1168  ring r=bucket->bucket_ring;
1169  int i;
1170  //PrintS("HHHHHHHHHHHHH");
1171  for (i=0;i<=MAX_BUCKET;i++)
1172  {
1173    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1174    //    PrintS("H2H2H2");
1175    if (i==0)
1176    {
1177      assume(bucket->buckets[i]==NULL);
1178    }
1179    if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]==NULL))
1180      return;
1181  }
1182  for (i=0;i<=MAX_BUCKET;i++)
1183  {
1184    //if ((bucket->buckets[i]!=NULL) && (bucket->coef[i]!=NULL))
1185    //    PrintS("H2H2H2");
1186    if (i==0)
1187    {
1188      assume(bucket->buckets[i]==NULL);
1189    }
1190    if ((bucket->buckets[i]!=NULL)
1191    && (nIsPseudoUnit(p_GetCoeff(bucket->coef[i],r),r)))
1192      return;
1193  }
1194  //return;
1195
1196  number coef=n_Init(0,r);
1197  //ATTENTION: will not work correct for GB over ring
1198  //if (TEST_OPT_PROT)
1199  //    PrintS("CCCCCCCCCCCCC");
1200  for (i=MAX_BUCKET;i>=0;i--)
1201  {
1202    if (i==0)
1203    {
1204      assume(bucket->buckets[i]==NULL);
1205    }
1206    if (bucket->buckets[i]!=NULL)
1207    {
1208      assume(bucket->coef[i]!=NULL);
1209      assume(!(n_IsZero(pGetCoeff(bucket->coef[i]),r)));
1210
1211      //in this way it should crash on programming errors, yeah
1212      number temp=nGcd(coef, pGetCoeff(bucket->coef[i]),r);
1213      n_Delete(&coef,r );
1214      coef=temp;
1215      if (nIsPseudoUnit(coef,r))
1216      {
1217        n_Delete(&coef,r);
1218        return;
1219      }
1220      assume(!(n_IsZero(coef,r)));
1221    }
1222  }
1223  if (n_IsZero(coef,r))
1224  {
1225    n_Delete(&coef,r);
1226    return;
1227  }
1228  if (TEST_OPT_PROT)
1229    PrintS("S");
1230  for(i=0;i<=MAX_BUCKET;i++)
1231  {
1232    if (bucket->buckets[i]!=NULL)
1233    {
1234      assume(!(n_IsZero(coef,r)));
1235      assume(bucket->coef[i]!=NULL);
1236      number lc=p_GetCoeff(bucket->coef[i],r);
1237      p_SetCoeff(bucket->coef[i], n_IntDiv(lc,coef,r),r);
1238      assume(!(n_IsZero(p_GetCoeff(bucket->coef[i],r),r)));
1239    }
1240  }
1241  n_Delete(&coef,r);
1242  #endif
1243}
1244
1245
1246poly kBucketExtractLmOfBucket(kBucket_pt bucket, int i)
1247{
1248  assume(bucket->buckets[i]!=NULL);
1249
1250  ring r=bucket->bucket_ring;
1251  poly p=bucket->buckets[i];
1252  bucket->buckets_length[i]--;
1253#ifdef USE_COEF_BUCKETS
1254  if (bucket->coef[i]!=NULL)
1255  {
1256    poly next=pNext(p);
1257    if (next==NULL)
1258    {
1259      MULTIPLY_BUCKET(bucket,i);
1260      p=bucket->buckets[i];
1261      bucket->buckets[i]=NULL;
1262      return p;
1263    }
1264    else
1265    {
1266      bucket->buckets[i]=next;
1267      number c=p_GetCoeff(bucket->coef[i],r);
1268      pNext(p)=NULL;
1269      p=p_Mult_nn(p,c,r);
1270      assume(p!=NULL);
1271      return p;
1272    }
1273  }
1274  else
1275#endif
1276  {
1277    bucket->buckets[i]=pNext(bucket->buckets[i]);
1278    pNext(p)=NULL;
1279    assume(p!=NULL);
1280    return p;
1281  }
1282}
Note: See TracBrowser for help on using the repository browser.