Return to toric.c CVS log

Up to [local] / OpenXM_contrib / TiGERS_0.9

Annotation of OpenXM_contrib/TiGERS_0.9/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:
                     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:   fprintf(stderr,"taking colon of J with %c:\n",'a'+lv);
                     38:
                     39:   /* set lex last variable to be lv */
                     40:   lold=ring_lv;
                     41:   ring_lv=lv;
                     42:
                     43:   /* compute grobner basis wrt to this deg rev lex term order */
                     44:   gset_rgb(g,monomial_rlexcomp);
                     45:
                     46:   /* now devide each binomial by highest power of lvar which divides it */
                     47:   for(ptr=gset_first(g); ptr!=0; ptr=binomial_next(ptr)){
                     48:     jk=min(binomial_lead(ptr)[lv],binomial_trail(ptr)[lv]);
                     49:     binomial_lead(ptr)[lv]-=jk;
                     50:     binomial_trail(ptr)[lv]-=jk;
                     51:   }
                     52:
                     53:   /* result is a grobner basis, reduce it to make it an rgb */
                     54:   gset_autoreduce(g);
                     55:
                     56:   /* restore original value of lvar*/
                     57:   ring_lv=lold;
                     58: }
                     59:
                     60: /*
                     61: ** gset gset_toric_ideal(int **M,int m, int n):
                     62: **  Given an mxn integer matrix M -- compute an rgb for the toric ideal I_M.
                     63: **  Uses repeated colon computations.
                     64: */
                     65: gset gset_toric_ideal(int **M,int m, int n){
                     66:    int i,j,jk,crk;
                     67:    binomial b;
                     68:    int **V,**S,**new_imatrix();
                     69:    gset g;
                     70:
                     71:    /* set ring dimension*/
                     72:    if (ring_N!=n){
                     73:     fprintf(stderr,"ERROR gset_toric_ideal():");
                     74:     fprintf(stderr," matrix dimensions incompatible with ring\n");
                     75:     return 0;
                     76:    }
                     77:    /* set wieghting of A */
                     78:    for(i=0;i<n;i++){
                     79:     ring_weight[i]=0;
                     80:     for(j=0;j<m;j++) ring_weight[i]+=M[j][i];
                     81:    }
                     82:
                     83:    /* reserve space for S and V */
                     84:    S=new_imatrix(n,m);
                     85:    V=new_imatrix(n,n);
                     86:
                     87:    /* copy transpose of M to S*/
                     88:    for(i=0;i<n;i++)for(j=0;j<m;j++) IMref(S,i,j)=IMref(M,j,i);
                     89:
                     90:    /* compute hermite normal form of M,S,V and corank of M*/
                     91:    crk=ihermite(S,V,m,n);
                     92:
                     93:    /* create new gset and read off equations corresponding to rows of V*/
                     94:    g=gset_new();
                     95:    for(i=crk;i<n;i++){
                     96:      b=binomial_new();
                     97:        for(j=0;j<n;j++){
                     98:          jk=IMref(V,i,j);
                     99:          if (jk>0){
                    100:              binomial_lead(b)[j]=jk; binomial_trail(b)[j]=0;}
                    101:          else {binomial_trail(b)[j]=-1*jk; binomial_lead(b)[j]=0;}
                    102:         }
                    103:      gset_insert(g,b);
                    104:    }
                    105:
                    106:    /* compute successive colon ideals */
                    107:    for(i=0;i<n;i++)gset_compute_colon(g,i);
                    108:
                    109:    /*unset weighting */
                    110:    for(i=0;i<n;i++) ring_weight[i]=1;
                    111:
                    112:    /* free space used by matrices */
                    113:    free_imatrix(M);
                    114:    free_imatrix(S);
                    115:    free_imatrix(V);
                    116:
                    117:    /* make sure result is the lex gbasis */
                    118:    gset_rgb(g,monomial_lexcomp);
                    119:
                    120: return g;
                    121: }
                    122:
                    123:
                    124:
                    125:
                    126:
                    127:
                    128:
                    129:
                    130:
                    131:
                    132:
                    133: