source: git/Singular/LIB/cimonom.lib @ 6fe9a5

spielwiese
Last change on this file since 6fe9a5 was 27d4bb, checked in by Hans Schoenemann <hannes@…>, 13 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13536 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 14.5 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version = "$Id$";
3category="Commutative Algebra";
4info="
5LIBRARY: cimonom.lib   Determines if the toric ideal of an affine monomial curve is a complete intersection
6
7AUTHORS: I.Bermejo,               ibermejo@ull.es
8@*       I.Garcia-Marco,          iggarcia@ull.es
9@*       J.-J.Salazar-Gonzalez,   jjsalaza@ull.es
10
11OVERVIEW:
12  A library for determining if the toric ideal of an affine monomial curve is a complete intersection with NO
13  NEED of computing explicitly a system of generators of such ideal. It also contains procedures to obtain the
14  minimum positive multiple of an integer which is in a semigroup of positive integers.
15  The procedures are based on a paper by Isabel Bermejo, Ignacio Garcia and Juan Jose Salazar-Gonzalez: 'An
16  algorithm to check whether the toric ideal of an affine monomial curve is a complete intersection', Preprint.
17
18SEE ALSO: Integer programming
19
20PROCEDURES:
21 BelongSemig(n,v[,sup]);    checks whether n is in the semigroup generated by v;
22 MinMult(a,b);              computes k, the minimum positive integer such that k*a is in the semigroup of
23                            positive integers generated by the elements in b.
24 CompInt(d);                checks wether I(d) is a complete intersection or not.
25";
26
27LIB "general.lib";
28
29///////////////////////////////////////////////////////////////////////////////////////////////////////////
30//
31proc BelongSemig(bigint n, intvec v, list #)
32"
33USAGE:   BelongSemig (n,v[,sup]);  n bigint, v and sup intvec
34RETURN:  In the default form, it returns 1 if n is in the semigroup generated by
35         the elements of v or 0 otherwise. If the argument sup is added and in case
36         n belongs to the semigroup generated by the elements of v, it returns
37         a monomial in the variables {x(i) | i in sup} of degree n if we set
38         deg(x(sup[j])) = v[j].
39ASSUME:  v and sup positive integer vectors of same size, sup has no
40         repeated entries, x(i) has to be an indeterminate in the current ring for
41         all i in sup.
42EXAMPLE: example BelongSemig; shows some examples
43"
44{
45//--------------------------- initialisation ---------------------------------
46  int i, j, num;
47  bigint PartialSum;
48  num = size(v);
49  int e = size(#);
50
51  if  (e > 0)
52  {
53    intvec sup = #[1];
54    poly mon;
55  }
56
57  for (i = 1; i <= nrows(v); i++)
58  {
59    if ((n % v[i]) == 0)
60    {
61      // ---- n is multiple of v[i]
62      if (e)
63      {
64        mon = x(sup[i])^(int(n/v[i]));
65        return(mon);
66      }
67      else
68      {
69        return (1);
70      }
71    }
72  }
73
74  if (num == 1)
75  {
76    // ---- num = 1 and n is not multiple of v[1] --> FALSE
77    return(0);
78  }
79
80  intvec counter;
81  counter[num] = 0;
82  PartialSum = 0;
83
84  intvec w = sort(v)[1];
85  intvec cambio = sort(v)[2];
86
87  // ---- Iterative procedure to determine if n is in the semigroup generated by v
88  while (1)
89  {
90    if (n  >= PartialSum)
91    {
92      if (((n - PartialSum) % w[1]) == 0)
93      {
94        // ---- n belongs to the semigroup generated by v,
95        if (e)
96        {
97          // ---- obtain the monomial.
98          mon = x(sup[cambio[1]])^(int((n - PartialSum) / w[1]));
99          for (j = 2; j <= num; j++)
100          {
101            mon = mon * x(sup[cambio[j]])^(counter[j]);
102          }
103          return(mon);
104        }
105        else
106        {
107          // ---- returns true.
108          return (1);
109        }
110      }
111    }
112    i = num;
113    while (!defined(end))
114    {
115      if (i == 1)
116      {
117        // ---- Stop, n is not in the semigroup
118        return(0);
119      }
120      if (i > 1)
121      {
122        // counters control
123        if (counter[i] >= ((n - PartialSum) / w[i]))
124        {
125          PartialSum = PartialSum - (counter[i]*w[i]);
126          counter[i] = 0;
127          i--;
128        }
129        else
130        {
131          counter[i] = counter[i] + 1;
132          PartialSum = PartialSum + w[i];
133          int end;
134        }
135      }
136    }
137    kill end;
138  }
139}
140example
141{ "EXAMPLE:";
142   ring r=0,x(1..5),dp;
143   int a = 125;
144   intvec v = 13,17,51;
145   intvec sup = 2,4,1;
146   BelongSemig(a,v,sup);
147   BelongSemig(a,v);
148}
149
150///////////////////////////////////////////////////////////////////////////////////////////////////////////
151///////////////////////////////////////////////////////////////////////////////////////////////////////////
152///////////////////////////////////////////////////////////////////////////////////////////////////////////
153
154proc MinMult(int a, intvec b)
155"
156USAGE:   MinMult (a, b); a integer, b integer vector.
157RETURN:  an integer k, the minimum positive integer such that ka belongs to the
158         semigroup generated by the integers in b.
159ASSUME:  a is a positive integer, b is a positive integers vector.
160EXAMPLE: example MinMult; shows some examples.
161"
162{
163//--------------------------- initialisation ---------------------------------
164  int i, j, min, max;
165  int n = nrows(b);
166
167  if (n == 1)
168  {
169    // ---- trivial case
170    return(b[1]/gcd(a,b[1]));
171  }
172
173  max = b[1];
174  for (i = 2; i <= n; i++)
175  {
176    if (b[i] > max)
177    {
178      max = b[i];
179    }
180  }
181  int NumNodes = a + max; //----Number of nodes in the graph
182
183  int dist = 1;
184  // ---- Auxiliary structures to obtain the shortest path between the nodes 1 and a+1 of this graph
185  intvec queue = 1;
186  intvec queue2;
187
188  // ---- Control vector:
189  //      control[i] = 0 -> node not reached yet
190  //      control[i] = 1 -> node in queue1
191  //      control[i] = 2 -> node in queue2
192  //      control[i] = 3 -> node already processed
193  intvec control;
194  control[1] = 3;         // Starting node
195  control[a + max] = 0;   // Ending node
196  int current = 1;        // Current node
197  int next;               // Node connected to corrent by arc (current, next)
198
199  int ElemQueue, ElemQueue2;
200  int PosQueue = 1;
201
202  // Algoritmo de Dijkstra
203  while (1)
204  {
205    if (current <= a)
206    {
207      // ---- current <= a, arcs are (current, current + b[i])
208      for (i = 1; i <= n; i++)
209      {
210        next = current + b[i];
211        if (next == a+1)
212        {
213          kill control;
214          kill queue;
215          kill queue2;
216          return (dist);
217        }
218        if ((control[next] == 0)||(control[next] == 2))
219        {
220          control[next] = 1;
221          queue = queue, next;
222        }
223      }
224    }
225    if (current > a)
226    {
227      // ---- current > a, the only possible ars is (current, current - a)
228      next = current - a;
229      if (control[next] == 0)
230      {
231        control[next] = 2;
232        queue2[nrows(queue2) + 1] = next;
233      }
234    }
235    PosQueue++;
236    if (PosQueue <= nrows(queue))
237    {
238      current = queue[PosQueue];
239    }
240    else
241    {
242      dist++;
243      if (control[a+1] == 2)
244      {
245        return(dist);
246      }
247      queue = queue2[2..nrows(queue2)];
248      current = queue[1];
249      PosQueue = 1;
250      queue2 = 0;
251    }
252    control[current] = 3;
253  }
254}
255example
256{ "EXAMPLE:";
257   "int a = 46;";
258   "intvec b = 13,17,59;";
259   "MinMult(a,b);";
260   int a = 46;
261   intvec b = 13,17,59;
262   MinMult(a,b);
263   "// 3*a = 8*b[1] + 2*b[2]"
264}
265
266///////////////////////////////////////////////////////////////////////////////////////////////////////////
267///////////////////////////////////////////////////////////////////////////////////////////////////////////
268///////////////////////////////////////////////////////////////////////////////////////////////////////////
269
270proc CompInt(intvec d)
271"
272USAGE:   CompInt(d); d intvec.
273RETURN:  1 if the toric ideal I(d) is a complete intersection or 0 otherwise.
274ASSUME:  d is a vector of positive integers.
275NOTE:    If printlevel > 0, additional info is displayed in case
276         I(d) is a complete intersection:
277         if printlevel >= 1, it displays a minimal set of generators of the toric
278         ideal formed by quasihomogeneous binomials. Moreover, if printlevel >= 2
279         and gcd(d) = 1, it also shows the Frobenius number of the semigroup
280         generated by the elements in d.
281EXAMPLE: example CompInt; shows some examples
282"
283{
284//--------------------------- initialisation ---------------------------------
285
286  int i,j,k,l,divide,equal,possible;
287
288  int n = nrows(d);
289  int max = 2*n - 1;
290  ring r = 0, x(1..n), dp;
291
292  int level = printlevel - voice + 2;
293  // ---- To decide how much extra information calculate and display
294  if (level > 1)
295  {
296    int e = d[1];
297    for (i = 2; i <= n; i++)
298    {
299      e = gcd(e,d[i]);
300    }
301    if (e <> 1)
302    {
303      print ("// Semigroup generated by d is not numerical!");
304    }
305  }
306  if (level > 0)
307  {
308    ideal id;
309    vector mon;
310    mon[max] = 0;
311    if ((level > 1)&&(e == 1))
312    {
313      bigint frob = 0;
314    }
315  }
316
317  // ---- Trivial cases: n = 1,2 (it is a complete intersection).
318  if (n == 1)
319  {
320    print("// Ideal is (0)");
321    return (1);
322  }
323
324  if (n == 2)
325  {
326    if (level > 0)
327    {
328      intvec d1 = d[1];
329      intvec d2 = d[2];
330      int f1 = MinMult(d[1],d2);
331      int f2 = MinMult(d[2],d1);
332      id = x(1)^(f1) - x(2)^(f2);
333      print ("// Toric ideal:");
334      id;
335      if ((level > 1)&&(e == 1))
336      {
337        frob = d[1]*f1 - d[1] - d[2];
338        print ("// Frobenius number of the numerical semigroup:");
339        frob;
340      }
341    }
342    return (1);
343  }
344
345  // ---- For n >= 3 (non-trivial cases)
346  matrix mat[max][n];
347  intvec using, bound, multiple;
348  multiple[max] = 0;
349  bound[max] = 0;
350  using[max] = 0;
351
352  for (i = 1; i <= n; i++)
353  {
354    using[i] = 1;
355    multiple[i] = 0;
356    mat[i,i] = 1;
357  }
358  if (level > 1)
359  {
360    if (e == 1)
361    {
362      for (i = 1; i <= n; i++)
363      {
364        frob = frob - d[i];
365      }
366    }
367  }
368
369  int new, new1, new2;
370  for (i = 1; i <= n; i++)
371  {
372    for (j = 1; j < i; j++)
373    {
374      if (i <> j)
375      {
376        new = gcd(d[i],d[j]);
377        new1 = d[j]/new;
378        new2 = d[i]/new;
379        if (!bound[i] ||(new1 < bound[i]))
380        {
381          bound[i] = new1;
382        }
383        if (!bound[j] ||(new2 < bound[j]))
384        {
385          bound[j] = new2;
386        }
387      }
388    }
389  }
390
391  // ---- Begins the inductive part
392  for (i = 1; i < n; i++)
393  {
394    // ---- n-1 stages
395    for (j = 1; j < n + i; j++)
396    {
397      if ((using[j])&&(multiple[j] == 0))
398      {
399        possible = 0;
400        for (k = 1; (k < n + i)&&(!possible); k++)
401        {
402          if ((using[k])&&(k != j)&&(bigint(bound[k])*d[k] == bigint(bound[j])*d[j]))
403          {
404            possible = 1;
405          }
406        }
407        if (possible)
408        {
409          // ---- If possible == 1, then c_j has to be computed
410          intvec aux;
411          // ---- auxiliary vector containing all d[l] in use except d[j]
412          k = 1;
413          for (l = 1; l < n + i; l++)
414          {
415            if (using[l] && (l != j))
416            {
417              aux[k] = d[l];
418              k++;
419            }
420          }
421
422          multiple[j] = MinMult(d[j], aux);
423          kill aux;
424
425          if (j <= n)
426          {
427            if (level > 0)
428            {
429              mon = mon + (x(j)^multiple[j])*gen(j);
430            }
431          }
432          else
433          {
434            // ---- if j > n, it has to be checked if c_j belongs to a certain semigroup
435            intvec aux, sup;
436            k = 1;
437            for (l = 1; l <= n; l++)
438            {
439              if (mat[j, l] <> 0)
440              {
441                sup[k] = l;
442                aux[k] = d[l];
443                k++;
444              }
445            }
446            if (level > 0)
447            {
448              mon = mon + (BelongSemig(bigint(multiple[j])*d[j], aux, sup))*gen(j);
449              if (mon[j] == 0)
450              {
451                // ---- multiple[j]*d[j] does not belong to the semigroup generated by aux,
452                // ---- then it is NOT a complete intersection
453                return (0);
454              }
455            }
456            else
457            {
458              if (!BelongSemig(bigint(multiple[j])*d[j], aux))
459              {
460                // ---- multiple[j]*d[j] does not belong to the semigroup generated by aux,
461                // ---- then it is NOT a complete intersection
462                return (0);
463              }
464            }
465            kill sup;
466            kill aux;
467          }
468
469          // ---- Searching if there exist k such that multiple[k]*d[k]= multiple[j]*d[j]
470          equal = 0;
471          for (k = 1; k < n+i; k++)
472          {
473            if ((k <> j) && multiple[k] && using[k])
474            {
475              if (d[j]*bigint(multiple[j]) == d[k]*bigint(multiple[k]))
476              {
477                // found
478                equal = k;
479                break;
480              }
481            }
482          }
483          // ---- if equal = 0 no coincidence
484          if (!equal)
485          {
486            if (j == n + i - 1)
487            {
488              // ---- All multiple[k]*d[k] in use are different -> NOT complete intersection
489              return (0);
490            }
491          }
492          else
493          {
494            // ---- Next stage is prepared
495            if (level > 0)
496            {
497              //---- New generator of the toric ideal
498              id[i] = mon[j] - mon[equal];
499              if ((level > 1)&&(e == 1))
500              {
501                frob = frob + bigint(multiple[j])*d[j];
502              }
503            }
504            //---- Two exponents are removed and one is added
505            using[j] = 0;
506            using[equal] = 0;
507            using[n + i] = 1;
508            d[n + i] = gcd(d[j], d[equal]);  //---- new exponent
509            for (l = 1; l <= n; l++)
510            {
511              mat[n + i, l] = mat[j, l] + mat[equal, l];
512            }
513
514            // Bounds are reestablished
515            for (l = 1; l < n+i; l++)
516            {
517              if (using[l])
518              {
519                divide = gcd(d[l],d[n+i]);
520                new = d[n+i] / divide;
521                if ((multiple[l])&&(multiple[l] > new))
522                {
523                  return (0);
524                }
525                if (new < bound[l])
526                {
527                  bound[l] = new;
528                }
529                new = d[l] / divide;
530                if  ( !bound[n+i] || (new < bound[n+i]))
531                {
532                  bound[n+i] = new;
533                }
534              }
535            }
536            break;
537          }
538        }
539      }
540      if (j == n + i - 1)
541      {
542        // ---- All multiple[k]*d[k] in use are different -> NOT complete intersection
543        return (0);
544      }
545    }
546  }
547  if (level > 0)
548  {
549    "// Toric ideal: ";
550    id;
551    if ((level > 1)&&(e == 1))
552    {
553      "// Frobenius number of the numerical semigroup: ";
554      frob;
555    }
556  }
557  return(1);
558}
559example
560{ "EXAMPLE:";
561   printlevel = 0;
562   intvec d = 14,15,10,21;
563   CompInt(d);
564   printlevel = 3;
565   d = 36,54,125,150,225;
566   CompInt(d);
567   d = 45,70,75,98,147;
568   CompInt(d);
569};
570///////////////////////////////////////////////////////////////////////////////
571///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.