source: git/kernel/hutil.cc @ fbc7cb

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