GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it

1
/*
2
 * Normaliz
3
 * Copyright (C) 2007-2014  Winfried Bruns, Bogdan Ichim, Christof Soeger
4
 * This program is free software: you can redistribute it and/or modify
5
 * it under the terms of the GNU General Public License as published by
6
 * the Free Software Foundation, either version 3 of the License, or
7
 * (at your option) any later version.
8
 *
9
 * This program is distributed in the hope that it will be useful,
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12
 * GNU General Public License for more details.
13
 *
14
 * You should have received a copy of the GNU General Public License
15
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
16
 *
17
 * As an exception, when this program is distributed through (i) the App Store
18
 * by Apple Inc.; (ii) the Mac App Store by Apple Inc.; or (iii) Google Play
19
 * by Google Inc., then that store may impose any digital rights management,
20
 * device limits and/or redistribution restrictions that are required by its
21
 * terms of service.
22
 */
23

24
//---------------------------------------------------------------------------
25

26
#include <boost/dynamic_bitset.hpp>
27

28
#include "libnormaliz/integer.h"
29
#include "libnormaliz/matrix.h"
30
#include "libnormaliz/nmz_nauty.h"
31
#include "libnormaliz/vector_operations.h"
32

33
#ifdef NMZ_NAUTY
34

35
// #define MAXN 5000    /* Define this before including nauty.h */
36
// we use dynamic allocation
37

38
#include <nauty/nauty.h>
39

