diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 8336df7664c3848e2cdefa342ef9d49138a15d29..41249d859ed923345d43425b186e3f63fa400219 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -270,7 +270,7 @@ std::string GEdge::getAdditionalInfoString()
     sstream << " transfinite (" << meshAttributes.nbPointsTransfinite;
     int type = meshAttributes.typeTransfinite;
     if(std::abs(type) == 1)
-      sstream << ", progression " << sign(type) * meshAttributes.coeffTransfinite;
+      sstream << ", progression " << gmsh_sign(type) * meshAttributes.coeffTransfinite;
     else if(std::abs(type) == 2)
       sstream << ", bump " << meshAttributes.coeffTransfinite;
     sstream << ")";
diff --git a/Geo/GModelIO_GEO.cpp b/Geo/GModelIO_GEO.cpp
index 69157892f269738a674862418ad08af3506cbf5f..b293e6103e3ba6b40986db90d3067239af2370db 100644
--- a/Geo/GModelIO_GEO.cpp
+++ b/Geo/GModelIO_GEO.cpp
@@ -375,7 +375,7 @@ int GModel::importGEOInternals()
         case MSH_PHYSICAL_SURFACE: ge = getFaceByTag(tag); break;
         case MSH_PHYSICAL_VOLUME:  ge = getRegionByTag(tag); break;
       }
