source: git/kernel/hutil.cc @ 5a9e7b

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