diff --git a/build/targets b/build/targets
index c25cfd1..0de391e 100644
--- a/build/targets
+++ b/build/targets
@@ -29,7 +29,6 @@
                                                 $out/fuzz/fuzz_main.o $
                                                 $out/skcms.o
 
-build $out/src/GaussNewton.o:      compile src/GaussNewton.c
 build $out/src/ICCProfile.o:       compile src/ICCProfile.c
 build $out/src/LinearAlgebra.o:    compile src/LinearAlgebra.c
 build $out/src/PortableMath.o:     compile src/PortableMath.c
diff --git a/skcms.c b/skcms.c
index 80ea143..6b1abf5 100644
--- a/skcms.c
+++ b/skcms.c
@@ -8,7 +8,6 @@
 // skcms.c is a unity build target for skcms, #including every other C source file.
 
 #include "src/Curve.c"
-#include "src/GaussNewton.c"
 #include "src/ICCProfile.c"
 #include "src/LinearAlgebra.c"
 #include "src/PortableMath.c"
diff --git a/skcms.gni b/skcms.gni
index 64eeb16..8577de3 100644
--- a/skcms.gni
+++ b/skcms.gni
@@ -6,8 +6,6 @@
 skcms_sources = [
   "src/Curve.c",
   "src/Curve.h",
-  "src/GaussNewton.c",
-  "src/GaussNewton.h",
   "src/ICCProfile.c",
   "src/LinearAlgebra.c",
   "src/LinearAlgebra.h",
diff --git a/src/GaussNewton.c b/src/GaussNewton.c
deleted file mode 100644
index 61aeddb..0000000
--- a/src/GaussNewton.c
+++ /dev/null
@@ -1,94 +0,0 @@
-/*
- * Copyright 2018 Google Inc.
- *
- * Use of this source code is governed by a BSD-style license that can be
- * found in the LICENSE file.
- */
-
-#include "../skcms.h"
-#include "GaussNewton.h"
-#include "LinearAlgebra.h"
-#include "PortableMath.h"
-#include "TransferFunction.h"
-#include <assert.h>
-#include <string.h>
-
-bool skcms_gauss_newton_step(float (*rg)(float x, const void*, const float P[3], float dfdP[3]),
-                             const void* ctx,
-                             float P[3],
-                             float x0, float dx, int N) {
-    // We'll sample x from the range [x0,x1] (both inclusive) N times with even spacing.
-    //
-    // We want to do P' = P + (Jf^T Jf)^-1 Jf^T r(P),
-    //   where r(P) is the residual vector
-    //   and Jf is the Jacobian matrix of f(), ∂r/∂P.
-    //
-    // Let's review the shape of each of these expressions:
-    //   r(P)   is [N x 1], a column vector with one entry per value of x tested
-    //   Jf     is [N x 3], a matrix with an entry for each (x,P) pair
-    //   Jf^T   is [3 x N], the transpose of Jf
-    //
-    //   Jf^T Jf   is [3 x N] * [N x 3] == [3 x 3], a 3x3 matrix,
-    //                                              and so is its inverse (Jf^T Jf)^-1
-    //   Jf^T r(P) is [3 x N] * [N x 1] == [3 x 1], a column vector with the same shape as P
-    //
-    // Our implementation strategy to get to the final ∆P is
-    //   1) evaluate Jf^T Jf,   call that lhs
-    //   2) evaluate Jf^T r(P), call that rhs
-    //   3) invert lhs
-    //   4) multiply inverse lhs by rhs
-    //
-    // This is a friendly implementation strategy because we don't have to have any
-    // buffers that scale with N, and equally nice don't have to perform any matrix
-    // operations that are variable size.
-    //
-    // Other implementation strategies could trade this off, e.g. evaluating the
-    // pseudoinverse of Jf ( (Jf^T Jf)^-1 Jf^T ) directly, then multiplying that by
-    // the residuals.  That would probably require implementing singular value
-    // decomposition, and would create a [3 x N] matrix to be multiplied by the
-    // [N x 1] residual vector, but on the upside I think that'd eliminate the
-    // possibility of this skcms_gauss_newton_step() function ever failing.
-
-    // 0) start off with lhs and rhs safely zeroed.
-    skcms_Matrix3x3 lhs = {{ {0,0,0}, {0,0,0}, {0,0,0} }};
-    skcms_Vector3   rhs = {  {0,0,0} };
-
-    // 1,2) evaluate lhs and evaluate rhs
-    //   We want to evaluate Jf only once, but both lhs and rhs involve Jf^T,
-    //   so we'll have to update lhs and rhs at the same time.
-    for (int i = 0; i < N; i++) {
-        float x = x0 + i*dx;
-
-        float dfdP[3] = {0,0,0};
-        float resid = rg(x,ctx,P, dfdP);
-
-        for (int r = 0; r < 3; r++) {
-            for (int c = 0; c < 3; c++) {
-                lhs.vals[r][c] += dfdP[r] * dfdP[c];
-            }
-            rhs.vals[r] += dfdP[r] * resid;
-        }
-    }
-
-    // If any of the 3 P parameters are unused, this matrix will be singular.
-    // Detect those cases and fix them up to indentity instead, so we can invert.
-    for (int k = 0; k < 3; k++) {
-        if (lhs.vals[0][k]==0 && lhs.vals[1][k]==0 && lhs.vals[2][k]==0 &&
-            lhs.vals[k][0]==0 && lhs.vals[k][1]==0 && lhs.vals[k][2]==0) {
-            lhs.vals[k][k] = 1;
-        }
-    }
-
-    // 3) invert lhs
-    skcms_Matrix3x3 lhs_inv;
-    if (!skcms_Matrix3x3_invert(&lhs, &lhs_inv)) {
-        return false;
-    }
-
-    // 4) multiply inverse lhs by rhs
-    skcms_Vector3 dP = skcms_MV_mul(&lhs_inv, &rhs);
-    P[0] += dP.vals[0];
-    P[1] += dP.vals[1];
-    P[2] += dP.vals[2];
-    return isfinitef_(P[0]) && isfinitef_(P[1]) && isfinitef_(P[2]);
-}
diff --git a/src/GaussNewton.h b/src/GaussNewton.h
deleted file mode 100644
index 5de9927..0000000
--- a/src/GaussNewton.h
+++ /dev/null
@@ -1,25 +0,0 @@
-/*
- * Copyright 2018 Google Inc.
- *
- * Use of this source code is governed by a BSD-style license that can be
- * found in the LICENSE file.
- */
-
-#pragma once
-
-#include <stdbool.h>
-
-// One Gauss-Newton step, tuning up to 3 parameters P to minimize [ r(x,ctx) ]^2.
-//
-//   rg:       residual function r(x,P) to minimize, and gradient at x in dfdP
-//   ctx:      arbitrary context argument passed to rg
-//   P:        in-out, both your initial guess for parameters of r(), and our updated values
-//   x0,dx,N:  N x-values to test with even dx spacing, [x0, x0+dx, x0+2dx, ...]
-//
-// If you have fewer than 3 parameters, set the unused P to zero, don't touch their dfdP.
-//
-// Returns true and updates P on success, or returns false on failure.
-bool skcms_gauss_newton_step(float (*rg)(float x, const void*, const float P[3], float dfdP[3]),
-                             const void* ctx,
-                             float P[3],
-                             float x0, float dx, int N);
diff --git a/src/TransferFunction.c b/src/TransferFunction.c
index 8f8fb36..0e3467c 100644
--- a/src/TransferFunction.c
+++ b/src/TransferFunction.c
@@ -7,7 +7,6 @@
 
 #include "../skcms.h"
 #include "Curve.h"
-#include "GaussNewton.h"
 #include "LinearAlgebra.h"
 #include "Macros.h"
 #include "PortableMath.h"
@@ -151,19 +150,15 @@
 //    ∂r/∂b =  g(ay + b)^(g-1)
 //          -  g(ad + b)^(g-1)
 
-typedef struct {
-    const skcms_Curve*            curve;
-    const skcms_TransferFunction* tf;
-} rg_nonlinear_arg;
-
 // Return the residual of roundtripping skcms_Curve(x) through f_inv(y) with parameters P,
 // and fill out the gradient of the residual into dfdP.
-static float rg_nonlinear(float x, const void* ctx, const float P[3], float dfdP[3]) {
-    const rg_nonlinear_arg* arg = (const rg_nonlinear_arg*)ctx;
+static float rg_nonlinear(float x,
+                          const skcms_Curve* curve,
+                          const skcms_TransferFunction* tf,
+                          const float P[3],
+                          float dfdP[3]) {
+    const float y = skcms_eval_curve(curve, x);
 
-    const float y = skcms_eval_curve(arg->curve, x);
-
-    const skcms_TransferFunction* tf = arg->tf;
     const float g = P[0],  a = P[1],  b = P[2],
                 c = tf->c, d = tf->d, f = tf->f;
 
@@ -230,6 +225,87 @@
     return lin_points;
 }
 
+static bool gauss_newton_step(const skcms_Curve* curve,
+                              const skcms_TransferFunction* tf,
+                              float P[3],
+                              float x0, float dx, int N) {
+    // We'll sample x from the range [x0,x1] (both inclusive) N times with even spacing.
+    //
+    // We want to do P' = P + (Jf^T Jf)^-1 Jf^T r(P),
+    //   where r(P) is the residual vector
+    //   and Jf is the Jacobian matrix of f(), ∂r/∂P.
+    //
+    // Let's review the shape of each of these expressions:
+    //   r(P)   is [N x 1], a column vector with one entry per value of x tested
+    //   Jf     is [N x 3], a matrix with an entry for each (x,P) pair
+    //   Jf^T   is [3 x N], the transpose of Jf
+    //
+    //   Jf^T Jf   is [3 x N] * [N x 3] == [3 x 3], a 3x3 matrix,
+    //                                              and so is its inverse (Jf^T Jf)^-1
+    //   Jf^T r(P) is [3 x N] * [N x 1] == [3 x 1], a column vector with the same shape as P
+    //
+    // Our implementation strategy to get to the final ∆P is
+    //   1) evaluate Jf^T Jf,   call that lhs
+    //   2) evaluate Jf^T r(P), call that rhs
+    //   3) invert lhs
+    //   4) multiply inverse lhs by rhs
+    //
+    // This is a friendly implementation strategy because we don't have to have any
+    // buffers that scale with N, and equally nice don't have to perform any matrix
+    // operations that are variable size.
+    //
+    // Other implementation strategies could trade this off, e.g. evaluating the
+    // pseudoinverse of Jf ( (Jf^T Jf)^-1 Jf^T ) directly, then multiplying that by
+    // the residuals.  That would probably require implementing singular value
+    // decomposition, and would create a [3 x N] matrix to be multiplied by the
+    // [N x 1] residual vector, but on the upside I think that'd eliminate the
+    // possibility of this gauss_newton_step() function ever failing.
+
+    // 0) start off with lhs and rhs safely zeroed.
+    skcms_Matrix3x3 lhs = {{ {0,0,0}, {0,0,0}, {0,0,0} }};
+    skcms_Vector3   rhs = {  {0,0,0} };
+
+    // 1,2) evaluate lhs and evaluate rhs
+    //   We want to evaluate Jf only once, but both lhs and rhs involve Jf^T,
+    //   so we'll have to update lhs and rhs at the same time.
+    for (int i = 0; i < N; i++) {
+        float x = x0 + i*dx;
+
+        float dfdP[3] = {0,0,0};
+        float resid = rg_nonlinear(x,curve,tf,P, dfdP);
+
+        for (int r = 0; r < 3; r++) {
+            for (int c = 0; c < 3; c++) {
+                lhs.vals[r][c] += dfdP[r] * dfdP[c];
+            }
+            rhs.vals[r] += dfdP[r] * resid;
+        }
+    }
+
+    // If any of the 3 P parameters are unused, this matrix will be singular.
+    // Detect those cases and fix them up to indentity instead, so we can invert.
+    for (int k = 0; k < 3; k++) {
+        if (lhs.vals[0][k]==0 && lhs.vals[1][k]==0 && lhs.vals[2][k]==0 &&
+            lhs.vals[k][0]==0 && lhs.vals[k][1]==0 && lhs.vals[k][2]==0) {
+            lhs.vals[k][k] = 1;
+        }
+    }
+
+    // 3) invert lhs
+    skcms_Matrix3x3 lhs_inv;
+    if (!skcms_Matrix3x3_invert(&lhs, &lhs_inv)) {
+        return false;
+    }
+
+    // 4) multiply inverse lhs by rhs
+    skcms_Vector3 dP = skcms_MV_mul(&lhs_inv, &rhs);
+    P[0] += dP.vals[0];
+    P[1] += dP.vals[1];
+    P[2] += dP.vals[2];
+    return isfinitef_(P[0]) && isfinitef_(P[1]) && isfinitef_(P[2]);
+}
+
+
 // Fit the points in [L,N) to the non-linear piece of tf, or return false if we can't.
 static bool fit_nonlinear(const skcms_Curve* curve, int L, int N, skcms_TransferFunction* tf) {
     float P[3] = { tf->g, tf->a, tf->b };
@@ -250,10 +326,9 @@
         assert (P[1] >= 0 &&
                 P[1] * tf->d + P[2] >= 0);
 
-        rg_nonlinear_arg arg = { curve, tf};
-        if (!skcms_gauss_newton_step(rg_nonlinear, &arg,
-                                     P,
-                                     L*dx, dx, N-L)) {
+        if (!gauss_newton_step(curve, tf,
+                               P,
+                               L*dx, dx, N-L)) {
             return false;
         }
     }