Draft version 1.6.1 by Hang

2026-03-25 21:57:21 +08:00
parent be8813bec1
commit bd6e268ce9
5 changed files with 7023 additions and 1601 deletions
--- a/predicates.cxx
+++ b/predicates.cxx
@@ -2190,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;
@@ -4703,7 +4706,623 @@ 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

 //==============================================================================