-      int pnum = CTX::instance()->geom.orientedPhysicals ? (sign(num) * p->Num) : p->Num;
+      int pnum = CTX::instance()->geom.orientedPhysicals ? (gmsh_sign(num) * p->Num) : p->Num;
       if(ge && std::find(ge->physicals.begin(), ge->physicals.end(), pnum) ==
           ge->physicals.end())
         ge->physicals.push_back(pnum);
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 16567d63de86a8a2f1f01ae72a2497352f839409..7af03725abc17a3efae25cef69b19a257c0ff5d4 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -4673,7 +4673,7 @@ void setSurfaceGeneratrices(Surface *s, List_T *loops)
          (i != 0 && iLoop < 0)){  // hole
         for(int j = 0; j < List_Nbr(el->Curves); j++) {
           List_Read(el->Curves, j, &ic);
-          ic *= sign(iLoop);
+          ic *= gmsh_sign(iLoop);
           if(i != 0) ic *= -1; // hole
           if(!(c = FindCurve(ic)))
             fromModel.push_back(ic);
@@ -4684,7 +4684,7 @@ void setSurfaceGeneratrices(Surface *s, List_T *loops)
       else{
         for(int j = List_Nbr(el->Curves)-1; j >= 0; j--) {
           List_Read(el->Curves, j, &ic);
-          ic *= sign(iLoop);
+          ic *= gmsh_sign(iLoop);
           if(i != 0) ic *= -1; // hole
           if(!(c = FindCurve(ic)))
             fromModel.push_back(ic);
@@ -4730,7 +4730,7 @@ void setVolumeSurfaces(Volume *v, List_T *loops)
           // create "negative" surfaces. So we just store the signs in
           // another list
           List_Add(v->Surfaces, &s);
-          int tmp = sign(is) * sign(il);
+          int tmp = gmsh_sign(is) * gmsh_sign(il);
           if(i > 0) tmp *= -1; // this is a hole
           List_Add(v->SurfacesOrientations, &tmp);
         }
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index 5185c65465ed37c32e118d87c1e1e391e5b66b74..af345dc87a3f7daf8cc316f8b918794826ec3966 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -79,7 +79,7 @@ void MPyramidN::getEdgeRep(bool curved, int num,
 
 int MPyramidN::getNumFacesRep(bool curved)
 {
-  return curved ? 6 * SQU(CTX::instance()->mesh.numSubEdges) : 6;
+  return curved ? 6 * gmsh_SQU(CTX::instance()->mesh.numSubEdges) : 6;
 }
 
 static void _myGetFaceRep(MPyramid *pyr, int num, double *x, double *y, double *z,
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 2071835ab26a31a2582672e39ca15b50e96ad092..11e9b1f8c5f58bd7f1283abd2955c0755a360190 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -99,7 +99,7 @@ std::string gmshEdge::getAdditionalInfoString()
       sstream << " transfinite (" << meshAttributes.nbPointsTransfinite;
       int type = meshAttributes.typeTransfinite;
       if(std::abs(type) == 1)
-        sstream << ", progression " << sign(type) * meshAttributes.coeffTransfinite;
+        sstream << ", progression " << gmsh_sign(type) * meshAttributes.coeffTransfinite;
       else if(std::abs(type) == 2)
         sstream << ", bump " << meshAttributes.coeffTransfinite;
       sstream << ")";
diff --git a/Geo/gmshRegion.cpp b/Geo/gmshRegion.cpp
index 7c31177df40741087696bf60f2a9532501f0aa1a..c064fb7dfffd4d0fcc8175eae2e9491f5887c508 100644
--- a/Geo/gmshRegion.cpp
+++ b/Geo/gmshRegion.cpp
@@ -32,7 +32,7 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
     GFace *f = m->getFaceByTag(abs(is));
     if(f){
       l_faces.push_back(f);
-      l_dirs.push_back(sign(is));
+      l_dirs.push_back(gmsh_sign(is));
       f->addRegion(this);
     }
     else
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 59f2a8894ba02fedfd4476fba52914b6ce32a801..752793b9b406262ed204fbe88777d43f14f66144 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -176,7 +176,7 @@ static double F_Transfinite(GEdge *ge, double t_)
 
     case 1: // Geometric progression ar^i; Sum of n terms = length = a (r^n-1)/(r-1)
       {
-        double r = (sign(type) >= 0) ? coef : 1. / coef;
+        double r = (gmsh_sign(type) >= 0) ? coef : 1. / coef;
         double a = length * (r - 1.) / (pow(r, nbpt - 1.) - 1.);
         int i = (int)(log(t * length / a * (r - 1.) + 1.) / log(r));
         val = d / (a * pow(r, (double)i));
diff --git a/Mesh/meshGFaceElliptic.cpp b/Mesh/meshGFaceElliptic.cpp
index 3589efa5e51e40cfdbad2bf5b2165b09f7dda99a..02e856468f9c7b23653441c8f757802a2faad3de 100644
--- a/Mesh/meshGFaceElliptic.cpp
+++ b/Mesh/meshGFaceElliptic.cpp
@@ -273,10 +273,10 @@ static void transfiniteSmoother(GFace* gf,
       for (int j = jStart; j<N-1; j++){
   	int jm1 = (j==0) ? N-2: j-1;
   	int jp1 = (isPeriodic ) ? (j+1)%(N-1) : j+1;
-  	double alpha = 0.25 * (SQU(uv(i,jp1).x() - uv(i,jm1).x()) +
-  			       SQU(uv(i,jp1).y() - uv(i,jm1).y())) ;
-  	double gamma = 0.25 * (SQU(uv(i+1,j).x() - uv(i-1,j).x()) +
-  			       SQU(uv(i+1,j).y() - uv(i-1,j).y()));
+  	double alpha = 0.25 * (gmsh_SQU(uv(i,jp1).x() - uv(i,jm1).x()) +
+  			       gmsh_SQU(uv(i,jp1).y() - uv(i,jm1).y())) ;
+  	double gamma = 0.25 * (gmsh_SQU(uv(i+1,j).x() - uv(i-1,j).x()) +
+  			       gmsh_SQU(uv(i+1,j).y() - uv(i-1,j).y()));
   	double beta = 0.0625 *
   	  ((uv(i+1,j).x() - uv(i-1,j).x()) * (uv(i,jp1).x() - uv(i,jm1).x()) +
   	   (uv(i+1,j).y() - uv(i-1,j).y()) * (uv(i,jp1).y() - uv(i,jm1).y()));
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 95b79f35d785e014a6032e09a1c65c6735c01cfd..e0145cea9df6bdcc91a3d80489c80334e347d361 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -36,15 +36,15 @@ double myacos(double a)
 }
 double norm2(double a[3][3]) {
   double norm2sq =
-    SQU(a[0][0])+
-    SQU(a[0][1])+
-    SQU(a[0][2])+
-    SQU(a[1][0])+
-    SQU(a[1][1])+
-    SQU(a[1][2])+
-    SQU(a[2][0])+
-    SQU(a[2][1])+
-    SQU(a[2][2]);
+    gmsh_SQU(a[0][0])+
+    gmsh_SQU(a[0][1])+
+    gmsh_SQU(a[0][2])+
+    gmsh_SQU(a[1][0])+
+    gmsh_SQU(a[1][1])+
+    gmsh_SQU(a[1][2])+
+    gmsh_SQU(a[2][0])+
+    gmsh_SQU(a[2][1])+
+    gmsh_SQU(a[2][2]);
   return sqrt(norm2sq);
 }
 
@@ -102,7 +102,7 @@ int sys2x2(double mat[2][2], double b[2], double res[2])
   double det, ud, norm;
   int i;
 
-  norm = SQU(mat[0][0]) + SQU(mat[1][1]) + SQU(mat[0][1]) + SQU(mat[1][0]);
+  norm = gmsh_SQU(mat[0][0]) + gmsh_SQU(mat[1][1]) + gmsh_SQU(mat[0][1]) + gmsh_SQU(mat[1][0]);
   det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
 
   // TOLERANCE ! WARNING WARNING
@@ -183,9 +183,9 @@ int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det)
 
   out = sys3x3(mat, b, res, det);
   norm =
-    SQU(mat[0][0]) + SQU(mat[0][1]) + SQU(mat[0][2]) +
-    SQU(mat[1][0]) + SQU(mat[1][1]) + SQU(mat[1][2]) +
-    SQU(mat[2][0]) + SQU(mat[2][1]) + SQU(mat[2][2]);
+    gmsh_SQU(mat[0][0]) + gmsh_SQU(mat[0][1]) + gmsh_SQU(mat[0][2]) +
+    gmsh_SQU(mat[1][0]) + gmsh_SQU(mat[1][1]) + gmsh_SQU(mat[1][2]) +
+    gmsh_SQU(mat[2][0]) + gmsh_SQU(mat[2][1]) + gmsh_SQU(mat[2][2]);
 
   // TOLERANCE ! WARNING WARNING
   if(norm == 0.0 || fabs(*det) / norm < 1.e-12) {
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 435a89932f460aac363dae45298e5bc7b1cb4e87..fe3d1de873eb7a86e9e4292dfc79ce4e9531fa52 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -13,9 +13,23 @@
 #include "SPoint3.h"
 #include "SVector3.h"
 
-#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b)))
-#define sign(x)      (((x)>=0)?1:-1)
-#define SQU(a)      ((a)*(a))
+//#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b)))
+
+template <class T> inline double myhypot(const T &a, const T& b)
+{
+  return sqrt((a)*(a)+(b)*(b));
+}
+
+
+template <class T> inline int gmsh_sign(const T &x)
+{
+  return (x<0)?-1:1;
+}
+
+template <class T> inline T gmsh_SQU(const T &x)
+{
+  return (x*x);
+}
 
 struct mean_plane
 {
@@ -34,6 +48,7 @@ inline double crossProd(double a[3], double b[3], int i)
   int i2 = (i+2) % 3;
   return a[i1]*b[i2] - a[i2]*b[i1];
 }
+
 inline double scalProd(double a[3], double b[3])
 {
   return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 0e03da27ff61649d95ae9d1b8c4fe01916eea698..671f6f0e607e09478918f4a159188e564a95839c 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -13,7 +13,6 @@
 #include "nodalBasis.h"
 #include <iostream>
 
-#define SQU(a)  ((a)*(a))
 
 inline double pow_int(const double &a, const int &n)
 {
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 0311486086615979d5276c4bc7a023168abf4665..afb491ba80422c41d45a5797e45a37a460e7b047 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -9826,23 +9826,23 @@ yyreduce:
           double d;
           List_Read((yyvsp[-4].l), i, &d);
           int j = (int)fabs(d);
-          for(int sign = -1; sign <= 1; sign += 2){
-            Curve *c = FindCurve(sign * j);
+          for(int sig = -1; sig <= 1; sig += 2){
+            Curve *c = FindCurve(sig * j);
             if(c){
               c->Method = MESH_TRANSFINITE;
               c->nbPointsTransfinite = ((yyvsp[-2].d) > 2) ? (int)(yyvsp[-2].d) : 2;
-              c->typeTransfinite = type * sign(d);
+              c->typeTransfinite = type * gmsh_sign(d);
               c->coeffTransfinite = coef;
             }
             else{
-              GEdge *ge = GModel::current()->getEdgeByTag(sign * j);
+              GEdge *ge = GModel::current()->getEdgeByTag(sig * j);
               if(ge){
                 ge->meshAttributes.method = MESH_TRANSFINITE;
                 ge->meshAttributes.nbPointsTransfinite = ((yyvsp[-2].d) > 2) ? (int)(yyvsp[-2].d) : 2;
-                ge->meshAttributes.typeTransfinite = type * sign(d);
+                ge->meshAttributes.typeTransfinite = type * gmsh_sign(d);
                 ge->meshAttributes.coeffTransfinite = coef;
               }
-              else if(sign > 0)
+              else if(sig > 0)
                 yymsg(0, "Unknown line %d", j);
             }
           }
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 0a1ad9c387cde09e125c8356226f0d34f9dc234d..aa03f28faa82fa324de7d6d1cc7bd2e7e09e2d03 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -4144,23 +4144,23 @@ Constraints :
           double d;
           List_Read($3, i, &d);
           int j = (int)fabs(d);
-          for(int sign = -1; sign <= 1; sign += 2){
-            Curve *c = FindCurve(sign * j);
+          for(int sig = -1; sig <= 1; sig += 2){
+            Curve *c = FindCurve(sig * j);
             if(c){
               c->Method = MESH_TRANSFINITE;
               c->nbPointsTransfinite = ($5 > 2) ? (int)$5 : 2;
-              c->typeTransfinite = type * sign(d);
+              c->typeTransfinite = type * gmsh_sign(d);
               c->coeffTransfinite = coef;
             }
             else{
-              GEdge *ge = GModel::current()->getEdgeByTag(sign * j);
+              GEdge *ge = GModel::current()->getEdgeByTag(sig * j);
               if(ge){
                 ge->meshAttributes.method = MESH_TRANSFINITE;
                 ge->meshAttributes.nbPointsTransfinite = ($5 > 2) ? (int)$5 : 2;
-                ge->meshAttributes.typeTransfinite = type * sign(d);
+                ge->meshAttributes.typeTransfinite = type * gmsh_sign(d);
                 ge->meshAttributes.coeffTransfinite = coef;
               }
-              else if(sign > 0)
+              else if(sig > 0)
                 yymsg(0, "Unknown line %d", j);
             }
           }
diff --git a/Plugin/FaultZone.cpp b/Plugin/FaultZone.cpp
index 08043fbd2a3d35dfb85449243fd87612b93395c6..a5ab23ac8f5b50bdd69f71de6610fb9ed8273b56 100644
--- a/Plugin/FaultZone.cpp
+++ b/Plugin/FaultZone.cpp
@@ -430,7 +430,7 @@ void GMSH_FaultZoneMesher::ComputeHeavisideFunction(){
           lsn = dot(vectsNor[j], vectsNor[i])*heav[i];
           assert(fabs(lsn) > tolerance || heav[i] == 0);
         }
-        heav[j] = sign(lsn);
+        heav[j] = gmsh_sign(lsn);
       }
 
       if (heav[i] != 0){
@@ -508,7 +508,7 @@ std::vector < int > GMSH_FaultZoneMesher::HeavisideFunc(MVertex* mVert, SPoint3&
       assert(dot(vectPoint, vectTan) > 0);
       lsn = 0;
     }
-    heav.push_back(sign(lsn));
+    heav.push_back(gmsh_sign(lsn));
   }
   else if (_nodesByJunctionNode.find( mVert ) != _nodesByJunctionNode.end()){
     // if it is a junction node
@@ -518,7 +518,7 @@ std::vector < int > GMSH_FaultZoneMesher::HeavisideFunc(MVertex* mVert, SPoint3&
       SVector3 vectNorm = _vectNormByFissure[fissures[i]];
       double lsn = dot(vectPoint, vectNorm);
       if (fabs(lsn) > tolerance) // tolerance seems to be ok
-        heav.push_back(sign(lsn));
+        heav.push_back(gmsh_sign(lsn));
       else
         heav.push_back(0);
     }
diff --git a/Post/shapeFunctions.h b/Post/shapeFunctions.h
index 98d68bacf474eee5ff07b0109ea38508bc4e844c..fd26320b9cabad56922d7625a8cfcf96c302a80f 100644
--- a/Post/shapeFunctions.h
+++ b/Post/shapeFunctions.h
@@ -9,7 +9,6 @@
 #include "Numeric.h"
 #include "GmshMessage.h"
 
-#define SQU(a)  ((a)*(a))
 
 class element{
 protected:
@@ -92,9 +91,9 @@ public:
         prodve(a, b, c);
         jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
       }
-      return sqrt(SQU(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
-                  SQU(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
-                  SQU(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));
+      return sqrt(gmsh_SQU(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
+                  gmsh_SQU(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
+                  gmsh_SQU(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));
     case 1:
       for(int i = 0; i < getNumNodes(); i++) {
         getGradShapeFunction(i, u, v, w, s);
@@ -114,7 +113,7 @@ public:
         jac[1][0] = b[0]; jac[1][1] = b[1]; jac[1][2] = b[2];
         jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
       }
-      return sqrt(SQU(jac[0][0])+SQU(jac[0][1])+SQU(jac[0][2]));
+      return sqrt(gmsh_SQU(jac[0][0])+gmsh_SQU(jac[0][1])+gmsh_SQU(jac[0][2]));
     default:
       jac[0][0] = jac[1][1] = jac[2][2] = 1.;
       return 1.;
@@ -246,7 +245,7 @@ public:
       double wn = uvw[2] +
         inv[0][2] * (xyz[0] - xn) + inv[1][2] * (xyz[1] - yn) + inv[2][2] * (xyz[2] - zn) ;
 
-      error = sqrt(SQU(un - uvw[0]) + SQU(vn - uvw[1]) + SQU(wn - uvw[2]));
+      error = sqrt(gmsh_SQU(un - uvw[0]) + gmsh_SQU(vn - uvw[1]) + gmsh_SQU(wn - uvw[2]));
       uvw[0] = un;
       uvw[1] = vn;
       uvw[2] = wn;
@@ -261,7 +260,7 @@ public:
     for(int i = 0; i < getNumEdges(); i++){
       int n1, n2;
       getEdge(i, n1, n2);
-      double d = sqrt(SQU(_x[n1]-_x[n2]) + SQU(_y[n1]-_y[n2]) + SQU(_z[n1]-_z[n2]));
+      double d = sqrt(gmsh_SQU(_x[n1]-_x[n2]) + gmsh_SQU(_y[n1]-_y[n2]) + gmsh_SQU(_z[n1]-_z[n2]));
       if(d > max) max = d;
     }
     return max;
@@ -936,16 +935,16 @@ public:
       switch(num) {
       case 0  : s[0] = 0.25 * ( -(1.-v) + v*w/(1.-w) ) ;
                 s[1] = 0.25 * ( -(1.-u) + u*w/(1.-w) ) ;
-                s[2] = 0.25 * ( -1.     + u*v/(1.-w) + u*v*w/SQU(1.-w) ) ; break ;
+                s[2] = 0.25 * ( -1.     + u*v/(1.-w) + u*v*w/gmsh_SQU(1.-w) ) ; break ;
       case 1  : s[0] = 0.25 * (  (1.-v) - v*w/(1.-w) ) ;
                 s[1] = 0.25 * ( -(1.+u) - u*w/(1.-w) ) ;
-                s[2] = 0.25 * ( -1.     - u*v/(1.-w) - u*v*w/SQU(1.-w) ) ; break ;
+                s[2] = 0.25 * ( -1.     - u*v/(1.-w) - u*v*w/gmsh_SQU(1.-w) ) ; break ;
       case 2  : s[0] = 0.25 * (  (1.+v) + v*w/(1.-w) ) ;
                 s[1] = 0.25 * (  (1.+u) + u*w/(1.-w) ) ;
-                s[2] = 0.25 * ( -1.     + u*v/(1.-w) + u*v*w/SQU(1.-w) ) ; break ;
+                s[2] = 0.25 * ( -1.     + u*v/(1.-w) + u*v*w/gmsh_SQU(1.-w) ) ; break ;
       case 3  : s[0] = 0.25 * ( -(1.+v) - v*w/(1.-w) ) ;
                 s[1] = 0.25 * (  (1.-u) - u*w/(1.-w) ) ;
-                s[2] = 0.25 * ( -1.     - u*v/(1.-w) - u*v*w/SQU(1.-w) ) ; break ;
+                s[2] = 0.25 * ( -1.     - u*v/(1.-w) - u*v*w/gmsh_SQU(1.-w) ) ; break ;
       case 4  : s[0] = 0. ;
                 s[1] = 0. ;
                 s[2] = 1. ; break ;