speed up normalize() and interp[12]d()

2 files changed, 15 insertions, 7 deletions

diff --git a/src/misc.c b/src/misc.c
index 3c004e4..f6e11dd 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -352,6 +352,7 @@ double interp2d(double x, double y, double *xp, double *yp, double *zp, size_t n
      * See https://en.wikipedia.org/wiki/Bilinear_interpolation. */
     size_t i, j;
     double q11, q12, q21, q22;
+    double idx, idy;
 
     if (x <= xp[0]) {
         fprintf(stderr, "x < xp[0]!\n");
@@ -363,14 +364,18 @@ double interp2d(double x, double y, double *xp, double *yp, double *zp, size_t n
         exit(1);
     }
 
-    i = (x-xp[0])/(xp[1]-xp[0]);
+    idx = 1.0/(xp[1]-xp[0]);
+
+    i = (x-xp[0])*idx;
 
     if (i > n1-2) {
         fprintf(stderr, "i > n1 - 2!\n");
         exit(1);
     }
 
-    j = (y-yp[0])/(yp[1]-yp[0]);
+    idy = 1.0/(yp[1]-yp[0]);
+
+    j = (y-yp[0])*idy;
 
     if (j > n2-2) {
         fprintf(stderr, "j > n2 - 2!\n");
@@ -382,7 +387,7 @@ double interp2d(double x, double y, double *xp, double *yp, double *zp, size_t n
     q21 = zp[(i+1)*n2+j];
     q22 = zp[(i+1)*n2+j+1];
 
-    return (q11*(xp[i+1]-x)*(yp[j+1]-y) + q21*(x-xp[i])*(yp[j+1]-y) + q12*(xp[i+1]-x)*(y-yp[j]) + q22*(x-xp[i])*(y-yp[j]))/((xp[i+1]-xp[i])*(yp[i+1]-yp[i]));
+    return (q11*(xp[i+1]-x)*(yp[j+1]-y) + q21*(x-xp[i])*(yp[j+1]-y) + q12*(xp[i+1]-x)*(y-yp[j]) + q22*(x-xp[i])*(y-yp[j]))*idx*idy;
 }
 
 double interp1d(double x, double *xp, double *yp, size_t n)
@@ -392,14 +397,17 @@ double interp1d(double x, double *xp, double *yp, size_t n)
      *
      * If x < xp[0] returns yp[0] and if x > xp[n-1] returns yp[n-1].  */
     size_t i;
+    double idx;
 
     if (x <= xp[0]) return yp[0];
 
-    i = (x-xp[0])/(xp[1]-xp[0]);
+    idx = 1.0/(xp[1]-xp[0]);
+
+    i = (x-xp[0])*idx;
 
     if (i > n-2) return yp[n-1];
 
-    return yp[i] + (yp[i+1]-yp[i])*(x-xp[i])/(xp[i+1]-xp[i]);
+    return yp[i] + (yp[i+1]-yp[i])*(x-xp[i])*idx;
 }
 
 int isclose(double a, double b, double rel_tol, double abs_tol)
diff --git a/src/vector.c b/src/vector.c
index d0e304a..061c3d6 100644
--- a/src/vector.c
+++ b/src/vector.c
@@ -4,8 +4,8 @@
 void normalize(double *a)
 {
     double c;
-    c = NORM(a);
-    DIV(a,c);
+    c = 1.0/NORM(a);
+    MUL(a,c);
 }
 
 void rotate(double *dest, double *v, double *n, double phi)