source: git/Singular/hutil.cc @ 4a81ec

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