#define rat43_N       15

#define rat43_P       4

static double rat43_x0a[rat43_P] = { 100.0, 10.0, 1.0, 1.0 };

static double rat43_x0b[rat43_P] = { 700.0, 5.0, 0.75, 1.3 };

static double rat43_epsrel = 1.0e-6;

static double rat43_sigma[rat43_P] = {

  1.6302297817E+01, 2.0828735829E+00,

  1.9566123451E-01, 6.8761936385E-01

};

static double rat43_F[rat43_N] = { 16.08, 33.83, 65.80, 97.20, 191.55,

                                   326.20, 386.87, 520.53, 590.03,

                                   651.92, 724.93, 699.56, 689.96,

                                   637.56, 717.41 };

static void

rat43_checksol(const double x[], const double sumsq,

                 const double epsrel, const char *sname,

                 const char *pname)

{

  size_t i;

  const double sumsq_exact = 8.7864049080E+03;

  const double rat43_x[rat43_P] = { 6.9964151270E+02,

                                    5.2771253025E+00,

                                    7.5962938329E-01,

                                    1.2792483859E+00 };

  gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",

               sname, pname);

  for (i = 0; i < rat43_P; ++i)

    {

      gsl_test_rel(x[i], rat43_x[i], epsrel, "%s/%s i=%zu",

                   sname, pname, i);

    }

}

static int

rat43_f (const gsl_vector * x, void *params, gsl_vector * f)

{

  double b[rat43_P];

  size_t i;

  for (i = 0; i < rat43_P; i++)

    {

      b[i] = gsl_vector_get(x, i);

    }

  for (i = 0; i < rat43_N; i++)

    {

      double xi = i + 1.0;

      double e = exp(b[1] - b[2]*xi);

      double yi = b[0] / pow(1.0 + e, 1.0 / b[3]);

      gsl_vector_set (f, i, yi - rat43_F[i]);

    }

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

  return GSL_SUCCESS;

}

static int

rat43_df (const gsl_vector * x, void *params, gsl_matrix * df)

{

  double b[rat43_P];

  size_t i;

  for (i = 0; i < rat43_P; i++)

    {

      b[i] = gsl_vector_get(x, i);

    }

  for (i = 0; i < rat43_N; i++)

    {

      double xi = i + 1.0;

      double e = exp(b[1] - b[2]*xi);

      double term1 = 1.0 + e;

      double term2 = pow(term1, -1.0 / b[3]);

      gsl_matrix_set (df, i, 0, term2);

      gsl_matrix_set (df, i, 1, -b[0] / b[3] * e * term2 / term1);

      gsl_matrix_set (df, i, 2, b[0] / b[3] * xi * e * term2 / term1);

      gsl_matrix_set (df, i, 3, b[0] / b[3] / b[3] * log(term1) * term2);

    }

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

  return GSL_SUCCESS;

}

static gsl_multifit_nlinear_fdf rat43_func =

{

  rat43_f,

  rat43_df,

  NULL, /* analytic expression too complex */

  rat43_N,

  rat43_P,

  NULL,

  0,

  0,

  0

};

static test_fdf_problem rat43a_problem =

{

  "nist-rat43a",

  rat43_x0a,

  NULL,

  rat43_sigma,

  &rat43_epsrel,

  &rat43_checksol,

  &rat43_func

};

static test_fdf_problem rat43b_problem =

{

  "nist-rat43b",

  rat43_x0b,

  NULL,

  rat43_sigma,

  &rat43_epsrel,

  &rat43_checksol,

  &rat43_func

};