Return to toric.c CVS log

Up to [local] / OpenXM / src / Ti

Annotation of OpenXM/src/Ti/toric.c, Revision 1.1.1.1

1.1       maekawa     1: /*
                      2: ** toric.c                                 Birk Huber, 4/99
                      3: ** -- support for computing first grobner bases of toric ideals.
                      4: ** -- includes functions gset_rlex_rgb() and gset_compute_colon()
                      5: **    which once tested and working will move back to the gset.c file.
                      6: **
                      7: **
                      8: ** TiGERS,  Toric Groebner Basis Enumeration by Reverse Search
                      9: ** copyright (c) 1999  Birk Huber
                     10: **
                     11: */
                     12: #include<stdio.h>
                     13: #include<stdlib.h>
                     14: #include "utils.h"
                     15: #include "gset.h"
                     16: #include "matrices.h"
                     17: #include "Rsimp.h"
                     18: extern int NT;
                     19:
                     20:
                     21:
                     22: /*
                     23: ** void gset_compute_colon(gset g,int lv):
                     24: **    Compute deg rev lex rgb (lv least) for colon of g with variable lv.
                     25: **    - set ring_lv=lv.
                     26: **    - convert g to rgb wrt deg_rev_lex as implemented in binomial.h.
                     27: **    - divide each binomial by highest power of var lv possible.
                     28: **    - autoreduce.
                     29: **
                     30: */
                     31: #define min(m,n) (((m)<(n))?(m):(n))
                     32: void gset_compute_colon(gset g,int lv){
                     33:   binomial ptr;
                     34:   int lold,jk;
                     35:
                     36:   /* progress report */
                     37:   if (NT == 1) {
                     38:   }else{
                     39:     fprintf(stderr,"taking colon of J with %c:\n",'a'+lv);
                     40:   }
                     41:
                     42:   /* set lex last variable to be lv */
                     43:   lold=ring_lv;
                     44:   ring_lv=lv;
                     45:
                     46:   /* compute grobner basis wrt to this deg rev lex term order */
                     47:   gset_rgb(g,monomial_rlexcomp);
                     48:
                     49:   /* now devide each binomial by highest power of lvar which divides it */
                     50:   for(ptr=gset_first(g); ptr!=0; ptr=binomial_next(ptr)){
                     51:     jk=min(binomial_lead(ptr)[lv],binomial_trail(ptr)[lv]);
                     52:     binomial_lead(ptr)[lv]-=jk;
                     53:     binomial_trail(ptr)[lv]-=jk;
                     54:   }
                     55:
                     56:   /* result is a grobner basis, reduce it to make it an rgb */
                     57:   gset_autoreduce(g);
                     58:
                     59:   /* restore original value of lvar*/
                     60:   ring_lv=lold;
                     61: }
                     62:
                     63: /*
                     64: ** gset gset_toric_ideal(int **M,int m, int n):
                     65: **  Given an mxn integer matrix M -- compute an rgb for the toric ideal I_M.
                     66: **  Uses repeated colon computations.
                     67: */
                     68: gset gset_toric_ideal(int **M,int m, int n){
                     69:    int i,j,jk,crk;
                     70:    binomial b;
                     71:    int **V,**S,**new_imatrix();
                     72:    gset g;
                     73:
                     74:    /* set ring dimension*/
                     75:    if (ring_N!=n){
                     76:     fprintf(stderr,"ERROR gset_toric_ideal():");
                     77:     fprintf(stderr," matrix dimensions incompatible with ring\n");
                     78:     return 0;
                     79:    }
                     80:    /* set wieghting of A */
                     81:    for(i=0;i<n;i++){
                     82:     ring_weight[i]=0;
                     83:     for(j=0;j<m;j++) ring_weight[i]+=M[j][i];
                     84:    }
                     85:
                     86:    /* reserve space for S and V */
                     87:    S=new_imatrix(n,m);
                     88:    V=new_imatrix(n,n);
                     89:
                     90:    /* copy transpose of M to S*/
                     91:    for(i=0;i<n;i++)for(j=0;j<m;j++) IMref(S,i,j)=IMref(M,j,i);
                     92:
                     93:    /* compute hermite normal form of M,S,V and corank of M*/
                     94:    crk=ihermite(S,V,m,n);
                     95:
                     96:    /* create new gset and read off equations corresponding to rows of V*/
                     97:    g=gset_new();
                     98:    for(i=crk;i<n;i++){
                     99:      b=binomial_new();
                    100:        for(j=0;j<n;j++){
                    101:          jk=IMref(V,i,j);
                    102:          if (jk>0){
                    103:              binomial_lead(b)[j]=jk; binomial_trail(b)[j]=0;}
                    104:          else {binomial_trail(b)[j]=-1*jk; binomial_lead(b)[j]=0;}
                    105:         }
                    106:      gset_insert(g,b);
                    107:    }
                    108:
                    109:    /* compute successive colon ideals */
                    110:    for(i=0;i<n;i++)gset_compute_colon(g,i);
                    111:
                    112:    /*unset weighting */
                    113:    for(i=0;i<n;i++) ring_weight[i]=1;
                    114:
                    115:    /* free space used by matrices */
                    116:    free_imatrix(M);
                    117:    free_imatrix(S);
                    118:    free_imatrix(V);
                    119:
                    120:    /* make sure result is the lex gbasis */
                    121:    gset_rgb(g,monomial_lexcomp);
                    122:
                    123: return g;
                    124: }
                    125:
                    126:
                    127:
                    128:
                    129:
                    130:
                    131:
                    132:
                    133:
                    134:
                    135:
                    136: