source: git/Singular/LIB/cimonom.lib @ 1288ef

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