Return to mvol.cpp CVS log

Up to [local] / OpenXM / src / ox_cdd

Annotation of OpenXM/src/ox_cdd/mvol.cpp, Revision 1.2

1.1       noro        1: /*     $Id: mvol.cpp,v 1.4 2004/11/12 04:37:15 fujiwara Exp $  */
                      2:
                      3: /*
                      4: * ##  origianl  #################################################
                      5: MixedVol: Polynomial system mixed volume calculation program.
                      6: * Version 1.0 Copyright (C) 2003  Tangan Gao, T. Y. Li, Xing Li, and Mengnien Wu.
                      7: * ###############################################################
                      8: */
                      9:
                     10: #include <cstdlib>
                     11: #include <string>
                     12: #include <ctime>
                     13: #include <fstream>
                     14: #include <iostream>
1.2     ! ohara      15: #include <limits.h>
1.1       noro       16:
                     17: #include "GetOption/getopt.h"
                     18: #include "PolyReader/PolynomialSystemReader.h"
                     19: #include "PolyReader/PolynomialException.h"
                     20: #include "PreProcess/Pre4MV.h"
                     21: #include "MixedVol/MixedVol.h"
                     22:
                     23: extern int debug_print;
                     24:
                     25: extern "C" int mvol(int nVar,int *Spt1stIdx,int **Spt)
                     26: {
                     27:        int i, j=0, k;
                     28:        int nPts = 0;
                     29:        int nSpt;
                     30:
                     31:        nPts = Spt1stIdx[nVar];
                     32:
                     33:        for (i=nVar-1; i>=0; i--){
                     34:                Spt1stIdx[i] = Spt1stIdx[i+1] - Spt1stIdx[i];
                     35:        }
                     36:
                     37:        // Quick return if 1-variable system or less than two terms
                     38:        k = -1;
                     39:        for (i=0; i<nVar; i++){
                     40:                if ( Spt1stIdx[i+1]-Spt1stIdx[i] < 2 )  {
                     41:                        k = i;
                     42:                        break;
                     43:                }
                     44:        }
                     45:        if ( k >= 0 ){
                     46:                if( debug_print >= 2 ){
                     47:                        cout << "The " << k+1 << "-th support has less than 2 points" << endl;
                     48:                }
                     49:                return 0;   // end of the case: too few terms
                     50:        }
                     51:
                     52:        if ( nVar == 1 ){
                     53:                int kmin, kmax;
                     54:                kmin = INT_MAX/2;
                     55:                kmax = -INT_MAX/2;
                     56:                for (i=0; i<Spt1stIdx[1]; i++){
                     57:                        kmax = max(kmax, Spt[i][0]);
                     58:                        kmin = min(kmin, Spt[i][0]);
                     59:                }
                     60:                if( debug_print >= 2 ){
                     61:                        cout << "A 1-dimensional support, its mixed volume is " << kmax-kmin << endl;
                     62:                }
                     63:                return (kmax-kmin);   // end of the case: 1-variable
                     64:        }
                     65:        for (i=0; i<Spt1stIdx[nVar]; i++){
                     66:                int kmin, kmax;
                     67:                kmin = INT_MAX/2;
                     68:                kmax = -INT_MAX/2;
                     69:                k=0;
                     70:                for (j=0; j<nVar; j++){
                     71:                        kmax = max(kmax, Spt[i][j]);
                     72:                        kmin = min(kmin, Spt[i][j]);
                     73:                        k += abs(Spt[i][j]);
                     74:                }
                     75:                if ( kmin != 0 || kmax > 1 || k > 1 ){
                     76:                        j = -1;
                     77:                        break;
                     78:                }
                     79:        }
                     80:        if ( j != -1 ){
                     81:                if( debug_print >= 2 ){
                     82:                        cout << "A linear support, its mixed volume <= 1" << endl;
                     83:                }
                     84:                return 1;
                     85:        }
                     86:        // end of quick return
                     87:
                     88:        // To preprocess the support
                     89:
                     90:        int* SptType = new int [nVar];
                     91:        int* SptVtx1stIdx = new int [nVar+1];
                     92:        int** SptVtx = new int* [nPts];
                     93:        SptVtx[0] = new int [nVar*nPts];
                     94:        for (i=1; i<nPts; i++) SptVtx[i] = SptVtx[i-1] + nVar;
                     95:        int* NuIdx2OldIdx = new int [nPts];
                     96:        nSpt = nVar;
                     97:
                     98:        Pre4MV(nVar,nSpt,SptType,Spt,Spt1stIdx,SptVtx,SptVtx1stIdx,NuIdx2OldIdx);
                     99:
                    100:        // To calculate the mixed volume of the support
                    101:
                    102:        srand( unsigned(time(0)) );
                    103:        double* lft = new double[SptVtx1stIdx[nSpt]];
                    104:        for (j=0; j<SptVtx1stIdx[nSpt]; j++){
                    105:                lft[j] = 1.0+(3*double(rand())-double(rand()))/RAND_MAX;
                    106:        }
                    107:
                    108:        int CellSize = nSpt;
                    109:        for (j=0; j<nSpt; j++ ){
                    110:                CellSize += SptType[j];
                    111:        }
                    112:        CellStack  MCells( CellSize );
                    113:        int MVol;
                    114:
                    115:        MixedVol(nVar,nSpt,SptType,SptVtx1stIdx,SptVtx,lft,MCells,MVol);
                    116:
                    117:        if( debug_print >= 2 ){
                    118:                cout << "The mixed volume is " << MVol << endl;
                    119:        }
                    120:        // Clean the memory
                    121:
                    122:        while( ! MCells.IsEmpty() ) MCells.Pop();
                    123:
                    124:        delete [] SptType;
                    125:        delete [] SptVtx1stIdx;
                    126:        delete [] SptVtx[0];
                    127:        delete [] SptVtx;
                    128:        delete [] NuIdx2OldIdx;
                    129:        delete [] lft;
                    130:
                    131:        return MVol;
                    132: }