text: #define watson_N         31
  #define watson_P         6
  
  #define watson_NTRIES    4
  
  static double watson_x0[watson_P] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  static double watson_epsrel = 1.0e-6;
  
  static void
  watson_checksol(const double x[], const double sumsq,
                  const double epsrel, const char *sname,
                  const char *pname)
  {
    size_t i;
    const double sumsq_exact = 2.287670053552372e-03;
    const double watson_x[watson_P] = {
      -1.572508640629858e-02,  1.012434869366059e+00, -2.329916259263380e-01,
       1.260430087686035e+00, -1.513728922580576e+00,  9.929964323646112e-01 };
  
    gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",
                 sname, pname);
  
    for (i = 0; i < watson_P; ++i)
      {
        gsl_test_rel(x[i], watson_x[i], epsrel, "%s/%s i=%zu",
                     sname, pname, i);
      }
  }
  
  static int
  watson_f (const gsl_vector * x, void *params, gsl_vector * f)
  {
    const double x1 = gsl_vector_get(x, 0);
    const double x2 = gsl_vector_get(x, 1);
    size_t i, j;
  
    for (i = 0; i < watson_N - 2; ++i)
      {
        double ti = (i + 1) / 29.0;
        double tjm1 = 1.0, tjm2 = 1.0;
        double sum1 = 0.0, sum2 = 0.0;
  
        for (j = 0; j < watson_P; ++j)
          {
            double xj = gsl_vector_get(x, j);
  
            sum1 += xj * tjm1;
            tjm1 *= ti;
  
            if (j > 0)
              {
                sum2 += j * xj * tjm2;
                tjm2 *= ti;
              }
          }
  
        gsl_vector_set (f, i, sum2 - sum1*sum1 - 1.0);
      }
  
    gsl_vector_set(f, watson_N - 2, x1);
    gsl_vector_set(f, watson_N - 1, x2 - x1*x1 - 1.0);
  
    return GSL_SUCCESS;
  }
  
  static int
  watson_df (const gsl_vector * x, void *params, gsl_matrix * J)
  {
    double x1 = gsl_vector_get (x, 0);
    size_t i, j;
  
    gsl_matrix_set_zero(J);
  
    for (i = 0; i < watson_N - 2; ++i)
      {
        double ti = (i + 1) / 29.0;
        double tjm1 = 1.0, tjm2 = 1.0;
        double sum1 = 0.0;
  
        for (j = 0; j < watson_P; ++j)
          {
            double xj = gsl_vector_get(x, j);
            sum1 += xj * tjm1;
            tjm1 *= ti;
          }
  
        tjm1 = 1.0;
        tjm2 = 1.0;
        for (j = 0; j < watson_P; ++j)
          {
            gsl_matrix_set(J, i, j, j * tjm2 - 2.0*sum1*tjm1);
            tjm1 *= ti;
  
            if (j > 0)
              tjm2 *= ti;
          }
      }
  
    gsl_matrix_set(J, watson_N - 2, 0, 1.0);
    gsl_matrix_set(J, watson_N - 1, 0, -2.0*x1);
    gsl_matrix_set(J, watson_N - 1, 1, 1.0);
  
    return GSL_SUCCESS;
  }
  
  static gsl_multifit_function_fdf watson_func =
  {
    &watson_f,
    &watson_df,
    NULL,
    watson_N,
    watson_P,
    NULL,
    0,
    0
  };
  
  static test_fdf_problem watson_problem =
  {
    "watson",
    watson_x0,
    NULL,
    &watson_epsrel,
    watson_NTRIES,
    &watson_checksol,
    &watson_func
  };
  