Changeset 987001 in git


Ignore:
Timestamp:
Jul 23, 2007, 2:50:07 PM (16 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
b529d9012a2f0139b1d63038125da10ba6fd7f90
Parents:
9866fc7d15f2e47654cbf50937ac78e9ab9c7a33
Message:
*hannes: doc., format


git-svn-id: file:///usr/local/Singular/svn/trunk@10227 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/cimonom.lib

    r9866fc r987001  
    1717  algorithm to check whether the toric ideal of an affine monomial curve is a complete intersection', Preprint.
    1818
     19SEE ALSO: Integer programming
     20
    1921PROCEDURES:
    2022 BelongSemig(n,v[,sup]);    checks whether n is in the semigroup generated by v;
    2123 MinMult(a,b);              computes k, the minimum positive integer such that k*a is in the semigroup of
    22                             possitive integers generated by the elements in b.
     24                            positive integers generated by the elements in b.
    2325 CompInt(d);                checks wether I(d) is a complete intersection or not.
    2426";
    … …  
    3032proc BelongSemig(bigint n, intvec v, list #)
    3133"
    32 USAGE:   BelongSemig (n,v[,sup]);  n big integer, v and sup integer vectors
     34USAGE:   BelongSemig (n,v[,sup]);  n bigint, v and sup intvec
    3335RETURN:  In the default form, it returns 1 if n is in the semigroup generated by
    3436         the elements of v or 0 otherwise. If the argument sup is added and in case
    … …  
    3638         a monomial in the variables {x(i) | i in sup} of degree n if we set
    3739         deg(x(sup[j])) = v[j].
    38 ASSUME:  n integer, v and sup possitive integer vectors of same size, sup has no
     40ASSUME:  v and sup positive integer vectors of same size, sup has no
    3941         repeated entries, x(i) has to be an indeterminate in the current ring for
    4042         all i in sup.
    … …  
    4345{
    4446//--------------------------- initialisation ---------------------------------
    45    int i, j, num;
    46    bigint PartialSum;
    47    num = size(v);
    48    int e = size(#);
    49 
    50    if  (e > 0)
    51    {
    52       intvec sup = #[1];
    53       poly mon;
    54    }
    55 
    56 
    57    for (i = 1; i <= nrows(v); i++) {
    58       if ((n % v[i]) == 0) {
    59          // ---- n is multiple of v[i]
    60          if (e) {
    61             mon = x(sup[i])^(int(n/v[i]));
    62             return(mon);
    63          }
    64          else {
    65             return (1);
    66          }
    67       }
    68    }
    69 
    70    if (num == 1) {
    71       // ---- num = 1 and n is not multiple of v[1] --> FALSE
    72       return(0);
    73    }
    74 
    75    intvec counter;
    76    counter[num] = 0;
    77    PartialSum = 0;
    78 
    79    intvec w = sort(v)[1];
    80    intvec cambio = sort(v)[2];
    81 
    82    // ---- Iterative procedure to determine if n is in the semigroup generated by v
    83    while (1) {
    84       if (n  >= PartialSum) {
    85          if (((n - PartialSum) % w[1]) == 0) {
    86             // ---- n belongs to the semigroup generated by v,
    87             if (e) {
    88                // ---- obtain the monomial.
    89                mon = x(sup[cambio[1]])^(int((n - PartialSum) / w[1]));
    90                for (j = 2; j <= num; j++) {
    91                   mon = mon * x(sup[cambio[j]])^(counter[j]);
    92                }
    93                return(mon);
    94             }
    95             else {
    96                // ---- returns true.
    97                return (1);
    98             }
    99          }
    100       }
    101       i = num;
    102       while (!defined(end)) {
    103          if (i == 1) {
    104             // ---- Stop, n is not in the semigroup
    105             return(0);
    106          }
    107          if (i > 1) {
    108             // counters control
    109             if (counter[i] >= ((n - PartialSum) / w[i])) {
    110                PartialSum = PartialSum - (counter[i]*w[i]);
    111                counter[i] = 0;
    112                i--;
    113             }
    114             else {
    115                counter[i] = counter[i] + 1;
    116                PartialSum = PartialSum + w[i];
    117                int end;
    118             }
    119          }
    120       }
    121       kill end;
    122    }
     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  }
    123140}
    124141example
    … …  
    139156"
    140157USAGE:   MinMult (a, b); a integer, b integer vector.
    141 RETURN:  an integer k, the minimum possitive integer such that ka belongs to the
     158RETURN:  an integer k, the minimum positive integer such that ka belongs to the
    142159         semigroup generated by the integers in b.
    143 ASSUME:  a is a possitive integer, b is a possitive integers vector.
     160ASSUME:  a is a positive integer, b is a positive integers vector.
    144161EXAMPLE: example MinMult; shows some examples.
    145162"
    146163{
    147164//--------------------------- initialisation ---------------------------------
    148    int i, j, min, max;
    149    int n = nrows(b);
    150 
    151    if (n == 1) {
    152       // ---- trivial case
    153       return(b[1]/gcd(a,b[1]));
    154    }
    155 
    156    max = b[1];
    157    for (i = 2; i <= n; i++) {
    158       if (b[i] > max) {
    159          max = b[i];
    160       }
    161    }
    162    int NumNodes = a + max; //----Number of nodes in the graph
    163 
    164    int dist = 1;
    165    // ---- Auxiliary structures to obtain the shortest path between the nodes 1 and a+1 of this graph
    166    intvec queue = 1;
    167    intvec queue2;
    168 
    169    // ---- Control vector:
    170    //      control[i] = 0 -> node not reached yet
    171    //      control[i] = 1 -> node in queue1
    172    //      control[i] = 2 -> node in queue2
    173    //      control[i] = 3 -> node already processed
    174    intvec control;
    175    control[1] = 3;         // Starting node
    176    control[a + max] = 0;   // Ending node
    177    int current = 1;        // Current node
    178    int next;               // Node connected to corrent by arc (current, next)
    179 
    180    int ElemQueue, ElemQueue2;
    181    int PosQueue = 1;
    182 
    183    // Algoritmo de Dijkstra
    184    while (1) {
    185       if (current <= a) {
    186          // ---- current <= a, arcs are (current, current + b[i])
    187          for (i = 1; i <= n; i++) {
    188             next = current + b[i];
    189             if (next == a+1) {
    190                kill control;
    191                kill queue;
    192                kill queue2;
    193                return (dist);
    194             }
    195             if ((control[next] == 0)||(control[next] == 2)) {
    196                control[next] = 1;
    197                queue = queue, next;
    198             }
    199          }
    200       }
    201       if (current > a) {
     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    {
    202228      // ---- current > a, the only possible ars is (current, current - a)
    203          next = current - a;
    204          if (control[next] == 0) {
    205             control[next] = 2;
    206             queue2[nrows(queue2) + 1] = next;
    207          }
    208       }
    209       PosQueue++;
    210       if (PosQueue <= nrows(queue)) {
    211          current = queue[PosQueue];
    212       }
    213       else {
    214          dist++;
    215          if (control[a+1] == 2) {
    216                 return(dist);
    217              }
    218              queue = queue2[2..nrows(queue2)];
    219              current = queue[1];
    220              PosQueue = 1;
    221              queue2 = 0;
    222       }
    223       control[current] = 3;
    224    }
     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  }
    225255}
    226256example
    … …  
    241271proc CompInt(intvec d)
    242272"
    243 USAGE:   CompInt(d); d integer vector.
     273USAGE:   CompInt(d); d intvec.
    244274RETURN:  1 if the toric ideal I(d) is a complete intersection or 0 otherwise.
    245 ASSUME:  d is a vector of possitive integers.
    246 NOTE:    If printlevel > 0 (default = 0), additional info is displayed in case
     275ASSUME:  d is a vector of positive integers.
     276NOTE:    If printlevel > 0, additional info is displayed in case
    247277         I(d) is a complete intersection:
    248278         if printlevel >= 1, it displays a minimal set of generators of the toric
    … …  
    255285//--------------------------- initialisation ---------------------------------
    256286
    257    int i,j,k,l,divide,equal,possible;
    258 
    259    int n = nrows(d);
    260    int max = 2*n - 1;
    261    ring r = 0, x(1..n), dp;
    262 
    263    int level = printlevel - voice + 2;
    264    // ---- To decide how much extra information calculate and display
    265    if (level > 1) {
    266       int e = d[1];
    267       for (i = 2; i <= n; i++) {
    268          e = gcd(e,d[i]);
    269       }
    270       if (e <> 1) {
    271          print ("// Semigroup generated by d is not numerical!");
    272       }
    273    }
    274    if (level > 0) {
    275       ideal id;
    276       vector mon;
    277       mon[max] = 0;
    278       if ((level > 1)&&(e == 1)) {
    279          bigint frob = 0;
    280       }
    281    }
    282 
    283    // ---- Trivial cases: n = 1,2 (it is a complete intersection).
    284    if (n == 1) {
    285       print("// Ideal is (0)");
    286       return (1);
    287    }
    288 
    289    if (n == 2) {
    290       if (level > 0) {
    291          intvec d1 = d[1];
    292          intvec d2 = d[2];
    293          int f1 = MinMult(d[1],d2);
    294          int f2 = MinMult(d[2],d1);
    295          id = x(1)^(f1) - x(2)^(f2);
    296          print ("// Toric ideal:");
    297          id;
    298          if ((level > 1)&&(e == 1)) {
    299             frob = d[1]*f1 - d[1] - d[2];
    300             print ("// Frobenius number of the numerical semigroup:");
    301             frob;
    302          }
    303       }
    304       return (1);
    305    }
    306 
    307    // ---- For n >= 3 (non-trivial cases)
    308    matrix mat[max][n];
    309    intvec using, bound, multiple;
    310    multiple[max] = 0;
    311    bound[max] = 0;
    312    using[max] = 0;
    313 
    314    for (i = 1; i <= n; i++) {
    315       using[i] = 1;
    316       multiple[i] = 0;
    317       mat[i,i] = 1;
    318    }
    319    if (level > 1) {
    320       if (e == 1) {
    321          for (i = 1; i <= n; i++) {
    322             frob = frob - d[i];
    323          }
    324       }
    325    }
    326 
    327 
    328    int new, new1, new2;
    329    for (i = 1; i <= n; i++) {
    330       for (j = 1; j < i; j++) {
    331          if (i <> j) {
    332             new = gcd(d[i],d[j]);
    333             new1 = d[j]/new;
    334             new2 = d[i]/new;
    335             if (!bound[i] ||(new1 < bound[i])) {
    336                bound[i] = new1;
    337             }
    338             if (!bound[j] ||(new2 < bound[j])) {
    339                bound[j] = new2;
    340             }
    341 
    342          }
    343       }
    344    }
    345 
    346 
    347    // ---- Begins the inductive part
    348    for (i = 1; i < n; i++) {
    349       // ---- n-1 stages
    350       for (j = 1; j < n + i; j++) {
    351          if ((using[j])&&(multiple[j] == 0)) {
    352                 possible = 0;
    353             for (k = 1; (k < n + i)&&(!possible); k++) {
    354                    if ((using[k])&&(k != j)&&(bigint(bound[k])*d[k] == bigint(bound[j])*d[j])) {
    355                           possible = 1;
    356                        }
     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);
    357525                }
    358             if (possible) {
    359                 // ---- If possible == 1, then c_j has to be computed
    360                intvec aux;
    361                // ---- auxiliary vector containing all d[l] in use except d[j]
    362                    k = 1;
    363                for (l = 1; l < n + i; l++) {
    364                   if (using[l] && (l != j)) {
    365                      aux[k] = d[l];
    366                      k++;
    367                   }
    368                }
    369 
    370                multiple[j] = MinMult(d[j], aux);
    371                kill aux;
    372 
    373                    if (j <= n) {
    374                       if (level > 0) {
    375                      mon = mon + (x(j)^multiple[j])*gen(j);
    376                   }
    377                    }
    378                else {
    379                // ---- if j > n, it has to be checked if c_j belongs to a certain semigroup
    380                   intvec aux, sup;
    381                           k = 1;
    382                   for (l = 1; l <= n; l++) {
    383                      if (mat[j, l] <> 0) {
    384                                 sup[k] = l;
    385                         aux[k] = d[l];
    386                         k++;
    387                              }
    388                   }
    389                   if (level > 0) {
    390                      mon = mon + (BelongSemig(bigint(multiple[j])*d[j], aux, sup))*gen(j);
    391                      if (mon[j] == 0) {
    392                         // ---- multiple[j]*d[j] does not belong to the semigroup generated by aux,
    393                         // ---- then it is NOT a complete intersection
    394                                 return (0);
    395                              }
    396                           }
    397                           else {
    398                              if (!BelongSemig(bigint(multiple[j])*d[j], aux)) {
    399                         // ---- multiple[j]*d[j] does not belong to the semigroup generated by aux,
    400                         // ---- then it is NOT a complete intersection
    401                                 return (0);
    402                              }
    403                           }
    404                   kill sup;
    405                   kill aux;
    406                    }
    407 
    408                    // ---- Searching if there exist k such that multiple[k]*d[k]= multiple[j]*d[j]
    409                equal = 0;
    410                for (k = 1; k < n+i; k++) {
    411                   if ((k <> j) && multiple[k] && using[k]) {
    412                      if (d[j]*bigint(multiple[j]) == d[k]*bigint(multiple[k])) {
    413                         // found
    414                         equal = k;
    415                         break;
    416                      }
    417                   }
    418                }
    419                // ---- if equal = 0 no coincidence
    420                if (!equal) {
    421                    if (j == n + i - 1) {
    422                       // ---- All multiple[k]*d[k] in use are different -> NOT complete intersection
    423                               return (0);
    424                    }
    425                    }
    426                else {
    427                   // ---- Next stage is prepared
    428                       if (level > 0) {
    429                          //---- New generator of the toric ideal
    430                      id[i] = mon[j] - mon[equal];
    431                      if ((level > 1)&&(e == 1)) {
    432                         frob = frob + bigint(multiple[j])*d[j];
    433                      }
    434                   }
    435                   //---- Two exponents are removed and one is added
    436                   using[j] = 0;
    437                   using[equal] = 0;
    438                   using[n + i] = 1;
    439                   d[n + i] = gcd(d[j], d[equal]);  //---- new exponent
    440                   for (l = 1; l <= n; l++) {
    441                      mat[n + i, l] = mat[j, l] + mat[equal, l];
    442                    }
    443 
    444                           // Bounds are reestablished
    445                           for (l = 1; l < n+i; l++) {
    446                              if (using[l]) {
    447                                 divide = gcd(d[l],d[n+i]);
    448                                 new = d[n+i] / divide;
    449                                    if ((multiple[l])&&(multiple[l] > new)) {
    450                                    return (0);
    451                                 }
    452                                 if (new < bound[l]) {
    453                                    bound[l] = new;
    454                                 }
    455                                 new = d[l] / divide;
    456                                 if  ( !bound[n+i] || (new < bound[n+i])) {
    457                                    bound[n+i] = new;
    458                                 }
    459                              }
    460                           }
    461                           break;
    462                }
    463             }
    464 
    465          }
    466          if (j == n + i - 1) {
    467             // ---- All multiple[k]*d[k] in use are different -> NOT complete intersection
    468                  return (0);
    469              }
    470       }
    471    }
    472    if (level > 0) {
    473       "// Toric ideal: ";
    474       id;
    475       if ((level > 1)&&(e == 1)) {
    476          "// Frobenius number of the numerical semigroup: ";
    477          frob;
    478       }
    479    }
    480    return(1);
    481 
     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);
    482559}
    483560example
Note: See TracChangeset for help on using the changeset viewer.