source: git/kernel/hutil.cc @ cd6ff49

spielwiese
Last change on this file since cd6ff49 was cd6ff49, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: fixed typo git-svn-id: file:///usr/local/Singular/svn/trunk@8350 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 16.3 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: hutil.cc,v 1.8 2005-06-09 08:30:02 Singular Exp $ */
5/*
6* ABSTRACT: Utilities for staircase operations
7*/
8
9#include "mod2.h"
10#include "structs.h"
11#include "febase.h"
12#include "omalloc.h"
13#include "ideals.h"
14#include "polys.h"
15#include "hutil.h"
16
17scfmon hexist, hstc, hrad, hwork;
18scmon hpure, hpur0;
19varset hvar, hsel;
20int  hNexist, hNstc, hNrad, hNvar, hNpure;
21int hisModule;
22monf stcmem, radmem;
23
24// Making a global "security" copy of the allocated exponent vectors
25// is a dirty fix for correct memory disallocation: It would be
26// better, if either the fields of heist are never touched
27// (i.e. changed) except in hInit, or, if hInit would return the
28// "security" copy as well. But then, all the relevant data is held in
29// global variables, so we might do that here, as well.
30static scfmon hsecure= NULL;
31
32scfmon hInit(ideal S, ideal Q, int *Nexist, ring tailRing)
33{
34  int  sl, ql, i, k = 0;
35  polyset si, qi, ss;
36  scfmon ex, ek;
37  if (tailRing != currRing)
38    hisModule = idRankFreeModule(S, currRing, tailRing);
39  else
40    hisModule = idRankFreeModule(S);
41  if (hisModule < 0)
42    hisModule = 0;
43  if (S!=NULL)
44  {
45    si = S->m;
46    sl = IDELEMS(S);
47  }
48  else
49  {
50    si = NULL;
51    sl = 0;
52  }
53  if (Q!=NULL)
54  {
55    qi = Q->m;
56    ql = IDELEMS(Q);
57  }
58  else
59  {
60    qi = NULL;
61    ql = 0;
62  }
63  if ((sl + ql) == 0)
64  {
65    *Nexist = 0;
66    return NULL;
67  }
68  ss = si;
69  for (i = sl; i>0; i--)
70  {
71    if (*ss!=0)
72      k++;
73    ss++;
74  }
75  ss = qi;
76  for (i = ql; i>0; i--)
77  {
78    if (*ss!=0)
79      k++;
80    ss++;
81  }
82  *Nexist = k;
83  if (k==0)
84    return NULL;
85  ek = ex = (scfmon)omAlloc0(k * sizeof(scmon));
86  hsecure = (scfmon) omAlloc0(k * sizeof(scmon));
87  for (i = sl; i>0; i--)
88  {
89    if (*si!=NULL)
90    {
91      *ek = (scmon) omAlloc((pVariables+1)*sizeof(int));
92      pGetExpV(*si, *ek);
93      ek++;
94    }
95    si++;
96  }
97  for (i = ql; i>0; i--)
98  {
99    if (*qi!=NULL)
100    {
101      *ek = (scmon) omAlloc((pVariables+1)*sizeof(int));
102      pGetExpV(*qi, *ek);
103      ek++;
104    }
105    qi++;
106  }
107  memcpy(hsecure, ex, k * sizeof(scmon));
108  return ex;
109}
110
111void hWeight()
112{
113  int i, k;
114  int x;
115
116  i = pVariables;
117  loop
118  {
119    if (pWeight(i) != 1) break;
120    i--;
121    if (i == 0) return;
122  }
123  for (i=pVariables; i>0; i--)
124  {
125    x = pWeight(i);
126    if (x != 1)
127    {
128      for (k=hNexist-1; k>=0; k--)
129      {
130        hexist[k][i] *= x;
131      }
132    }
133  }
134}
135
136void hDelete(scfmon ev, int ev_length)
137{
138  int i;
139
140  if (ev_length>0)
141  {
142    for (i=ev_length-1;i>=0;i--)
143      omFreeSize(hsecure[i],(pVariables+1)*sizeof(int));
144    omFreeSize(hsecure, ev_length*sizeof(scmon));
145    omFreeSize(ev,  ev_length*sizeof(scmon));
146  }
147}
148
149
150void hComp(scfmon exist, int Nexist, int ak, scfmon stc, int *Nstc)
151{
152  int k = 0;
153  scfmon ex = exist, co = stc;
154  int i;
155
156  for (i = Nexist; i>0; i--)
157  {
158    if (((**ex) == 0) || ((**ex) == ak))
159    {
160      *co = *ex;
161      co++;
162      k++;
163    }
164    ex++;
165  }
166  *Nstc = k;
167}
168
169
170void hSupp(scfmon stc, int Nstc, varset var, int *Nvar)
171{
172  int  nv, i0, i1, i, j;
173  nv = i0 = *Nvar;
174  i1 = 0;
175  for (i = 1; i <= nv; i++)
176  {
177    j = 0;
178    loop
179    {
180      if (stc[j][i]>0)
181      {
182        i1++;
183        var[i1] = i;
184        break;
185      }
186      j++;
187      if (j == Nstc)
188      {
189        var[i0] = i;
190        i0--;
191        break;
192      }
193    }
194  }
195  *Nvar = i1;
196}
197
198void hOrdSupp(scfmon stc, int Nstc, varset var, int Nvar)
199{
200  int  i, i1, j, jj, k, l;
201  int  x;
202  scmon temp, count;
203  float o, h, g, *v1;
204
205  v1 = (float *)omAlloc(Nvar * sizeof(float));
206  temp = (int *)omAlloc(Nstc * sizeof(int));
207  count = (int *)omAlloc(Nstc * sizeof(int));
208  for (i = 1; i <= Nvar; i++)
209  {
210    i1 = var[i];
211    *temp = stc[0][i1];
212    *count = 1;
213    jj = 1;
214    for (j = 1; j < Nstc; j++)
215    {
216      x = stc[j][i1];
217      k = 0;
218      loop
219      {
220        if (x > temp[k])
221        {
222          k++;
223          if (k == jj)
224          {
225            temp[k] = x;
226            count[k] = 1;
227            jj++;
228            break;
229          }
230        }
231        else if (x < temp[k])
232        {
233          for (l = jj; l > k; l--)
234          {
235            temp[l] = temp[l-1];
236            count[l] = count[l-1];
237          }
238          temp[k] = x;
239          count[k] = 1;
240          jj++;
241          break;
242        }
243        else
244        {
245          count[k]++;
246          break;
247        }
248      }
249    }
250    h = 0.0;
251    o = (float)Nstc/(float)jj;
252    for(j = 0; j < jj; j++)
253    {
254       g = (float)count[j];
255       if (g > o)
256         g -= o;
257       else
258         g = o - g;
259       if (g > h)
260         h = g;
261    }
262    v1[i-1] = h * (float)jj;
263  }
264  omFreeSize((ADDRESS)count, Nstc * sizeof(int));
265  omFreeSize((ADDRESS)temp, Nstc * sizeof(int));
266  for (i = 1; i < Nvar; i++)
267  {
268    i1 = var[i+1];
269    h = v1[i];
270    j = 0;
271    loop
272    {
273      if (h > v1[j])
274      {
275        for (l = i; l > j; l--)
276        {
277          v1[l] = v1[l-1];
278          var[l+1] = var[l];
279        }
280        v1[j] = h;
281        var[j+1] = i1;
282        break;
283      }
284      j++;
285      if (j == i)
286        break;
287    }
288  }
289  omFreeSize((ADDRESS)v1, Nvar * sizeof(float));
290}
291
292
293static void hShrink(scfmon co, int a, int Nco)
294{
295  while ((co[a]!=NULL) && (a<Nco)) a++;
296  int i = a;
297  int j;
298  for (j = a; j < Nco; j++)
299  {
300    if (co[j]!=NULL)
301    {
302      co[i] = co[j];
303      i++;
304    }
305  }
306}
307
308
309void hStaircase(scfmon stc, int *Nstc, varset var, int Nvar)
310{
311  int  nc = *Nstc;
312  if (nc < 2)
313    return;
314  int z = 0;
315  int i = 0;
316  int j = 1;
317  scmon n = stc[1 /*j*/];
318  scmon o = stc[0];
319  int k = Nvar;
320  loop
321  {
322    int k1 = var[k];
323    if (o[k1] > n[k1])
324    {
325      loop
326      {
327        k--;
328        if (k==0)
329        {
330          stc[i] = NULL;
331          z++;
332          break;
333        }
334        else
335        {
336          k1 = var[k];
337          if (o[k1] < n[k1])
338            break;
339        }
340      }
341      k = Nvar;
342    }
343    else if (o[k1] < n[k1])
344    {
345      loop
346      {
347        k--;
348        if (k==0)
349        {
350          stc[j] = NULL;
351          z++;
352          break;
353        }
354        else
355        {
356          k1 = var[k];
357          if (o[k1] > n[k1])
358            break;
359        }
360      }
361      k = Nvar;
362    }
363    else
364    {
365      k--;
366      if (k==0)
367      {
368        stc[j] = NULL;
369        z++;
370        k = Nvar;
371      }
372    }
373    if (k == Nvar)
374    {
375      if (stc[j]==NULL)
376        i = j - 1;
377      loop
378      {
379        i++;
380        if (i == j)
381        {
382          i = -1;
383          j++;
384          if (j < nc)
385            n = stc[j];
386          else
387          {
388            if (z!=0)
389            {
390              *Nstc -= z;
391              hShrink(stc, 0, nc);
392            }
393            return;
394          }
395        }
396        else if (stc[i]!=NULL)
397        {
398          o = stc[i];
399          break;
400        }
401      }
402    }
403  }
404}
405
406
407void hRadical(scfmon rad, int *Nrad, int Nvar)
408{
409  int  nc = *Nrad, z = 0, i, j, k;
410  scmon n, o;
411  if (nc < 2)
412    return;
413  i = 0;
414  j = 1;
415  n = rad[j];
416  o = rad[0];
417  k = Nvar;
418  loop
419  {
420    if ((o[k]!=0) && (n[k]==0))
421    {
422      loop
423      {
424        k--;
425        if (k==0)
426        {
427          rad[i] = NULL;
428          z++;
429          break;
430        }
431        else
432        {
433          if ((o[k]==0) && (n[k]!=0))
434            break;
435        }
436      }
437      k = Nvar;
438    }
439    else if (!o[k] && n[k])
440    {
441      loop
442      {
443        k--;
444        if (!k)
445        {
446          rad[j] = NULL;
447          z++;
448          break;
449        }
450        else
451        {
452          if (o[k] && !n[k])
453            break;
454        }
455      }
456      k = Nvar;
457    }
458    else
459    {
460      k--;
461      if (!k)
462      {
463        rad[j] = NULL;
464        z++;
465        k = Nvar;
466      }
467    }
468    if (k == Nvar)
469    {
470      if (!rad[j])
471        i = j - 1;
472      loop
473      {
474        i++;
475        if (i == j)
476        {
477          i = -1;
478          j++;
479          if (j < nc)
480            n = rad[j];
481          else
482          {
483            if (z)
484            {
485              *Nrad -= z;
486              hShrink(rad, 0, nc);
487            }
488            return;
489          }
490        }
491        else if (rad[i])
492        {
493          o = rad[i];
494          break;
495        }
496      }
497    }
498  }
499}
500
501
502void hLexS(scfmon stc, int Nstc, varset var, int Nvar)
503{
504  if (Nstc < 2)
505    return;
506  int  j = 1, i = 0;
507  scmon n = stc[j];
508  scmon o = stc[0];
509  int k = Nvar;
510  loop
511  {
512    int k1 = var[k];
513    if (o[k1] < n[k1])
514    {
515      i++;
516      if (i < j)
517      {
518        o = stc[i];
519        k = Nvar;
520      }
521      else
522      {
523        j++;
524        if (j < Nstc)
525        {
526          i = 0;
527          o = stc[0];
528          n = stc[j];
529          k = Nvar;
530        }
531        else
532          return;
533      }
534    }
535    else if (o[k1] > n[k1])
536    {
537      int tmp_k;
538      for (tmp_k = j; tmp_k > i; tmp_k--)
539        stc[tmp_k] = stc[tmp_k - 1];
540      stc[i] = n;
541      j++;
542      if (j < Nstc)
543      {
544        i = 0;
545        o = stc[0];
546        n = stc[j];
547        k = Nvar;
548      }
549      else
550        return;
551    }
552    else
553    {
554      k--;
555      if (k<=0) return;
556    }
557  }
558}
559
560
561void hLexR(scfmon rad, int Nrad, varset var, int Nvar)
562{
563  int  j = 1, i = 0, k, k1;
564  scmon n, o;
565  if (Nrad < 2)
566    return;
567  n = rad[j];
568  o = rad[0];
569  k = Nvar;
570  loop
571  {
572    k1 = var[k];
573    if (!o[k1] && n[k1])
574    {
575      i++;
576      if (i < j)
577      {
578        o = rad[i];
579        k = Nvar;
580      }
581      else
582      {
583        j++;
584        if (j < Nrad)
585        {
586          i = 0;
587          o = rad[0];
588          n = rad[j];
589          k = Nvar;
590        }
591        else
592          return;
593      }
594    }
595    else if (o[k1] && !n[k1])
596    {
597      for (k = j; k > i; k--)
598        rad[k] = rad[k - 1];
599      rad[i] = n;
600      j++;
601      if (j < Nrad)
602      {
603        i = 0;
604        o = rad[0];
605        n = rad[j];
606        k = Nvar;
607      }
608      else
609        return;
610    }
611    else
612      k--;
613  }
614}
615
616
617void hPure(scfmon stc, int a, int *Nstc, varset var, int Nvar,
618 scmon pure, int *Npure)
619{
620  int  nc = *Nstc, np = 0, nq = 0, j, i, i1, c, l;
621  scmon x;
622  for (j = a; j < nc; j++)
623  {
624    x = stc[j];
625    i = Nvar;
626    c = 2;
627    l = 0;
628    loop
629    {
630      i1 = var[i];
631      if (x[i1])
632      {
633        c--;
634        if (!c)
635        {
636          l = 0;
637          break;
638        }
639        else if (c == 1)
640          l = i1;
641      }
642      i--;
643      if (!i)
644        break;
645    }
646    if (l)
647    {
648      if (!pure[l])
649      {
650        np++;
651        pure[l] = x[l];
652      }
653      else if (x[l] < pure[l])
654        pure[l] = x[l];
655      stc[j] = NULL;
656      nq++;
657    }
658  }
659  *Npure = np;
660  if (nq!=0)
661  {
662    *Nstc -= nq;
663    hShrink(stc, a, nc);
664  }
665}
666
667
668void hElimS(scfmon stc, int *e1, int a2, int e2, varset var, int Nvar)
669{
670  int  nc = *e1, z = 0, i, j, k, k1;
671  scmon n, o;
672  if (!nc || (a2 == e2))
673    return;
674  j = 0;
675  i = a2;
676  o = stc[i];
677  n = stc[0];
678  k = Nvar;
679  loop
680  {
681    k1 = var[k];
682    if (o[k1] > n[k1])
683    {
684      k = Nvar;
685      i++;
686      if (i < e2)
687        o = stc[i];
688      else
689      {
690        j++;
691        if (j < nc)
692        {
693          i = a2;
694          o = stc[i];
695          n = stc[j];
696        }
697        else
698        {
699          if (z!=0)
700          {
701            *e1 -= z;
702            hShrink(stc, 0, nc);
703          }
704          return;
705        }
706      }
707    }
708    else
709    {
710      k--;
711      if (k==0)
712      {
713        stc[j] = NULL;
714        z++;
715        j++;
716        if (j < nc)
717        {
718          i = a2;
719          o = stc[i];
720          n = stc[j];
721          k = Nvar;
722        }
723        else
724        {
725          if (z!=0)
726          {
727            *e1 -= z;
728            hShrink(stc, 0, nc);
729          }
730          return;
731        }
732      }
733    }
734  }
735}
736
737
738void hElimR(scfmon rad, int *e1, int a2, int e2, varset var, int Nvar)
739{
740  int  nc = *e1, z = 0, i, j, k, k1;
741  scmon n, o;
742  if (!nc || (a2 == e2))
743    return;
744  j = 0;
745  i = a2;
746  o = rad[i];
747  n = rad[0];
748  k = Nvar;
749  loop
750  {
751    k1 = var[k];
752    if (o[k1] && !n[k1])
753    {
754      k = Nvar;
755      i++;
756      if (i < e2)
757        o = rad[i];
758      else
759      {
760        j++;
761        if (j < nc)
762        {
763          i = a2;
764          o = rad[i];
765          n = rad[j];
766        }
767        else
768        {
769          if (z!=0)
770          {
771            *e1 -= z;
772            hShrink(rad, 0, nc);
773          }
774          return;
775        }
776      }
777    }
778    else
779    {
780      k--;
781      if (!k)
782      {
783        rad[j] = NULL;
784        z++;
785        j++;
786        if (j < nc)
787        {
788          i = a2;
789          o = rad[i];
790          n = rad[j];
791          k = Nvar;
792        }
793        else
794        {
795          if (z!=0)
796          {
797            *e1 -= z;
798            hShrink(rad, 0, nc);
799          }
800          return;
801        }
802      }
803    }
804  }
805}
806
807
808void hLex2S(scfmon rad, int e1, int a2, int e2, varset var,
809 int Nvar, scfmon w)
810{
811  int  j0 = 0, j = 0, i = a2, k, k1;
812  scmon n, o;
813  if (!e1)
814  {
815    for (; i < e2; i++)
816      rad[i - a2] = rad[i];
817    return;
818  }  else if (i == e2)
819    return;
820  n = rad[j];
821  o = rad[i];
822  loop
823  {
824    k = Nvar;
825    loop
826    {
827      k1 = var[k];
828      if (o[k1] < n[k1])
829      {
830        w[j0] = o;
831        j0++;
832        i++;
833        if (i < e2)
834        {
835          o = rad[i];
836          break;
837        }
838        else
839        {
840          for (; j < e1; j++)
841          {
842            w[j0] = rad[j];
843            j0++;
844          }
845          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
846          return;
847        }
848      }
849      else if (o[k1] > n[k1])
850      {
851        w[j0] = n;
852        j0++;
853        j++;
854        if (j < e1)
855        {
856          n = rad[j];
857          break;
858        }
859        else
860        {
861          for (; i < e2; i++)
862          {
863            w[j0] = rad[i];
864            j0++;
865          }
866          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
867          return;
868        }
869      }
870      k--;
871    }
872  }
873}
874
875
876void hLex2R(scfmon rad, int e1, int a2, int e2, varset var,
877 int Nvar, scfmon w)
878{
879  int  j0 = 0, j = 0, i = a2, k, k1;
880  scmon n, o;
881  if (!e1)
882  {
883    for (; i < e2; i++)
884      rad[i - a2] = rad[i];
885    return;
886  }
887  else if (i == e2)
888    return;
889  n = rad[j];
890  o = rad[i];
891  loop
892  {
893    k = Nvar;
894    loop
895    {
896      k1 = var[k];
897      if (!o[k1] && n[k1])
898      {
899        w[j0] = o;
900        j0++;
901        i++;
902        if (i < e2)
903        {
904          o = rad[i];
905          break;
906        }
907        else
908        {
909          for (; j < e1; j++)
910          {
911            w[j0] = rad[j];
912            j0++;
913          }
914          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
915          return;
916        }
917      }
918      else if (o[k1] && !n[k1])
919      {
920        w[j0] = n;
921        j0++;
922        j++;
923        if (j < e1)
924        {
925          n = rad[j];
926          break;
927        }
928        else
929        {
930          for (; i < e2; i++)
931          {
932            w[j0] = rad[i];
933            j0++;
934          }
935          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
936          return;
937        }
938      }
939      k--;
940    }
941  }
942}
943
944
945void hStepS(scfmon stc, int Nstc, varset var, int Nvar, int *a, int *x)
946{
947  int  k1, i;
948  int  y;
949  k1 = var[Nvar];
950  y = *x;
951  i = *a;
952  loop
953  {
954    if (y < stc[i][k1])
955    {
956      *a = i;
957      *x = stc[i][k1];
958      return;
959    }
960    i++;
961    if (i == Nstc)
962    {
963      *a = i;
964      return;
965    }
966  }
967}
968
969
970void hStepR(scfmon rad, int Nrad, varset var, int Nvar, int *a)
971{
972  int  k1, i;
973  k1 = var[Nvar];
974  i = 0;
975  loop
976  {
977    if (rad[i][k1])
978    {
979      *a = i;
980      return;
981    }
982    i++;
983    if (i == Nrad)
984    {
985      *a = i;
986      return;
987    }
988  }
989}
990
991
992monf hCreate(int Nvar)
993{
994  monf xmem;
995  int  i;
996  xmem = (monf)omAlloc((Nvar + 1) * sizeof(monp));
997  for (i = Nvar; i>0; i--)
998  {
999    xmem[i] = (monp)omAlloc(LEN_MON);
1000    xmem[i]->mo = NULL;
1001  }
1002  return xmem;
1003}
1004
1005
1006void hKill(monf xmem, int Nvar)
1007{
1008  int  i;
1009  for (i = Nvar; i!=0; i--)
1010  {
1011    if (xmem[i]->mo!=NULL)
1012      omFreeSize((ADDRESS)xmem[i]->mo, xmem[i]->a * sizeof(scmon));
1013    omFreeSize((ADDRESS)xmem[i], LEN_MON);
1014  }
1015  omFreeSize((ADDRESS)xmem, (Nvar + 1) * sizeof(monp));
1016}
1017
1018
1019scfmon hGetmem(int lm, scfmon old, monp monmem)
1020{
1021  scfmon x = monmem->mo;
1022  int  lx = monmem->a;
1023  if ((x==NULL) || (lm > lx))
1024  {
1025    if ((x!=NULL)&&(lx>0)) omFreeSize((ADDRESS)x, lx * sizeof(scmon));
1026    monmem->mo = x = (scfmon)omAlloc(lm * sizeof(scmon));
1027    monmem->a = lm;
1028  }
1029  memcpy(x, old, lm * sizeof(scmon));
1030  return x;
1031}
1032
1033/*
1034* a bug in Metrowerks with "lifetime analysis"
1035*scmon hGetpure(scmon p)
1036*{
1037*  scmon p1, pn;
1038*  p1 = p + 1;
1039*  pn = p1 + pVariables;
1040*  memcpy(pn, p1, pVariables * sizeof(int));
1041*  return pn - 1;
1042*}
1043*/
1044scmon hGetpure(scmon p)
1045{
1046  scmon p1 = p;
1047  scmon pn;
1048  p1++;
1049  pn = p1;
1050  pn += pVariables;
1051  memcpy(pn, p1, pVariables * sizeof(int));
1052  return pn - 1;
1053}
1054
Note: See TracBrowser for help on using the repository browser.