diff --git a/CMakeLists.txt b/CMakeLists.txt
index 167fbf3993ae0a9b20b106c1c5a205df7f586588..0831fd83477efe0f9e9976a00149d555988d68e5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -83,7 +83,7 @@ set(GMSH_API
Post/PView.h Post/PViewData.h Plugin/PluginManager.h
Graphics/drawContext.h
contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h
- contrib/kbipack/gmp_blas.h
+ contrib/kbipack/gmp_blas.h contrib/kbipack/mpz.h
contrib/DiscreteIntegration/DILevelset.h)
execute_process(COMMAND date "+%Y%m%d" OUTPUT_VARIABLE DATE
@@ -427,9 +427,11 @@ endif(ENABLE_GMM)
if(ENABLE_KBIPACK)
find_library(GMP_LIB NAMES libgmp.a libgmp.dll.a)
if(GMP_LIB)
- add_subdirectory(contrib/kbipack)
+ message("-- Found GMP")
set(HAVE_KBIPACK TRUE)
list(APPEND CONFIG_OPTIONS "Kbipack")
+ add_subdirectory(contrib/kbipack)
+ set(HAVE_GMP TRUE)
list(APPEND EXTERNAL_LIBRARIES ${GMP_LIB})
endif(GMP_LIB)
endif(ENABLE_KBIPACK)
diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index 97c42ba0b6ac0d49c95c38a9528055e6e146bc0b..55b80babe31f3410416433fa32238d4a5ff295f5 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -15,6 +15,7 @@
#cmakedefine HAVE_FL_TREE
#cmakedefine HAVE_FOURIER_MODEL
#cmakedefine HAVE_GMM
+#cmakedefine HAVE_GMP
#cmakedefine HAVE_KBIPACK
#cmakedefine HAVE_LAPACK
#cmakedefine HAVE_LIBCGNS
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 87a8515e670f9008227eceee1283d44dae7c867e..51a67ebd3c01f6acc130fc0342aa96b3bd102ea0 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -348,7 +348,7 @@ int CellComplex::coreduction(Cell* generator, int omitted){
int CellComplex::reduction(int dim, int omitted){
if(dim < 1 || dim > 3) return 0;
- std::map<Cell*, int, Less_Cell > cbd_c;
+
int count = 0;
bool reduced = true;
@@ -377,90 +377,38 @@ int CellComplex::reduction(int dim, int omitted){
}
return count;
}
-/*
-int CellComplex::reduction(Cell* generator){
-
- int coreductions = 0;
- std::queue<Cell*> Q;
- std::set<Cell*, Less_Cell> Qset;
-
- Q.push(generator);
- Qset.insert(generator);
-
-
- std::list<Cell*> cbd_s;
- std::list<Cell*> bd_c;
- Cell* s;
- int round = 0;
- while( !Q.empty() ){
- round++;
- //printf("%d ", round);
-
- s = Q.front();
- Q.pop();
- removeCellQset(s, Qset);
-
- cbd_s = s->getCoboundary();
- if( cbd_s.size() == 1 && inSameDomain(s, cbd_s.front()) ){
- removeCell(s);
- bd_c = bd_s.front()->getBoundary();
- enqueueCells(cbd_c, Q, Qset);
- removeCell(cbd_s.front());
- reductions++;
- }
- else if(cbd_s.empty()){
- bd_c = s->getBoundary();
- enqueueCells(bd_c, Q, Qset);
- }
-
-
- }
- //printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
- return reductions;
-}
-*/
-/*
-int CellComplex::reduction(int dim){
+int CellComplex::coreduction(int dim, int omitted){
if(dim < 1 || dim > 3) return 0;
- std::list<Cell*> cbd_c;
- std::list<Cell*> bd_c;
+
int count = 0;
-
- std::queue<Cell*> Q;
- std::set<Cell*, Less_Cell> Qset;
-
-
- for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
-
- Cell* cell = *cit;
- Q.push(cell);
- Qset.insert(cell);
-
- while(Q.size() != 0){
-
- Cell* s = Q.front();
- Q.pop();
- cbd_c = s->getCoboundary();
- if( cbd_c.size() == 1 && inSameDomain(s, cbd_c.front()) ){
-
- removeCell(cbd_c.front());
- removeCell(s);
+
+ bool reduced = true;
+ bool ignoreCells = true;
+ while (reduced){
+
+ reduced = false;
+ citer cit = firstCell(dim);
+ while(cit != lastCell(dim)){
+ Cell* cell = *cit;
+ if( cell->getBoundarySize() == 1
+ && inSameDomain(cell, cell->firstBoundary()->first)){
+ ++cit;
+ if(dim-1 == 0 && omitted > 0){
+ _store.at(omitted-1).insert(cell->firstBoundary()->first);
+ }
+ removeCell(cell->firstBoundary()->first);
+ removeCell(cell);
count++;
- bd_c = cbd_c.front()->getBoundary();
- enqueueCells(bd_c, Q, Qset);
-
- cit = firstCell(dim-1);
+ reduced = true;
}
-
- removeCellQset(s, Qset);
-
+
+ if(getSize(dim) == 0 || getSize(dim-1) == 0) break;
+ cit++;
}
}
-
return count;
}
-*/
int CellComplex::reduceComplex(){
@@ -552,18 +500,21 @@ int CellComplex::coreduceComplex()
void CellComplex::computeBettiNumbers(){
+ removeSubdomain();
+
for(int i = 0; i < 4; i++){
- _betti[i] = 0;
- Msg::Debug("Betti number computation process: step %d of 4 \n", i+1);
+ if (getSize(i) != 0) _betti[i] = -1;
+ else _betti[i] = 0;
+ Msg::Debug("Betti number computation process: step %d of 4 \n", i+1);
while (getSize(i) != 0){
citer cit = firstCell(i);
Cell* cell = *cit;
- while(!cell->inSubdomain() && cit != lastCell(i)){
+ /*while(!cell->inSubdomain() && cit != lastCell(i)){
cell = *cit;
cit++;
}
- if(!cell->inSubdomain()) _betti[i] = _betti[i] + 1;
+ if(!cell->inSubdomain())*/ _betti[i] = _betti[i] + 1;
removeCell(cell, false);
coreduction(cell);
}
@@ -705,86 +656,6 @@ int CellComplex::combine(int dim){
return count;
}
-/*
-int CellComplex::combine(int dim){ // version of combine that doesn't have a queue
- double t1 = Cpu();
- printf("Cell complex before combining: %d volumes, %d faces, %d edges and %d vertices.\n",
- getSize(3), getSize(2), getSize(1), getSize(0));
-
- if(dim < 1 || dim > 3) return 0;
-
- std::list< std::pair<int, Cell*> > cbd_c;
- std::list<Cell*> bd_c;
- int count = 0;
- bool combined = true;
- while(combined){
- combined = false;
- for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
- Cell* s = *cit;
- cbd_c = s->getCoboundary();
-
- if(s->getCoboundarySize() == 2 && !(*(cbd_c.front().second) == *(cbd_c.back().second))
- && inSameDomain(s, cbd_c.front().second) && inSameDomain(s, cbd_c.back().second)
- && cbd_c.front().second->getNumVertices() < getSize(dim) // heuristics for mammoth cell birth control
- && cbd_c.back().second->getNumVertices() < getSize(dim)){
- int or1 = cbd_c.front().first;
- int or2 = cbd_c.back().first;
- Cell* c1 = cbd_c.front().second;
- Cell* c2 = cbd_c.back().second;
-
- removeCell(s);
-
- CombinedCell* newCell = new CombinedCell(c1, c2, (or1 != or2) );
- removeCell(c1);
- removeCell(c2);
- std::pair<citer, bool> insertInfo = _cells[dim].insert(newCell);
- if(!insertInfo.second) printf("Warning: Combined cell not inserted! \n");
- cit = firstCell(dim-1);
- count++;
- combined = true;
- }
- }
-
- }
-
- double t2 = Cpu();
- printf("Cell complex after combining: %d volumes, %d faces, %d edges and %d vertices (%g s).\n",
- getSize(3), getSize(2), getSize(1), getSize(0), t2 - t1);
-
-
- return count;
-}
-*/
-
-/*
-void CellComplex::swapSubdomain(){
-
- for(int i = 0; i < 4; i++){
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- if(cell->onDomainBoundary() && cell->inSubdomain()) cell->setInSubdomain(false);
- else if(cell->onDomainBoundary() && !cell->inSubdomain()) cell->setInSubdomain(true);
- }
- }
-
- // make subdomain closed
- for(int i = 0; i < 4; i++){
- for(citer cit = firstCell(i); cit != lastCell(i); cit++){
- Cell* cell = *cit;
- if(cell->inSubdomain()){
- std::list<Cell*> boundary = cell->getBoundary();
- for(std::list<Cell*>::iterator it = boundary.begin(); it != boundary.end(); it++){
- Cell* bdCell = *it;
- bdCell->setInSubdomain(true);
- }
- }
- }
- }
-
- return;
-}*/
-
-
void CellComplex::printComplex(int dim){
for (citer cit = firstCell(dim); cit != lastCell(dim); cit++){
Cell* cell = *cit;
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index d1a59bbb0eabe35c49dcc53abb8b4525a6970b34..058d85e57a4743fe8877557df8520c4faf8ced4a 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -121,11 +121,12 @@ class CellComplex
return ( (!c1->inSubdomain() && !c2->inSubdomain())
|| (c1->inSubdomain() && c2->inSubdomain()) ); }
- // reduction of this cell complex
- // removes reduction pairs of cell of dimension dim and dim-1
+ // (co)reduction of this cell complex
+ // removes (co)reduction pairs of cell of dimension dim and dim-1
int reduction(int dim, int omitted=0);
-
- // useful functions for (co)reduction of cell complex
+ int coreduction(int dim, int omitted=0);
+
+ // full (co)reduction of this cell complex (all dimensions)
int reduceComplex();
int coreduceComplex();
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index 9a9d75d6a6173575a36de3dd089592ebadf1196d..0d600737302583c0af45514b3a307ec02c121f20 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -68,10 +68,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
}
else{
-
mpz_t elem;
mpz_init(elem);
-
_HMatrix[dim] = create_gmp_matrix_zero(rows, cols);
for( std::set<Cell*, Less_Cell>::iterator cit =
cellComplex->firstCell(dim);
@@ -106,8 +104,8 @@ ChainComplex::ChainComplex(CellComplex* cellComplex){
}
}
}
- }
- mpz_clear(elem);
+ }
+ mpz_clear(elem);
}
_kerH[dim] = NULL;
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 5b85d8a1d867c64bbf055297ebaf5d7d88704a10..20546dea01cc3e8c6a2a08eb480fc1503c404f16 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -25,9 +25,9 @@
#include "GVertex.h"
#include "CellComplex.h"
-#include "gmp.h"
+#include "mpz.h"
extern "C" {
- #include "gmp_normal_form.h" // perhaps make c++ headers instead?
+ #include "gmp_normal_form.h"
}
// A class representing a chain complex of a cell complex.
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index 1ca2f803e7dbc39e20ad73ab3a7dfec99270ffe7..9d54dffef6eda4e0b87fc01297e6cf516161747f 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -199,8 +199,11 @@ void Homology::findDualGenerators(std::string fileName){
int omitted = _cellComplex->coreduceComplex();
_cellComplex->cocombine(0);
+ _cellComplex->coreduction(1);
_cellComplex->cocombine(1);
+ _cellComplex->coreduction(2);
_cellComplex->cocombine(2);
+ _cellComplex->coreduction(3);
_cellComplex->checkCoherence();
double t2 = Cpu();
diff --git a/contrib/kbipack/CMakeLists.txt b/contrib/kbipack/CMakeLists.txt
index 36e436e5c65a8034446bb6f928974e58b6831ed7..7ac9260016af0013be1e66df62416831ae3f5400 100644
--- a/contrib/kbipack/CMakeLists.txt
+++ b/contrib/kbipack/CMakeLists.txt
@@ -8,6 +8,7 @@ set(SRC
gmp_matrix_io.c
gmp_matrix.c
gmp_blas.c
+ mpz.c
)
file(GLOB_RECURSE HDR RELATIVE ${CMAKE_SOURCE_DIR}/contrib/kbipack *.h)
diff --git a/contrib/kbipack/gmp_blas.c b/contrib/kbipack/gmp_blas.c
index 501c4eba543f02980b4138c039bbbb8e920913db..6003888d5b67598c71dd23a44d50e5702f2c0feb 100644
--- a/contrib/kbipack/gmp_blas.c
+++ b/contrib/kbipack/gmp_blas.c
@@ -164,7 +164,7 @@ size_t gmp_blas_iamax(size_t n, const mpz_t* x, size_t incx)
mpz_t max_so_far;
mpz_init(max_so_far);
- mpz_set(max_so_far, 0);
+ mpz_set_si(max_so_far, 0);
ind_so_far = 0;
for(ind = 0; ind < n; ind++)
diff --git a/contrib/kbipack/gmp_blas.h b/contrib/kbipack/gmp_blas.h
index 6c9c88b2f8454539f9ba694239d52232b33e34df..b2800e5c2fa930627c3b74366d25c16d38fbd8a1 100644
--- a/contrib/kbipack/gmp_blas.h
+++ b/contrib/kbipack/gmp_blas.h
@@ -31,7 +31,7 @@
#define __GMP_BLAS_H__
#include <stdio.h>
-#include <gmp.h>
+#include <mpz.h>
/***********/
/* Level 1 */
@@ -39,7 +39,7 @@
/* x <-> y */
void
-gmp_blas_swap(size_t n, mpz_t * x, size_t incx, mpz_t * y, size_t incy);
+gmp_blas_swap(size_t n, mpz_t * x, size_t incx, mpz_t * y, size_t incy);
/* x <- ax */
void
diff --git a/contrib/kbipack/gmp_matrix.c b/contrib/kbipack/gmp_matrix.c
index 531466abd1372dc38931dd69152b3cb2c7486dca..c690292c6c70a4c07757ea8a974761aae76f4938 100644
--- a/contrib/kbipack/gmp_matrix.c
+++ b/contrib/kbipack/gmp_matrix.c
@@ -132,7 +132,7 @@ create_gmp_matrix_identity(size_t dim)
gmp_matrix *
create_gmp_matrix_zero(size_t rows, size_t cols)
-{
+{
gmp_matrix * new_matrix;
size_t ind;
@@ -156,7 +156,6 @@ create_gmp_matrix_zero(size_t rows, size_t cols)
{
mpz_init_set_si (new_matrix -> storage[ind], 0);
}
-
return new_matrix;
}
@@ -374,7 +373,6 @@ gmp_matrix_negate_col(size_t col, gmp_matrix * m)
{
return EXIT_FAILURE;
}
-
mpz_init(minus_one);
mpz_set_si(minus_one, -1);
gmp_blas_scal(m -> rows, minus_one,
diff --git a/contrib/kbipack/gmp_matrix.h b/contrib/kbipack/gmp_matrix.h
index 5d0e7c57dc853238729492548f1da24eb634b837..ea0c9af66b339b88684c3c438442bb554db5c0c1 100644
--- a/contrib/kbipack/gmp_matrix.h
+++ b/contrib/kbipack/gmp_matrix.h
@@ -28,7 +28,7 @@
#ifndef __GMP_MATRIX_H__
#define __GMP_MATRIX_H__
-#include"gmp_blas.h"
+#include "gmp_blas.h"
typedef struct
{
diff --git a/contrib/kbipack/mpz.c b/contrib/kbipack/mpz.c
new file mode 100755
index 0000000000000000000000000000000000000000..09c40df9610e6c24c6c802f68adb5833161268e9
--- /dev/null
+++ b/contrib/kbipack/mpz.c
@@ -0,0 +1,238 @@
+/*
+ This is an interface to masquerade ordinary long integers in C
+ language behind an interface mimicing arbitrary precision arithmetic
+ interger (mpz) interface in The GNU MP Bignum library.
+
+ Only the functions needed by Kbipack are available.
+
+ THIS INTERFACE IS NOT BERFORMING ANY ARBITRARY PRECISION ARITHMETIC.
+ You should always use the GMP library (http://gmplib.org/) when possible.
+
+ Some basic techniques to detect integer overflows are implemented.
+ However, they don't interfere the program flow IN ANY WAY, they just write
+ an error message to stdout.
+
+ */
+
+#include "mpz.h"
+
+// disabled for now
+#if 0
+
+#include "limits.h"
+
+void overflow()
+{
+ printf("ERROR: Integer overflow detected! Compile with GMP library to fix this. \n");
+}
+
+long int addcheck(long int a, long int b){
+
+ long int c = a + b;
+ if (b >= 0){
+ if (c < a) overflow();
+ }
+ else if (b < 0){
+ if (c > a) overflow();
+ }
+ return c;
+}
+
+long int multcheck(long int a, long int b){
+ long long int c = (long long int)a*b;
+ if(c > LONG_MAX || c < LONG_MIN) overflow();
+ return c;
+}
+
+void mpz_init(mpz_ptr x)
+{
+ x->z = 0;
+}
+
+void mpz_init_set(mpz_ptr rop, mpz_ptr op)
+{
+ rop->z = op->z;
+}
+
+void mpz_init_set_si(mpz_ptr rop, signed long int op)
+{
+ rop->z = op;
+}
+
+void mpz_set_ui(mpz_ptr rop, unsigned long int op)
+{
+ rop->z = op;
+}
+
+void mpz_set(mpz_ptr rop, mpz_ptr op)
+{
+ rop->z = op->z;
+}
+void mpz_set_si(mpz_ptr rop, signed long int op)
+{
+ rop->z = op;
+}
+
+void mpz_clear(mpz_ptr x)
+{
+
+}
+
+// arithmethic
+void mpz_add(mpz_ptr rop, mpz_ptr op1, mpz_ptr op2)
+{
+ rop->z = addcheck(op1->z, op2->z);
+}
+
+void mpz_mul(mpz_ptr rop, mpz_ptr op1, mpz_ptr op2)
+{
+ rop->z = multcheck(op1->z, op2->z);
+}
+
+void mpz_addmul(mpz_ptr rop, mpz_ptr op1, mpz_ptr op2)
+{
+ rop->z = addcheck(rop->z,multcheck(op1->z,op2->z));
+}
+
+void mpz_neg(mpz_ptr rop, mpz_ptr op)
+{
+ rop->z = -op->z;
+}
+
+
+// division
+void mpz_divexact(mpz_ptr q, mpz_ptr n, mpz_ptr d)
+{
+ div_t temp;
+ temp = div(n->z, d->z);
+ q->z = temp.quot;
+}
+
+void mpz_cdiv_q(mpz_ptr q, mpz_ptr n, mpz_ptr d)
+{
+ div_t temp;
+ temp = div(n->z, d->z);
+ q->z = temp.quot;
+}
+void mpz_cdiv_qr(mpz_ptr q, mpz_ptr r, mpz_ptr n, mpz_ptr d)
+{
+ div_t temp;
+ temp = div(n->z, d->z);
+ q->z = temp.quot;
+ r->z = temp.rem;
+}
+
+void mpz_tdiv_r(mpz_ptr r, mpz_ptr n, mpz_ptr d)
+{
+ div_t temp;
+ temp = div(n->z, d->z);
+ if(n->z < 0) r->z = -temp.rem;
+ else r->z = temp.rem;
+}
+
+// compare
+int mpz_cmp_si(mpz_ptr op1, signed long int op2)
+{
+ if(op1->z > op2) return 1;
+ else if(op1->z < op2) return -1;
+ else return 0;
+ return 0;
+}
+
+int mpz_cmpabs(mpz_ptr op1, mpz_ptr op2)
+{
+ if(abs(op1->z) > abs(op2->z)) return 1;
+ else if(abs(op1->z) < abs(op2->z)) return -1;
+ else return 0;
+ return 0;
+}
+
+int mpz_sgn(mpz_ptr op)
+{
+ if(op->z > 0) return 1;
+ else if(op->z < 0) return -1;
+ else return 0;
+}
+
+signed long int mpz_get_si(mpz_ptr op)
+{
+ signed long int temp = op->z;
+ return temp;
+}
+
+long int gcd(long int a, long int b)
+{
+ long int temp = 0;
+ long int remainder = 0;
+ if (a < b) {
+ temp = a;
+ a = b;
+ b = temp;
+ }
+ remainder = a % b;
+ if(remainder == 0) return b;
+ else return gcd(b, remainder);
+}
+
+void extended_gcd(long int* g, long int* s, long int* t,
+ long int a, long int b)
+{
+ long int x = 0;
+ long int lastx = 1;
+ long int y = 1;
+ long int lasty = 0;
+ while (b != 0){
+ div_t divt = div(a,b);
+ long int quotient = divt.quot;
+
+ long int temp = b;
+ b = a % b;
+ a = temp;
+
+ temp = x;
+ x = addcheck(lastx, multcheck(-quotient,x));
+ lastx = temp;
+
+ temp = y;
+ y = addcheck(lasty, multcheck(-quotient,y));
+ lasty = temp;
+ }
+ *s = lastx;
+ *t = lasty;
+ *g = a;
+}
+
+void mpz_gcdext(mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_ptr a, mpz_ptr b)
+{
+
+ extended_gcd(&g->z, &s->z, &t->z, a->z, b->z);
+ /*printf("g: %ld, s: %ld, t: %ld, a: %ld, b: %ld. \n",
+ g->z, s->z, t->z, a->z, b->z); */
+}
+
+size_t mpz_out_str(FILE *stream, int base, mpz_ptr op)
+{
+ return 0;
+}
+
+/*
+main()
+{
+
+ mpz_t a;
+ mpz_t b;
+ mpz_t c;
+ mpz_init(a);
+ mpz_init(b);
+ mpz_init(c);
+
+ mpz_set_si(a, 2147483647);
+ mpz_set_si(b, 2);
+ mpz_mul(c, a, b);
+
+ printf("a: %ld, b: %ld, c: %ld. \n", a->z, b->z, c->z);
+
+ return EXIT_SUCCESS;
+
+ }*/
+#endif
diff --git a/contrib/kbipack/mpz.h b/contrib/kbipack/mpz.h
new file mode 100755
index 0000000000000000000000000000000000000000..d0d508583475802942300ffb6c41b2b533e96c7a
--- /dev/null
+++ b/contrib/kbipack/mpz.h
@@ -0,0 +1,75 @@
+/*
+ This is an interface to masquerade ordinary long integers in C
+ language behind an interface mimicing arbitrary precision arithmetic
+ interger (mpz) interface in The GNU MP Bignum library.
+
+ Only the functions needed by Kbipack are available.
+
+ THIS INTERFACE IS NOT BERFORMING ANY ARBITRARY PRECISION ARITHMETIC.
+ You should always use the GMP library (http://gmplib.org/) when possible.
+
+ Some basic techniques to detect integer overflows are implemented.
+ However, they don't interfere the program flow IN ANY WAY, they just write
+ an error message to stdout.
+
+ */
+
+#ifndef __MPZ_H__
+#define __MPZ_H__
+
+#include "gmp.h"
+
+// disabled for now
+#if 0
+#include <stdlib.h>
+#include <stdio.h>
+
+typedef struct {
+ long int z;
+} mpz_struct;
+
+typedef mpz_struct mpz_t[1];
+typedef mpz_struct* mpz_ptr;
+
+// init
+void mpz_init(mpz_ptr x);
+void mpz_init_set(mpz_ptr rop, mpz_ptr op);
+void mpz_init_set_si(mpz_ptr rop, signed long int op);
+
+// set
+void mpz_set(mpz_ptr rop, mpz_ptr op);
+void mpz_set_si(mpz_ptr rop, signed long int op);
+void mpz_set_ui(mpz_ptr rop, unsigned long int op);
+
+// clear
+void mpz_clear(mpz_ptr x);
+
+
+// arithmethic
+void mpz_add(mpz_ptr rop, mpz_ptr op1, mpz_ptr op2);
+void mpz_mul(mpz_ptr rop, mpz_ptr op1, mpz_ptr op2);
+void mpz_addmul(mpz_ptr rop, mpz_ptr op1, mpz_ptr op2);
+void mpz_neg(mpz_ptr rop, mpz_ptr op);
+
+// division
+void mpz_divexact(mpz_ptr q, mpz_ptr n, mpz_ptr d);
+void mpz_cdiv_q(mpz_ptr q, mpz_ptr n, mpz_ptr d);
+void mpz_cdiv_qr(mpz_ptr q, mpz_ptr r, mpz_ptr n, mpz_ptr d);
+void mpz_tdiv_r(mpz_ptr r, mpz_ptr n, mpz_ptr d);
+
+// compare
+int mpz_cmp_si(mpz_ptr op1, signed long int op2);
+int mpz_cmpabs(mpz_ptr op1, mpz_ptr op2);
+int mpz_sgn(mpz_ptr op);
+
+// extended Euclid's algorithm
+void mpz_gcdext(mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_ptr a, mpz_ptr b);
+
+// conversion
+signed long int mpz_get_si(mpz_ptr op);
+
+// io
+size_t mpz_out_str(FILE *stream, int base, mpz_ptr op);
+#endif
+#endif
+
diff --git a/tutorial/t10.geo b/tutorial/t10.geo
index 6e8606054f0a3495490e4da8ee6056bc0c835745..b80d7b1218e2f8c90a9c2a31336dd98be5fd2d4f 100644
--- a/tutorial/t10.geo
+++ b/tutorial/t10.geo
@@ -129,15 +129,23 @@ Physical Surface(75) = {46, 63, 66, 52, 50, 48, 54, 60, 58, 56};
// Create a mesh of the model
Mesh 3;
-// Find generators of relative homology spaces of the domain modulo the four terminals.
+// Find generators of relative homology spaces of the domain modulo the
+// four terminals.
// Save the generator chains to t10_hom.msh.
HomGen("t10_hom.msh") = {{69}, {70, 71, 72, 73}};
-// Find the corresponding thin cuts.
+// Find the corresponding thin cuts,
+// generators of relative homology spaces modulo the
+// non-terminal domain surface.
// Save the cut chains to t10_hom.msh.
HomGen("t10_hom.msh") = {{69}, {75}};
-// Only find and print the ranks of the relative homology spaces (Betti numbers).
+// Find the corresponding thick cuts.
+// Save the cut chains to t10_hom.msh.
+ HomCut("t10_hom.msh") = {{69}, {70, 71, 72, 73}};
+
+// Only find and print the ranks of the relative homology spaces
+// (Betti numbers).
HomRank {{69},{70, 71, 72, 73}};
// More examples (uncomment):