Add #include <vector> for compilation on linux

Reviewed-on: https://codeberg.org/TetGen/TetGen/pulls/10

.github/issue_template.md

@@ -0,0 +1,22 @@
+Thank you very much for your interest in TetGen! Reporting your issue helps to improve the code.
+
+### Legal
+By opening this issue, you agree to the terms of the Contributor License Agreement (CLA) stated in [CONTRIBUTING.md](https://codeberg.org/TetGen/TetGen/src/branch/main/CONTRIBUTING.md).
+
+### Kind of issue (bug, feature request...):
+
+### Issue description:
+
+### Expected behavior: 
+
+### Observed behavior: 
+
+### Environment:
+   - TeGen version:
+   - OS (e.g. MacOS 14):
+   - Architecture (e.g. x64, aarch64):
+   - Compiler vendor + version:
+
+### Minimum working example:
+
+Please attach input data which allow to reproduce the problem.

.github/pull_request_template.md

@@ -0,0 +1,12 @@
+Thank you very much for helping to improve TetGen!
+
+### Legal:
+By opening this pull request, you agree to the terms of the Contributor License Agreement (CLA) stated in [CONTRIBUTING.md](https://codeberg.org/TetGen/TetGen/src/branch/main/CONTRIBUTING.md).
+
+- [ ] I have read and agree to the Contributor License Agreement in CONTRIBUTING.md
+
+### Pull request description:
+
+
+### Test data:
+Please provide input data which help to test the feature implemented in the pull request.

CHANGELOG.md

@@ -1,7 +1,28 @@
 # TetGen: Release Notes
+## Version 1.6.0 (August 31, 2020)
+
+-   Improved the speed of the Bowyer-Watson point insertion algorithm
+    for creating Delaunay tetrahedralization.
+-   Improved the robustness of the boundary recovery algorithm (-p -Y
+    options) for creating constrained tetrahedralizations. Now the
+    default algorithm for the \`\`-p\" option (boundary recobvery) is
+    the constrained tetrahedralization algorithm. It is robust and it
+    uses less Steiner points than the constrained Delaunay
+    tetrahedralization (CDT) algorothm uses. The option to use the CDT
+    algorithm is \`\`-p -D\".
+-   A new implementation of the constrained Delaunay refinement
+    algorithm (-q option). It uses new Steiner points insertion schemes
+    to remove badly-shaped elements.
+-   Implemented a new set of mesh smoothing and mesh improvement
+    operations (-O -o) for optimizing the quality of the meshes. The
+    overall mesh quality has been improved.
+-   (Change of the -d option) To detect self-intersection in the input
+    surface mesh is now directly done in the constrained
+    tetrahedralization algorithm (the -p option).
 
 ## Version 1.5.0 (November 4, 2013)
 
+-   License switched to AGPLv3
 -   Improved the efficiency of the mesh data structure
     (tetrahedron-based).
 -   Implemented a new edge flip algorithm that does recursive

CMakeLists.txt

@@ -1,14 +1,25 @@
 # Set  the minimum  required version  of cmake  for a  project.
-cmake_minimum_required(VERSION 2.6)
+cmake_minimum_required(VERSION 3.6)
+
+include_directories(./)
+
+set(INC
+  tetgen.h
+)
+
+set(SRC
+  tetgen.cxx
+  predicates.cxx
+)
 
 # Add an executable to the project using the specified source files.
-add_executable(tetgen tetgen.cxx predicates.cxx)
+add_executable(tetgen ${INC} ${SRC})
 
 #Add a library to the project using the specified source files. 
 # In Linux/Unix, it will creates the libtet.a
-add_library(tet STATIC tetgen.cxx predicates.cxx)
+#add_library(tet STATIC tetgen.cxx predicates.cxx)
 
 #Set properties on a target. 
 #We use this here to set -DTETLIBRARY for when compiling the
 #library
-set_target_properties(tet PROPERTIES "COMPILE_DEFINITIONS" TETLIBRARY)
+#set_target_properties(tet PROPERTIES "COMPILE_DEFINITIONS" TETLIBRARY)

CONTRIBUTING.md

@@ -0,0 +1,56 @@
+# Contributor License Agreement (CLA)
+
+A "Contributor" is any person or entity that intentionally submits copyrighted material to the project,
+including but not limited to code, documentation, or other material, via a pull request, issue, or other means.
+
+By submitting a contribution Contributor agrees to the following terms:
+
+## Dual AGPLv3/MIT License Grant
+Contributor agrees that their contribution is provided under the GNU Affero General Public License v3 (AGPLv3) reproduced
+in the file [LICENSE](LICENSE) or any later version and, in addition, grants the contribution under the MIT License terms below:
+
+```text
+Permission is hereby granted, free of charge, to any person obtaining a copy  
+of this software and associated documentation files (the "Software"), to deal  
+in the Software without restriction, including without limitation the rights  
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell  
+copies of the Software, and to permit persons to whom the Software is  
+furnished to do so, subject to the following conditions:  
+
+The above copyright notice and this permission notice shall be included in all  
+copies or substantial portions of the Software.  
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR  
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,  
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE  
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER  
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  
+SOFTWARE.  
+```
+
+## Annotating contributions to source code
+
+### Trivial changes
+Trivial changes (e.g., 1-2 line fixes or typo fixes) inherit the AGPLv3 license of the project and without further annotation are also licensed under MIT as specified above.
+
+### Functions and nontrivial snippets of code 
+Functions and nontrivial snippets of code by default inherit the AGPLv3 license of the project and are also licensed under MIT as specified above. They shall be marked in a PR with an SPDX header as in this example:
+  ```
+  # SPDX-License-Identifier: AGPL-3.0-or-later OR MIT
+  # Copyright (c) [Year] [Contributor Name]
+  void frobnicate(int ntet)
+  {
+  ...
+  }
+  ```
+
+## Representations and Warranties
+Contributor represents that:
+- They are legally entitled to grant the above licenses;
+- Their contributions are original works and do not violate third-party rights;
+- They are not required to grant licenses under third-party terms (e.g., employer agreements).
+
+## No Obligation
+The Project is under no obligation to accept or use any contribution.
+

LICENSE

@@ -1,23 +1,3 @@
-TetGen License
---------------
-
-TetGen is distributed under a dual licensing scheme. You can
-redistribute it and/or modify it under the terms of the GNU Affero
-General Public License as published by the Free Software Foundation,
-either version 3 of the License, or (at your option) any later
-version. A copy of the GNU Affero General Public License is reproduced
-below.
-
-If the terms and conditions of the AGPL v.3. would prevent you from
-using TetGen, please consider the option to obtain a commercial
-license for a fee. These licenses are offered by the Weierstrass
-Institute for Applied Analysis and Stochastics (WIAS). As a rule,
-licenses are provided "as-is", unlimited in time for a one time
-fee. Please send corresponding requests to:
-tetgen@wias-berlin.de. Please do not forget to include some
-description of your company and the realm of its activities.
-
-=====================================================================
 GNU AFFERO GENERAL PUBLIC LICENSE
 
 Version 3, 19 November 2007

README

@@ -1,25 +0,0 @@
-This is TetGen version 1.5 (released on November 4, 2013)
-
-Please see the documentation of TetGen for compiling and using TetGen.
-It is available at the following link:
-
-            http://www.tetgen.org
-
-For more information on this product, contact :
-
-  Hang Si
-  Research Group of Numerical Mathematics and Scientific Computing
-  Weierstrass Institute for Applied Analysis and Stochastics
-  Mohrenstr. 39
-  10117 Berlin, Germany
-
- Phone: +49 (0) 30-20372-446   Fax: +49 (0) 30-2044975
- EMail: <si@wias-berlin.de>
- Web Site: http://www.wias-berlin.de/~si
-
-------------------- IMPORTANCE NOTICE -----------------------------
-
-BEFORE INTALLING OR USING TetGen(R) READ the 
-GENERAL LICENSE TERMS AND CONDITIONS
-
--------------------------------------------------------------------

README.md

@@ -0,0 +1,35 @@
+TetGen
+======
+
+TetGen is a program to generate tetrahedral meshes of any 3D polyhedral domains.  TetGen generates exact constrained Delaunay tetrahedralizations, boundary conforming Delaunay meshes, and Voronoi partitions.
+
+TetGen versions up to v1.6.0 have been developed headed by [Hang Si](https://github.com/sihang0592) during his affiliation with [Weierstrass Institute for Applied Analysis and Stochastics, (WIAS) Berlin](https://www.wias-berlin.de/software/tetgen).
+
+The primary development repository is [https://codeberg.org/TetGen/TetGen](https://codeberg.org/TetGen/TetGen). Please
+open any issues in this repository. The github repository [https://github.com/TetGen/Tetgen](https://github.com/TetGen/Tetgen)
+is a mirror of the codeberg repository.
+
+This repository provides the source code of the most recent versions of TetGen.
+The following releases are provided:
+- [v1.4.3, 2011](https://codeberg.org/TetGen/TetGen/archive/v1.4.3.tar.gz): [MIT license with noncommercial clause](https://raw.githubusercontent.com/TetGen/TetGen/refs/tags/v1.4.3/LICENSE).
+- [v1.5.0, 2013](https://codeberg.org/TetGen/TetGen/archive/v1.5.0.tar.gz): [AGPLv3 license](https://www.gnu.org/licenses/agpl-3.0.html).
+- [v1.5.1, 2018](https://codeberg.org/TetGen/TetGen/archive/v1.5.1.tar.gz): [AGPLv3 license](https://www.gnu.org/licenses/agpl-3.0.html), recommended as most stable version.
+- [v1.6.0, 2020](https://codeberg.org/TetGen/TetGen/archive/v1.6.0.tar.gz): [AGPLv3 license](https://www.gnu.org/licenses/agpl-3.0.html), most recent version with some rough edges.
+
+It is planned to resume active development of TetGen in cooperation between the main developer Hang Si and  the Weierstrass Institute. Future versions of TetGen will be published via this repository.
+
+## Documentation
+- Manuals are available in https://codeberg.org/TetGen/Manuals.
+- A technical paper about TetGen is available at [Hang Si, "TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator". ACM Trans. on Mathematical Software. 41 (2), 2015](http://doi.acm.org/10.1145/2629697)
+  
+  
+## Commercial licensing
+
+TetGen versions developed at WIAS are distributed under a dual licensing scheme.  As an alternative to the use of TetGen under the AGPLv3 license, consider the option to obtain a commercial license for a fee.  These licenses are offered by the Weierstrass Institute for Applied Analysis and Stochastics (WIAS). As a rule, licenses are provided "as-is", unlimited in time for a one time fee.  Please send corresponding requests to: `tetgen at wias-berlin.de`.  Please do not forget to include some description of your company and the realm of its activities.
+
+Details about the extension of this dual licensing scheme for future versions of TetGen will be announced in due time.
+
+## Contributing
+
+By submitting a contribution (including but not limited to code, documentation, or other materials) via a pull request, issue, or other means, a contributor agrees to the terms of the contributor license agreement stated in [CONTRIBUTING.md](CONTRIBUTING.md). It is planned to revise this agreement, once the extension of the dual licensing scheme has been set up.
+

example.poly

@@ -1,84 +0,0 @@
-28  3  0  1
-1  0  0  0  1
-2  2  0  0  1
-3  2  2  0  1
-4  0  2  0  1
-5  0  0  4  9
-6  2  0  4  9
-7  2  2  3  9
-8  0  2  3  9
-9  0  0  5  2
-10  2  0  5  2
-11  2  2  5  2
-12  0  2  5  2
-13  0.25  0.25  0.5  4
-14  1.75  0.25  0.5  4
-15  1.75  1.5  0.5  4
-16  0.25  1.5  0.5  4
-17  0.25  0.25  1  4
-18  1.75  0.25  1  4
-19  1.75  1.5  1  4
-20  0.25  1.5  1  4
-21  0.25  0  2  4
-22  1.75  0  2  4
-23  1.75  1.5  2  4
-24  0.25  1.5  2  4
-25  0.25  0  2.5  4
-26  1.75  0  2.5  4
-27  1.75  1.5  2.5  4
-28  0.25  1.5  2.5  4
-23  1
-1  0  1  # 1
-4    1  2  3  4
-1  0  9  # 2
-4    5  6  7  8
-2  1  3  # 3
-4    1  2  6  5
-4    21  22  26  25
-1  1  0  2.25
-1  0  3  # 4
-4    2  3  7  6
-1  0  3  # 5
-4    3  4  8  7
-1  0  3  # 6
-4    4  1  5  8
-1  0  2  # 7
-4    9  10  11  12
-1  0  3  # 8
-4    9  10  6  5
-1  0  3  # 9
-4    10  11  7  6
-1  0  3  # 10
-4    11  12  8  7
-1  0  3  # 11
-4    12  9  5  8
-1  0  4  # 12
-4    13  14  15  16
-1  0  4  # 13
-4    17  18  19  20
-1  0  4  # 14
-4    13  14  18  17
-1  0  4  # 15
-4    14  15  19  18
-1  0  4  # 16
-4    15  16  20  19
-1  0  4  # 17
-4    16  13  17  20
-1  0  4  # 18
-4    21  22  23  24
-1  0  4  # 19
-4    25  26  27  28
-1  0  4  # 20
-4    21  22  26  25
-1  0  4  # 21
-4    22  23  27  26
-1  0  4  # 22
-4    23  24  28  27
-1  0  4  # 23
-4    24  21  25  28
-2
-1  1  0.4  2.25
-2  1  0.4  0.75
-2
-1  1  0.25  0.1  10  0.001
-2  1  0.5  4  20  0.01

predicates.cxx

@@ -399,7 +399,7 @@ static REAL ispstaticfilter;
 //          http://www.math.utah.edu/~beebe/software/ieee/
 // The original program was "fpinfo2.c".
 
-double fppow2(int n)
+static double fppow2(int n)
 {
   double x, power;
   x = (n < 0) ? ((double)1.0/(double)2.0) : (double)2.0;
@@ -412,12 +412,12 @@ double fppow2(int n)
 
 #ifdef SINGLE
 
-float fstore(float x)
+static float fstore(float x)
 {
   return (x);
 }
 
-int test_float(int verbose)
+static int test_float(int verbose)
 {
   float x;
   int pass = 1;
@@ -467,12 +467,12 @@ int test_float(int verbose)
 
 # else
 
-double dstore(double x)
+static double dstore(double x)
 {
   return (x);
 }
 
-int test_double(int verbose)
+static int test_double(int verbose)
 {
   double x;
   int pass = 1;
@@ -629,10 +629,6 @@ void exactinit(int verbose, int noexact, int nofilter, REAL maxx, REAL maxy,
   // Added by H. Si, 2012-08-23.
 
   // Sort maxx < maxy < maxz. Re-use 'half' for swapping.
-  assert(maxx > 0);
-  assert(maxy > 0);
-  assert(maxz > 0);
-
   if (maxx > maxz) {
     half = maxx; maxx = maxz; maxz = half;
   }
@@ -885,7 +881,7 @@ int expansion_sum_zeroelim2(int elen, REAL *e, int flen, REAL *f, REAL *h)
 /*                                                                           */
 /*****************************************************************************/
 
-int fast_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
+static int fast_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
 /* h cannot be e or f. */
 {
   REAL Q;
@@ -957,7 +953,7 @@ int fast_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
 /*  properties.                                                              */
 /*                                                                           */
 /*****************************************************************************/
-
+static
 int fast_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, REAL *h)
 /* h cannot be e or f. */
 {
@@ -1169,7 +1165,7 @@ int linear_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f,
 /*  will h.)                                                                 */
 /*                                                                           */
 /*****************************************************************************/
-
+static
 int scale_expansion(int elen, REAL *e, REAL b, REAL *h)
 /* e and h cannot be the same. */
 {
@@ -1215,7 +1211,7 @@ int scale_expansion(int elen, REAL *e, REAL b, REAL *h)
 /*  will h.)                                                                 */
 /*                                                                           */
 /*****************************************************************************/
-
+static
 int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
 /* e and h cannot be the same. */
 {
@@ -1267,7 +1263,7 @@ int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
 /*  nonadjacent expansion.                                                   */
 /*                                                                           */
 /*****************************************************************************/
-
+static
 int compress(int elen, REAL *e, REAL *h)
 /* e and h may be the same. */
 {
@@ -1310,7 +1306,7 @@ int compress(int elen, REAL *e, REAL *h)
 /*  See either version of my paper for details.                              */
 /*                                                                           */
 /*****************************************************************************/
-
+static
 REAL estimate(int elen, REAL *e)
 {
   REAL Q;
@@ -1348,7 +1344,7 @@ REAL estimate(int elen, REAL *e)
 /*  nearly so.                                                               */
 /*                                                                           */
 /*****************************************************************************/
-
+static
 REAL orient2dfast(REAL *pa, REAL *pb, REAL *pc)
 {
   REAL acx, bcx, acy, bcy;
@@ -1360,6 +1356,7 @@ REAL orient2dfast(REAL *pa, REAL *pb, REAL *pc)
   return acx * bcy - acy * bcx;
 }
 
+//static
 REAL orient2dexact(REAL *pa, REAL *pb, REAL *pc)
 {
   INEXACT REAL axby1, axcy1, bxcy1, bxay1, cxay1, cxby1;
@@ -2193,6 +2190,9 @@ REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
 
 REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
 {
+//#ifdef USING_GMP
+//  return orient3d_gmp(pa, pb, pc, pd);
+//#endif
   REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
   REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
   REAL det;
@@ -4185,7 +4185,7 @@ REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
 /*  See my Robust Predicates paper for details.                              */
 /*                                                                           */
 /*****************************************************************************/
-
+//static
 REAL orient4dexact(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
                    REAL aheight, REAL bheight, REAL cheight, REAL dheight, 
                    REAL eheight)
@@ -4401,6 +4401,7 @@ REAL orient4dexact(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
   return deter[deterlen - 1];
 }
 
+static
 REAL orient4dadapt(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
                    REAL aheight, REAL bheight, REAL cheight, REAL dheight, 
                    REAL eheight, REAL permanent)
@@ -4704,3 +4705,625 @@ REAL orient4d(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
 
 
 
+
+//==============================================================================
+
+static REAL det3x3(REAL adx, REAL ady, REAL adz,
+                   REAL bdx, REAL bdy, REAL bdz,
+                   REAL cdx, REAL cdy, REAL cdz)
+{
+  return adx * (bdy * cdz - bdz * cdy)
+       + bdx * (cdy * adz - cdz * ady)
+       + cdx * (ady * bdz - adz * bdy);
+}
+
+void tetgenmesh::pre_calculate_insphere(point pa, point pb, point pc, point pd,
+                                        REAL* dets)
+{
+  if (pd != dummypoint) {
+    REAL ba_x = pb[0] - pa[0];
+    REAL ba_y = pb[1] - pa[1];
+    REAL ba_z = pb[2] - pa[2];
+    REAL ba_norm = ba_x*ba_x + ba_y*ba_y + ba_z*ba_z;
+    
+    REAL ca_x = pc[0] - pa[0];
+    REAL ca_y = pc[1] - pa[1];
+    REAL ca_z = pc[2] - pa[2];
+    REAL ca_norm = ca_x*ca_x + ca_y*ca_y + ca_z*ca_z;
+    
+    REAL da_x = pd[0] - pa[0];
+    REAL da_y = pd[1] - pa[1];
+    REAL da_z = pd[2] - pa[2];
+    REAL da_norm = da_x*da_x + da_y*da_y + da_z*da_z;
+
+    dets[0] = det3x3(ba_y, ba_z, ba_norm,
+                     ca_y, ca_z, ca_norm,
+                     da_y, da_z, da_norm);
+    
+    dets[1] = det3x3(ba_x, ba_z, ba_norm,
+                     ca_x, ca_z, ca_norm,
+                     da_x, da_z, da_norm);
+
+    dets[2] = det3x3(ba_x, ba_y, ba_norm,
+                     ca_x, ca_y, ca_norm,
+                     da_x, da_y, da_norm);
+    
+    dets[3] = det3x3(ba_x, ba_y, ba_z,
+                     ca_x, ca_y, ca_z,
+                     da_x, da_y, da_z);
+  } else {
+    double ab[4],ac[4];
+    double* a = pa; // mesh->vertices[Node[0]].coord;
+    double* b = pb; //mesh->vertices[Node[1]].coord;
+    double* c = pc; //mesh->vertices[Node[2]].coord;
+    unsigned i;
+    for (i=0; i<3; i++)
+    {
+      ab[i]=b[i]-a[i]; //AB
+      ac[i]=c[i]-a[i]; //AC
+    }
+    dets[0] = ac[1]*ab[2] - ac[2]*ab[1];
+    dets[1] = ac[2]*ab[0] - ac[0]*ab[2];
+    dets[2] = ac[0]*ab[1] - ac[1]*ab[0];
+    dets[3] = dets[0]*dets[0] + dets[1]*dets[1] + dets[2]*dets[2];
+  }
+}
+
+REAL tetgenmesh::insphere_use_subdets(tetrahedron *tet, REAL* pe)
+{
+  REAL *dets = get_polar(tet);
+
+  if (dets[3] == 0) { // Only calculate once.
+    point *pts = (point *) tet;
+    pre_calculate_insphere(pts[4], pts[5], pts[6], pts[7], dets);
+  }
+
+  point pa = (point) tet[4];
+  
+  if (((point) tet[7]) == dummypoint) {
+    double aex = pe[0] - pa[0];
+    double aey = pe[1] - pa[1];
+    double aez = pe[2] - pa[2];
+    double det = aex*dets[0]+aey*dets[1]+aez*dets[2];
+    if(fabs(det) > o3dstaticfilter) return det;
+    point *pts = (point *) tet;
+    det = orient3d(pts[4],pts[4],pts[6],pe);
+    return det;
+  }
+  
+  REAL ea_x = pe[0] - pa[0];
+  REAL ea_y = pe[1] - pa[1];
+  REAL ea_z = pe[2] - pa[2];
+  REAL ea_norm = ea_x * ea_x + ea_y * ea_y + ea_z * ea_z;
+
+  REAL det = -ea_x * dets[0]
+           +  ea_y * dets[1]
+           -  ea_z * dets[2]
+           +  ea_norm * dets[3];
+
+  if (fabs(det) < ispstaticfilter) {
+    point *pts = (point *) tet;
+    det = insphere_s(pts[4], pts[5], pts[6], pts[7], pe);
+  }
+  return det;
+}
+
+#ifdef USING_GMP
+
+#include <gmp.h>
+
+//============================================================================//
+//                                                                            //
+// orient3d_gmp()                                                             //
+//                                                                            //
+//============================================================================//
+
+REAL orient3d_gmp(REAL *aa, REAL *bb, REAL *cc, REAL *dd)
+{
+  mpz_t pa[3], pb[3], pc[3], pd[3];
+  
+  for (int i = 0; i < 3; i++) {
+    mpz_init_set_d(pa[i], aa[i]*1.e+15);
+    mpz_init_set_d(pb[i], bb[i]*1.e+15);
+    mpz_init_set_d(pc[i], cc[i]*1.e+15);
+    mpz_init_set_d(pd[i], dd[i]*1.e+15);
+  }
+  
+  mpz_t adx, bdx, cdx;
+  mpz_t ady, bdy, cdy;
+  mpz_t adz, bdz, cdz;
+  
+  //adx = pa[0] - pd[0];
+  mpz_init(adx);
+  mpz_sub(adx, pa[0], pd[0]);
+  //bdx = pb[0] - pd[0];
+  mpz_init(bdx);
+  mpz_sub(bdx, pb[0], pd[0]);
+  //cdx = pc[0] - pd[0];
+  mpz_init(cdx);
+  mpz_sub(cdx, pc[0], pd[0]);
+  //ady = pa[1] - pd[1];
+  mpz_init(ady);
+  mpz_sub(ady, pa[1], pd[1]);
+  //bdy = pb[1] - pd[1];
+  mpz_init(bdy);
+  mpz_sub(bdy, pb[1], pd[1]);
+  //cdy = pc[1] - pd[1];
+  mpz_init(cdy);
+  mpz_sub(cdy, pc[1], pd[1]);
+  //adz = pa[2] - pd[2];
+  mpz_init(adz);
+  mpz_sub(adz, pa[2], pd[2]);
+  //bdz = pb[2] - pd[2];
+  mpz_init(bdz);
+  mpz_sub(bdz, pb[2], pd[2]);
+  //cdz = pc[2] - pd[2];
+  mpz_init(cdz);
+  mpz_sub(cdz, pc[2], pd[2]);
+
+  /*
+  return adx * (bdy * cdz - bdz * cdy)
+       + bdx * (cdy * adz - cdz * ady)
+       + cdx * (ady * bdz - adz * bdy);
+   */
+  mpz_t bdy_cdz;
+  mpz_init(bdy_cdz);
+  mpz_mul(bdy_cdz, bdy, cdz);
+  mpz_t bdz_cdy;
+  mpz_init(bdz_cdy);
+  mpz_mul(bdz_cdy, bdz, cdy);
+  
+  mpz_t cdy_adz;
+  mpz_init(cdy_adz);
+  mpz_mul(cdy_adz, cdy, adz);
+  mpz_t cdz_ady;
+  mpz_init(cdz_ady);
+  mpz_mul(cdz_ady, cdz, ady);
+  
+  mpz_t ady_bdz;
+  mpz_init(ady_bdz);
+  mpz_mul(ady_bdz, ady, bdz);
+  mpz_t adz_bdy;
+  mpz_init(adz_bdy);
+  mpz_mul(adz_bdy, adz, bdy);
+  
+  mpz_t bdy_cdz_bdz_cdy;
+  mpz_init(bdy_cdz_bdz_cdy);
+  mpz_sub(bdy_cdz_bdz_cdy, bdy_cdz, bdz_cdy);
+  
+  mpz_t cdy_adz_cdz_ady;
+  mpz_init(cdy_adz_cdz_ady);
+  mpz_sub(cdy_adz_cdz_ady, cdy_adz, cdz_ady);
+  
+  mpz_t ady_bdz_adz_bdy;
+  mpz_init(ady_bdz_adz_bdy);
+  mpz_sub(ady_bdz_adz_bdy, ady_bdz, adz_bdy);
+  
+  mpz_t adx_;
+  mpz_init(adx_);
+  mpz_mul(adx_, adx, bdy_cdz_bdz_cdy);
+  mpz_t bdx_;
+  mpz_init(bdx_);
+  mpz_mul(bdx_, bdx, cdy_adz_cdz_ady);
+  mpz_t cdx_;
+  mpz_init(cdx_);
+  mpz_mul(cdx_, cdx, ady_bdz_adz_bdy);
+  
+  mpz_t det;
+  mpz_init(det);
+  mpz_add(det, adx_, bdx_);
+  mpz_add(det, det, cdx_);
+  
+  bool debug_flag = false;
+  if (debug_flag) { // Debug only
+    char str[1024];
+    mpz_get_str(str, 10, det);
+    printf("\ndet_str = %s\n", str);
+    double detd = mpz_get_d(det);
+    printf ("\ndetd = %.17g\n", detd);
+  }
+  
+  int sign = mpz_sgn(det);
+  
+  for (int i = 0; i < 3; i++) {
+    mpz_clear(pa[i]);
+    mpz_clear(pb[i]);
+    mpz_clear(pc[i]);
+    mpz_clear(pd[i]);
+  }
+  
+  mpz_clear(adx); mpz_clear(bdx); mpz_clear(cdx);
+  mpz_clear(ady); mpz_clear(bdy); mpz_clear(cdy);
+  mpz_clear(adz); mpz_clear(bdz); mpz_clear(cdz);
+  
+  mpz_clear(bdy_cdz);
+  mpz_clear(bdz_cdy);
+  mpz_clear(cdy_adz);
+  mpz_clear(cdz_ady);
+  mpz_clear(ady_bdz);
+  mpz_clear(adz_bdy);
+  mpz_clear(bdy_cdz_bdz_cdy);
+  mpz_clear(cdy_adz_cdz_ady);
+  mpz_clear(ady_bdz_adz_bdy);
+  mpz_clear(adx_);
+  mpz_clear(bdx_);
+  mpz_clear(cdx_);
+  mpz_clear(det);
+  
+  return double(sign);
+}
+
+//============================================================================//
+//                                                                            //
+// insphere_gmp()                                                             //
+//                                                                            //
+//============================================================================//
+
+REAL insphere_gmp(REAL *aa, REAL *bb, REAL *cc, REAL *dd, REAL *ee)
+{
+  mpz_t pa[3], pb[3], pc[3], pd[3], pe[3];
+
+  for (int i = 0; i < 3; i++) {
+    mpz_init_set_d(pa[i], aa[i]*1.e+15);
+    mpz_init_set_d(pb[i], bb[i]*1.e+15);
+    mpz_init_set_d(pc[i], cc[i]*1.e+15);
+    mpz_init_set_d(pd[i], dd[i]*1.e+15);
+    mpz_init_set_d(pe[i], ee[i]*1.e+15);
+  }
+
+  mpz_t aex, bex, cex, dex;
+  mpz_t aey, bey, cey, dey;
+  mpz_t aez, bez, cez, dez;
+
+  //aex = pa[0] - pe[0];
+  mpz_init(aex);
+  mpz_sub(aex, pa[0], pe[0]);
+  //bex = pb[0] - pe[0];
+  mpz_init(bex);
+  mpz_sub(bex, pb[0], pe[0]);
+  //cex = pc[0] - pe[0];
+  mpz_init(cex);
+  mpz_sub(cex, pc[0], pe[0]);
+  //dex = pd[0] - pe[0];
+  mpz_init(dex);
+  mpz_sub(dex, pd[0], pe[0]);
+  //aey = pa[1] - pe[1];
+  mpz_init(aey);
+  mpz_sub(aey, pa[1], pe[1]);
+  //bey = pb[1] - pe[1];
+  mpz_init(bey);
+  mpz_sub(bey, pb[1], pe[1]);
+  //cey = pc[1] - pe[1];
+  mpz_init(cey);
+  mpz_sub(cey, pc[1], pe[1]);
+  //dey = pd[1] - pe[1];
+  mpz_init(dey);
+  mpz_sub(dey, pd[1], pe[1]);
+  //aez = pa[2] - pe[2];
+  mpz_init(aez);
+  mpz_sub(aez, pa[2], pe[2]);
+  //bez = pb[2] - pe[2];
+  mpz_init(bez);
+  mpz_sub(bez, pb[2], pe[2]);
+  //cez = pc[2] - pe[2];
+  mpz_init(cez);
+  mpz_sub(cez, pc[2], pe[2]);
+  //dez = pd[2] - pe[2];
+  mpz_init(dez);
+  mpz_sub(dez, pd[2], pe[2]);
+
+  mpz_t aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey;
+  mpz_t aexcey, cexaey, bexdey, dexbey;
+  mpz_t alift, blift, clift, dlift;
+  mpz_t ab, bc, cd, da, ac, bd;
+  mpz_t abc, bcd, cda, dab;
+  mpz_t det;
+
+  //aexbey = aex * bey;
+  mpz_init(aexbey);
+  mpz_mul(aexbey, aex, bey);
+  
+  //bexaey = bex * aey;
+  mpz_init(bexaey);
+  mpz_mul(bexaey, bex, aey);
+  
+  //ab = aexbey - bexaey;
+  mpz_init(ab);
+  mpz_sub(ab, aexbey, bexaey);
+  
+  //bexcey = bex * cey;
+  mpz_init(bexcey);
+  mpz_mul(bexcey, bex, cey);
+  
+  //cexbey = cex * bey;
+  mpz_init(cexbey);
+  mpz_mul(cexbey, cex, bey);
+  
+  //bc = bexcey - cexbey;
+  mpz_init(bc);
+  mpz_sub(bc, bexcey, cexbey);
+  
+  //cexdey = cex * dey;
+  mpz_init(cexdey);
+  mpz_mul(cexdey, cex, dey);
+  
+  //dexcey = dex * cey;
+  mpz_init(dexcey);
+  mpz_mul(dexcey, dex, cey);
+  
+  //cd = cexdey - dexcey;
+  mpz_init(cd);
+  mpz_sub(cd, cexdey, dexcey);
+  
+  //dexaey = dex * aey;
+  mpz_init(dexaey);
+  mpz_mul(dexaey, dex, aey);
+  
+  //aexdey = aex * dey;
+  mpz_init(aexdey);
+  mpz_mul(aexdey, aex, dey);
+  
+  //da = dexaey - aexdey;
+  mpz_init(da);
+  mpz_sub(da, dexaey, aexdey);
+  
+  //aexcey = aex * cey;
+  mpz_init(aexcey);
+  mpz_mul(aexcey, aex, cey);
+  
+  //cexaey = cex * aey;
+  mpz_init(cexaey);
+  mpz_mul(cexaey, cex, aey);
+  
+  //ac = aexcey - cexaey;
+  mpz_init(ac);
+  mpz_sub(ac, aexcey, cexaey);
+  
+  //bexdey = bex * dey;
+  mpz_init(bexdey);
+  mpz_mul(bexdey, bex, dey);
+  
+  //dexbey = dex * bey;
+  mpz_init(dexbey);
+  mpz_mul(dexbey, dex, bey);
+  
+  //bd = bexdey - dexbey;
+  mpz_init(bd);
+  mpz_sub(bd, bexdey, dexbey);
+  
+  //abc = aez * bc - bez * ac + cez * ab;
+  mpz_t cez_ab;
+  mpz_init(cez_ab);
+  mpz_mul(cez_ab, cez, ab);
+  mpz_t bez_ac;
+  mpz_init(bez_ac);
+  mpz_mul(bez_ac, bez, ac);
+  mpz_t aez_bc;
+  mpz_init(aez_bc);
+  mpz_mul(aez_bc, aez, bc);
+  mpz_init_set(abc, aez_bc);
+  mpz_sub(abc, abc, bez_ac);
+  mpz_add(abc, abc, cez_ab);
+  
+  //bcd = bez * cd - cez * bd + dez * bc;
+  mpz_t bez_cd;
+  mpz_init(bez_cd);
+  mpz_mul(bez_cd, bez, cd);
+  mpz_t cez_bd;
+  mpz_init(cez_bd);
+  mpz_mul(cez_bd, cez, bd);
+  mpz_t dez_bc;
+  mpz_init(dez_bc);
+  mpz_mul(dez_bc, dez, bc);
+  mpz_init_set(bcd, bez_cd);
+  mpz_sub(bcd, bcd, cez_bd);
+  mpz_add(bcd, bcd, dez_bc);
+  
+  //cda = cez * da + dez * ac + aez * cd;
+  mpz_t cez_da;
+  mpz_init(cez_da);
+  mpz_mul(cez_da, cez, da);
+  mpz_t dez_ac;
+  mpz_init(dez_ac);
+  mpz_mul(dez_ac, dez, ac);
+  mpz_t aez_cd;
+  mpz_init(aez_cd);
+  mpz_mul(aez_cd, aez, cd);
+  mpz_init_set(cda, cez_da);
+  mpz_add(cda, cda, dez_ac);
+  mpz_add(cda, cda, aez_cd);
+  
+  //dab = dez * ab + aez * bd + bez * da;
+  mpz_t dez_ab;
+  mpz_init(dez_ab);
+  mpz_mul(dez_ab, dez, ab);
+  mpz_t aez_bd;
+  mpz_init(aez_bd);
+  mpz_mul(aez_bd, aez, bd);
+  mpz_t bez_da;
+  mpz_init(bez_da);
+  mpz_mul(bez_da, bez, da);
+  mpz_init_set(dab, dez_ab);
+  mpz_add(dab, dab, aez_bd);
+  mpz_add(dab, dab, bez_da);
+  
+  //alift = aex * aex + aey * aey + aez * aez;
+  mpz_t aex_aex;
+  mpz_init(aex_aex);
+  mpz_mul(aex_aex, aex, aex);
+  mpz_t aey_aey;
+  mpz_init(aey_aey);
+  mpz_mul(aey_aey, aey, aey);
+  mpz_t aez_aez;
+  mpz_init(aez_aez);
+  mpz_mul(aez_aez, aez, aez);
+  mpz_init_set(alift, aex_aex);
+  mpz_add(alift, alift, aey_aey);
+  mpz_add(alift, alift, aez_aez);
+  
+  //blift = bex * bex + bey * bey + bez * bez;
+  mpz_t bex_bex;
+  mpz_init(bex_bex);
+  mpz_mul(bex_bex, bex, bex);
+  mpz_t bey_bey;
+  mpz_init(bey_bey);
+  mpz_mul(bey_bey, bey, bey);
+  mpz_t bez_bez;
+  mpz_init(bez_bez);
+  mpz_mul(bez_bez, bez, bez);
+  mpz_init_set(blift, bex_bex);
+  mpz_add(blift, blift, bey_bey);
+  mpz_add(blift, blift, bez_bez);
+  
+  //clift = cex * cex + cey * cey + cez * cez;
+  mpz_t cex_cex;
+  mpz_init(cex_cex);
+  mpz_mul(cex_cex, cex, cex);
+  mpz_t cey_cey;
+  mpz_init(cey_cey);
+  mpz_mul(cey_cey, cey, cey);
+  mpz_t cez_cez;
+  mpz_init(cez_cez);
+  mpz_mul(cez_cez, cez, cez);
+  mpz_init_set(clift, cex_cex);
+  mpz_add(clift, clift, cey_cey);
+  mpz_add(clift, clift, cez_cez);
+  
+  //dlift = dex * dex + dey * dey + dez * dez;
+  mpz_t dex_dex;
+  mpz_init(dex_dex);
+  mpz_mul(dex_dex, dex, dex);
+  mpz_t dey_dey;
+  mpz_init(dey_dey);
+  mpz_mul(dey_dey, dey, dey);
+  mpz_t dez_dez;
+  mpz_init(dez_dez);
+  mpz_mul(dez_dez, dez, dez);
+  mpz_init_set(dlift, dex_dex);
+  mpz_add(dlift, dlift, dey_dey);
+  mpz_add(dlift, dlift, dez_dez);
+  
+  //det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
+  mpz_t dlift_abc;
+  mpz_init(dlift_abc);
+  mpz_mul(dlift_abc, dlift, abc);
+  mpz_t clift_dab;
+  mpz_init(clift_dab);
+  mpz_mul(clift_dab, clift, dab);
+  mpz_t blift_cda;
+  mpz_init(blift_cda);
+  mpz_mul(blift_cda, blift, cda);
+  mpz_t alift_bcd;
+  mpz_init(alift_bcd);
+  mpz_mul(alift_bcd, alift, bcd);
+  mpz_init_set(det, dlift_abc);
+  mpz_sub(det, det, clift_dab);
+  mpz_add(det, det, blift_cda);
+  mpz_sub(det, det, alift_bcd);
+  
+  bool debug_flag = false;
+  if (debug_flag) { // Debug only
+    char str[1024];
+    mpz_get_str(str, 10, det);
+    printf("\ndet_str = %s\n", str);
+    double detd = mpz_get_d(det);
+    printf ("\ndetd = %.17g\n", detd);
+  }
+  
+  int sign = mpz_sgn(det);
+  
+  // Clear memory
+  for (int i = 0; i < 3; i++) {
+    mpz_clear(pa[i]);
+    mpz_clear(pb[i]);
+    mpz_clear(pc[i]);
+    mpz_clear(pd[i]);
+    mpz_clear(pe[i]);
+  }
+  
+  mpz_clear(aex);
+  mpz_clear(bex);
+  mpz_clear(cex);
+  mpz_clear(dex);
+  mpz_clear(aey);
+  mpz_clear(bey);
+  mpz_clear(cey);
+  mpz_clear(dey);
+  mpz_clear(aez);
+  mpz_clear(bez);
+  mpz_clear(cez);
+  mpz_clear(dez);
+  
+  mpz_clear(aexbey);
+  mpz_clear(bexaey);
+  mpz_clear(bexcey);
+  mpz_clear(cexbey);
+  mpz_clear(cexdey);
+  mpz_clear(dexcey);
+  mpz_clear(dexaey);
+  mpz_clear(aexdey);
+  mpz_clear(aexcey);
+  mpz_clear(cexaey);
+  mpz_clear(bexdey);
+  mpz_clear(dexbey);
+  mpz_clear(alift);
+  mpz_clear(blift);
+  mpz_clear(clift);
+  mpz_clear(dlift);
+  mpz_clear(ab);
+  mpz_clear(bc);
+  mpz_clear(cd);
+  mpz_clear(da);
+  mpz_clear(ac);
+  mpz_clear(bd);
+  mpz_clear(abc);
+  mpz_clear(bcd);
+  mpz_clear(cda);
+  mpz_clear(dab);
+  mpz_clear(det);
+  
+  mpz_clear(cez_ab);
+  mpz_clear(bez_ac);
+  mpz_clear(aez_bc);
+  
+  mpz_clear(bez_cd);
+  mpz_clear(cez_bd);
+  mpz_clear(dez_bc);
+  
+  mpz_clear(cez_da);
+  mpz_clear(dez_ac);
+  mpz_clear(aez_cd);
+  
+  mpz_clear(dez_ab);
+  mpz_clear(aez_bd);
+  mpz_clear(bez_da);
+  
+  mpz_clear(aex_aex);
+  mpz_clear(aey_aey);
+  mpz_clear(aez_aez);
+  
+  mpz_clear(bex_bex);
+  mpz_clear(bey_bey);
+  mpz_clear(bez_bez);
+  
+  mpz_clear(cex_cex);
+  mpz_clear(cey_cey);
+  mpz_clear(cez_cez);
+  
+  mpz_clear(dex_dex);
+  mpz_clear(dey_dey);
+  mpz_clear(dez_dez);
+  
+  mpz_clear(dlift_abc);
+  mpz_clear(clift_dab);
+  mpz_clear(blift_cda);
+  mpz_clear(alift_bcd);
+  
+  return double(sign);
+}
+
+#endif
+
+//==============================================================================
+
+