#define thurber_N       37
#define thurber_P       7

static double thurber_x0a[thurber_P] = { 1000.0, 1000.0, 400.0, 40.0,
                                         0.7, 0.3, 0.03 };
static double thurber_x0b[thurber_P] = { 1300.0, 1500.0, 500.0, 75.0,
                                         1.0, 0.4, 0.05 };
static double thurber_epsrel = 1.0e-6;

static double thurber_sigma[thurber_P] = {
  4.6647963344E+00, 3.9571156086E+01, 2.8698696102E+01,
  5.5675370270E+00, 3.1333340687E-02, 1.4984928198E-02,
  6.5842344623E-03
};

static double thurber_X[thurber_N] = {
  -3.067, -2.981, -2.921, -2.912, -2.840, -2.797, -2.702,
  -2.699, -2.633, -2.481, -2.363, -2.322, -1.501, -1.460,
  -1.274, -1.212, -1.100, -1.046, -0.915, -0.714, -0.566,
  -0.545, -0.400, -0.309, -0.109, -0.103, 0.010, 0.119,
  0.377, 0.790, 0.963, 1.006, 1.115, 1.572, 1.841,
  2.047, 2.200
};

static double thurber_F[thurber_N] = {
  80.574, 84.248, 87.264, 87.195, 89.076, 89.608, 89.868, 90.101,
  92.405, 95.854, 100.696, 101.060, 401.672, 390.724, 567.534,
  635.316, 733.054, 759.087, 894.206, 990.785, 1090.109, 1080.914,
  1122.643, 1178.351, 1260.531, 1273.514, 1288.339, 1327.543, 1353.863,
  1414.509, 1425.208, 1421.384, 1442.962, 1464.350, 1468.705,
  1447.894, 1457.628
};

static void
thurber_checksol(const double x[], const double sumsq,
                 const double epsrel, const char *sname,
                 const char *pname)
{
  size_t i;
  const double sumsq_exact = 5.6427082397E+03;
  const double thurber_x[thurber_P] = {
    1.2881396800E+03, 1.4910792535E+03, 5.8323836877E+02,
    7.5416644291E+01, 9.6629502864E-01, 3.9797285797E-01,
    4.9727297349E-02 };

  gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",
               sname, pname);

  for (i = 0; i < thurber_P; ++i)
    {
      gsl_test_rel(x[i], thurber_x[i], epsrel, "%s/%s i=%zu",
                   sname, pname, i);
    }
}


static int
thurber_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double b[thurber_P];
  size_t i;

  for (i = 0; i < thurber_P; i++)
    {
      b[i] = gsl_vector_get(x, i);
    }

  for (i = 0; i < thurber_N; i++)
    {
      double xi = thurber_X[i];
      double yi;

      yi = b[0] + b[1]*xi + b[2]*xi*xi + b[3]*xi*xi*xi;
      yi /= 1.0 + b[4]*xi + b[5]*xi*xi + b[6]*xi*xi*xi;

      gsl_vector_set (f, i, yi - thurber_F[i]);
    }

  (void)params; /* avoid unused parameter warning */

  return GSL_SUCCESS;
}

static int
thurber_df (const gsl_vector * x, void *params, gsl_matrix * df)
{
  double b[thurber_P];
  size_t i;

  for (i = 0; i < thurber_P; i++)
    {
      b[i] = gsl_vector_get(x, i);
    }

  for (i = 0; i < thurber_N; i++)
    {
      double xi = thurber_X[i];
      double d, n, d_sq;

      n = b[0] + b[1]*xi + b[2]*xi*xi + b[3]*xi*xi*xi;
      d = 1.0 + b[4]*xi + b[5]*xi*xi + b[6]*xi*xi*xi;
      d_sq = d * d;

      gsl_matrix_set (df, i, 0, 1.0 / d);
      gsl_matrix_set (df, i, 1, xi / d);
      gsl_matrix_set (df, i, 2, (xi * xi) / d);
      gsl_matrix_set (df, i, 3, (xi * xi * xi) / d);
      gsl_matrix_set (df, i, 4, -xi * n / d_sq);
      gsl_matrix_set (df, i, 5, -xi * xi * n / d_sq);
      gsl_matrix_set (df, i, 6, -xi * xi * xi * n / d_sq);
    }

  (void)params; /* avoid unused parameter warning */

  return GSL_SUCCESS;
}

static gsl_multifit_nlinear_fdf thurber_func =
{
  thurber_f,
  thurber_df,
  NULL, /* analytic expression too complex */
  thurber_N,
  thurber_P,
  NULL,
  0,
  0,
  0
};

static test_fdf_problem thurbera_problem =
{
  "nist-thurbera",
  thurber_x0a,
  NULL,
  thurber_sigma,
  &thurber_epsrel,
  &thurber_checksol,
  &thurber_func
};

static test_fdf_problem thurberb_problem =
{
  "nist-thurberb",
  thurber_x0b,
  NULL,
  thurber_sigma,
  &thurber_epsrel,
  &thurber_checksol,
  &thurber_func
};