source: git/Singular/hutil.cc @ 18dd47

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