regress: add linear regression with two independent variables

regress.c

@@ -653,3 +653,57 @@ RGR_FindBestRobustRegression
 }
 
 /* ================================================== */
+/* This routine performs linear regression with two independent variables.
+   It returns non-zero status if there were enough data points and there
+   was a solution. */
+
+int
+RGR_MultipleRegress
+(double *x1,                    /* first independent variable */
+ double *x2,                    /* second independent variable */
+ double *y,                     /* measured data */
+
+ int n,                         /* number of data points */
+
+ /* The results */
+ double *b2                     /* estimated second slope */
+                                /* other values are not needed yet */
+)
+{
+  double Sx1, Sx2, Sx1x1, Sx1x2, Sx2x2, Sx1y, Sx2y, Sy;
+  double U, V, V1, V2, V3;
+  int i;
+
+  if (n < 4)
+    return 0;
+
+  Sx1 = Sx2 = Sx1x1 = Sx1x2 = Sx2x2 = Sx1y = Sx2y = Sy = 0.0;
+
+  for (i = 0; i < n; i++) {
+    Sx1 += x1[i];
+    Sx2 += x2[i];
+    Sx1x1 += x1[i] * x1[i];
+    Sx1x2 += x1[i] * x2[i];
+    Sx2x2 += x2[i] * x2[i];
+    Sx1y += x1[i] * y[i];
+    Sx2y += x2[i] * y[i];
+    Sy += y[i];
+  }
+
+  U = n * (Sx1x2 * Sx1y - Sx1x1 * Sx2y) +
+      Sx1 * Sx1 * Sx2y - Sx1 * Sx2 * Sx1y +
+      Sy * (Sx2 * Sx1x1 - Sx1 * Sx1x2);
+
+  V1 = n * (Sx1x2 * Sx1x2 - Sx1x1 * Sx2x2);
+  V2 = Sx1 * Sx1 * Sx2x2 + Sx2 * Sx2 * Sx1x1;
+  V3 = -2.0 * Sx1 * Sx2 * Sx1x2;
+  V = V1 + V2 + V3;
+
+  /* Check if there is a (numerically stable) solution */
+  if (fabs(V) * 1.0e10 <= -V1 + V2 + fabs(V3))
+    return 0;
+
+  *b2 = U / V;
+
+  return 1;
+}

regress.h

@@ -119,4 +119,16 @@ RGR_FindBestRobustRegression
  int *n_runs,
  int *best_start);
 
+int
+RGR_MultipleRegress
+(double *x1,                    /* first independent variable */
+ double *x2,                    /* second independent variable */
+ double *y,                     /* measured data */
+
+ int n,                         /* number of data points */
+
+ /* The results */
+ double *b2                     /* estimated second slope */
+);
+
 #endif /* GOT_REGRESS_H */