40
namespace libnormaliz {
41
using namespace std;
42

43
vector<vector<long> > CollectedAutoms;
44

45
void getmyautoms(int count, int *perm, int *orbits,
46
               int numorbits, int stabvertex, int n){
47
    int i;
48
    vector<long> this_perm(n);
49
    for (i = 0; i < n; ++i) this_perm[i] = perm[i];
50
    CollectedAutoms.push_back(this_perm);
51
}
52

53
template<typename Integer>
54
vector<vector<long> > compute_automs_by_nauty(const vector<vector<Integer> >& Generators, size_t nr_special_gens, 
55
                                              const vector<vector<Integer> >& LinForms, 
56
                                              const size_t nr_special_linforms, mpz_class& group_order,
57
                                              BinaryMatrix<Integer>& CanType){
58
    CollectedAutoms.clear();
59
    
60
    DYNALLSTAT(graph,g,g_sz);
61
    DYNALLSTAT(graph,cg,cg_sz);
62
    DYNALLSTAT(int,lab,lab_sz);
63
    DYNALLSTAT(int,ptn,ptn_sz);
64
    DYNALLSTAT(int,orbits,orbits_sz);
65
    static DEFAULTOPTIONS_GRAPH(options);
66
    statsblk stats;
67
    
68
    options.userautomproc = getmyautoms;
69
    options.getcanon = TRUE;
70
    
71
    int n,m;
72
    
73
    options.writeautoms = FALSE;
74
    options.defaultptn = FALSE;
75
    
76
    size_t mm=Generators.size();
77
    size_t nn=LinForms.size();
78
    
79
    BinaryMatrix<Integer> MM(mm,nn);
80
    Matrix<Integer> MVal(mm,nn);
81
    
82
    key_t i,j,k;
83

84
    bool first=true;
85
    Integer mini=0;
86
    for(i=0;i<mm; ++i){
87
        for(j=0;j<nn;++j){
88
            MVal[i][j]=v_scalar_product(Generators[i],LinForms[j]);
89
            if(MVal[i][j]<mini || first){
90
                mini=MVal[i][j];
91
                first=false;
92
            }
93
        }
94
    }
95
        
96
    MM.set_offset(mini);
97
    for(i=0;i<mm; ++i){
98
        for(j=0;j<nn;++j)
99
            MM.insert(MVal[i][j]-mini,i,j);
100
    }
101
        
102
    size_t ll=MM.nr_layers();
103
        
104
    size_t layer_size=mm+nn;
105
    n=ll*layer_size;
106
    m = SETWORDSNEEDED(n);
107
    
108
    nauty_check(WORDSIZE,m,n,NAUTYVERSIONID);
109
    
110
    DYNALLOC2(graph,g,g_sz,m,n,"malloc");
111
    DYNALLOC2(graph,cg,cg_sz,n,m,"malloc");
112
    DYNALLOC1(int,lab,lab_sz,n,"malloc");
113
    DYNALLOC1(int,ptn,ptn_sz,n,"malloc");
114
    DYNALLOC1(int,orbits,orbits_sz,n,"malloc");
115
    
116
    EMPTYGRAPH(g,m,n);
117
    
118
    for(i=0;i<layer_size;++i){   // make vertical edges over all layers
119
        for(k=1;k<ll;++k)
120
            ADDONEEDGE(g,(k-1)*layer_size+i,k*layer_size+i,m);
121
    }
122
    
123
    for(i=0;i<mm;++i){   // make horizontal edges layer by layer
124
        for(j=0;j<nn;++j){
125
            for(k=0;k<ll;++k){
126
                if(MM.test(i,j,k))  // k is the number of layers below the current one
127
                    ADDONEEDGE(g,k*layer_size+i,k*layer_size+mm+j,m);
128
            }
129
        }
130
    }           
131
    
132
    for(int ii=0;ii<n;++ii){ // prepare partitions
133
        lab[ii]=ii;
134
        ptn[ii]=1;
135
    }
136
    
137
    for(k=0;k<ll;++k){ // make partitions layer by layer
138
        ptn[k*layer_size+ mm-1]=0; // row vertices in one partition
139
        for(size_t s=0; s< nr_special_gens;++s) // speciall generators in extra partitions (makes them fixed points)
140
            ptn[k*layer_size+s]=0;
141
        ptn[(k+1)*layer_size-1]=0; // column indices in the next
142
        for(size_t s=0; s< nr_special_linforms;++s) // special linear forms in extra partitions
143
            ptn[(k+1)*layer_size-2-s]=0;            
144
    } 
145

146
    densenauty(g,lab,ptn,orbits,&options,&stats,m,n,cg);
147
    
148
    vector<vector<long> > AutomsAndOrbits(2*CollectedAutoms.size());
149
    AutomsAndOrbits.reserve(2*CollectedAutoms.size()+3);
150

151
    for(k=0;k<CollectedAutoms.size();++k){
152
        vector<long> GenPerm(mm);
153
        for(i=0;i<mm;++i)
154
            GenPerm[i]=CollectedAutoms[k][i];
155
        AutomsAndOrbits[k]=GenPerm;
156
        vector<long> LFPerm(nn-nr_special_linforms);  // we remove the special linear forms here
157
        for(i=mm;i<mm+nn-nr_special_linforms;++i)
158
            LFPerm[i-mm]=CollectedAutoms[k][i]-mm;
159
        AutomsAndOrbits[k+CollectedAutoms.size()]=LFPerm;
160
        AutomsAndOrbits[k+CollectedAutoms.size()];
161
        
162
    }    
163
    
164
    vector<long> GenOrbits(mm);
165
    for(i=0;i<mm;++i)
166
        GenOrbits[i]=orbits[i];
167
    AutomsAndOrbits.push_back(GenOrbits);
168
    
169
    vector<long> LFOrbits(nn-nr_special_linforms); // we remove the special linear forms here
170
    for(i=0;i<nn-nr_special_linforms;++i)
171
        LFOrbits[i]=orbits[i+mm]-mm;
172
    AutomsAndOrbits.push_back(LFOrbits);
173
 
174
    group_order=mpz_class(stats.grpsize1);
175
    
176
    vector<long> row_order(mm), col_order(nn);
177
    for(key_t i=0;i<mm;++i)
178
        row_order[i]=lab[i];
179
    for(key_t i=0;i<nn;++i)
180
        col_order[i]=lab[mm+i]-mm;
181
    
182
    AutomsAndOrbits.push_back(row_order);
183
    
184
    nauty_freedyn();
185
    
186
    CanType=MM.reordered(row_order,col_order);
187
    
188
    return AutomsAndOrbits;
189
        
190
}
191

192
#ifndef NMZ_MIC_OFFLOAD  //offload with long is not supported
193
template vector<vector<long> > compute_automs_by_nauty(const vector<vector<long> >& Generators, size_t nr_special_gens, 
194
                        const vector<vector<long> >& LinForms,const size_t nr_special_linforms,   
195
                        mpz_class& group_order, BinaryMatrix<long>& CanType);
196
#endif // NMZ_MIC_OFFLOAD
197
template vector<vector<long> > compute_automs_by_nauty(const vector<vector<long long> >& Generators, size_t nr_special_gens, 
198
                        const vector<vector<long long> >& LinForms,const size_t nr_special_linforms,   
199
                        mpz_class& group_order, BinaryMatrix<long long>& CanType);
200
template vector<vector<long> > compute_automs_by_nauty(const vector<vector<mpz_class> >& Generators, size_t nr_special_gens, 
201
                        const vector<vector<mpz_class> >& LinForms,const size_t nr_special_linforms,   
202
                        mpz_class& group_order,BinaryMatrix<mpz_class>& CanType);
203

204
} // namespace
205

206
#endif // NMZ_NAUTY
207

208

209