test.c

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

  {
    gsl_complex z0, z1;

    int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, -20.0, &z0, &z1);

    gsl_test (n != 2,
              "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = 20");
    gsl_test_rel (GSL_REAL (z0), -2.0, 1e-9, "z0.real, 5 x^2 = 20");
    gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, 5 x^2 = 20");
    gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, 5 x^2 = 20");
    gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, 5 x^2 = 20");
  }

  {
    gsl_complex z0, z1;

    int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, 20.0, &z0, &z1);

    gsl_test (n != 2,
              "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = -20");
    gsl_test_rel (GSL_REAL (z0), 0.0, 1e-9, "z0.real, 5 x^2 = -20");
    gsl_test_rel (GSL_IMAG (z0), -2.0, 1e-9, "z0.imag, 5 x^2 = -20");
    gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, 5 x^2 = -20");
    gsl_test_rel (GSL_IMAG (z1), 2.0, 1e-9, "z1.imag, 5 x^2 = -20");
  }

  /* Cubic with complex roots */

  {
    gsl_complex z0, z1, z2;

    int n = gsl_poly_complex_solve_cubic (0.0, 0.0, -27.0, &z0, &z1, &z2);

    gsl_test (n != 3, "gsl_poly_complex_solve_cubic, three root, x^3 = 27");
    gsl_test_rel (GSL_REAL (z0), -1.5, 1e-9, "z0.real, x^3 = 27");
    gsl_test_rel (GSL_IMAG (z0), -1.5 * sqrt (3.0), 1e-9,
                  "z0.imag, x^3 = 27");
    gsl_test_rel (GSL_REAL (z1), -1.5, 1e-9, "z1.real, x^3 = 27");
    gsl_test_rel (GSL_IMAG (z1), 1.5 * sqrt (3.0), 1e-9, "z1.imag, x^3 = 27");
    gsl_test_rel (GSL_REAL (z2), 3.0, 1e-9, "z2.real, x^3 = 27");
    gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, x^3 = 27");
  }

  {
    gsl_complex z0, z1, z2;

    int n = gsl_poly_complex_solve_cubic (-1.0, 1.0, 39.0, &z0, &z1, &z2);

    gsl_test (n != 3,
              "gsl_poly_complex_solve_cubic, three root, (x+3)(x^2-4x+13) = 0");
    gsl_test_rel (GSL_REAL (z0), -3.0, 1e-9, "z0.real, (x+3)(x^2+1) = 0");
    gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+3)(x^2+1) = 0");
    gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, (x+3)(x^2+1) = 0");
    gsl_test_rel (GSL_IMAG (z1), -3.0, 1e-9, "z1.imag, (x+3)(x^2+1) = 0");
    gsl_test_rel (GSL_REAL (z2), 2.0, 1e-9, "z2.real, (x+3)(x^2+1) = 0");
    gsl_test_rel (GSL_IMAG (z2), 3.0, 1e-9, "z2.imag, (x+3)(x^2+1) = 0");
  }

  {
    gsl_complex z0, z1, z2;

    int n =
      gsl_poly_complex_solve_cubic (-51.0, 867.0, -4913.0, &z0, &z1, &z2);

    gsl_test (n != 3,
              "gsl_poly_complex_solve_cubic, three roots, (x-17)^3=0");
    gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)^3=0");
    gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)^3=0");
    gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)^3=0");
    gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)^3=0");
    gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x-17)^3=0");
    gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)^3=0");
  }

  {
    gsl_complex z0, z1, z2;

    int n =
      gsl_poly_complex_solve_cubic (-57.0, 1071.0, -6647.0, &z0, &z1, &z2);

    gsl_test (n != 3,
              "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-17)(x-23)=0");
    gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-17)(x-23)=0");
    gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-17)(x-23)=0");
    gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)(x-17)(x-23)=0");
    gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-17)(x-23)=0");
    gsl_test_rel (GSL_REAL (z2), 23.0, 1e-9, "z2.real, (x-17)(x-17)(x-23)=0");
    gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-17)(x-23)=0");
  }

  {
    gsl_complex z0, z1, z2;

    int n =
      gsl_poly_complex_solve_cubic (-11.0, -493.0, +6647.0, &z0, &z1, &z2);

    gsl_test (n != 3,
              "gsl_poly_complex_solve_cubic, three roots, (x+23)(x-17)(x-17)=0");
    gsl_test_rel (GSL_REAL (z0), -23.0, 1e-9,
                  "z0.real, (x+23)(x-17)(x-17)=0");
    gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+23)(x-17)(x-17)=0");
    gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x+23)(x-17)(x-17)=0");
    gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x+23)(x-17)(x-17)=0");
    gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x+23)(x-17)(x-17)=0");
    gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x+23)(x-17)(x-17)=0");
  }


  {
    gsl_complex z0, z1, z2;

    int n =
      gsl_poly_complex_solve_cubic (-143.0, 5087.0, -50065.0, &z0, &z1, &z2);

    gsl_test (n != 3,
              "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-31)(x-95)=0");
    gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-31)(x-95)=0");
    gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-31)(x-95)=0");
    gsl_test_rel (GSL_REAL (z1), 31.0, 1e-9, "z1.real, (x-17)(x-31)(x-95)=0");
    gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-31)(x-95)=0");
    gsl_test_rel (GSL_REAL (z2), 95.0, 1e-9, "z2.real, (x-17)(x-31)(x-95)=0");
    gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-31)(x-95)=0");
  }


  {
    /* Wilkinson polynomial: y = (x-1)(x-2)(x-3)(x-4)(x-5) */

    double a[6] = { -120, 274, -225, 85, -15, 1 };
    double z[6*2];

    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (6);

    int status = gsl_poly_complex_solve (a, 6, w, z);

    gsl_poly_complex_workspace_free (w);

    gsl_test (status,
              "gsl_poly_complex_solve, 5th-order Wilkinson polynomial");
    gsl_test_rel (z[0], 1.0, 1e-9, "z0.real, 5th-order polynomial");
    gsl_test_rel (z[1], 0.0, 1e-9, "z0.imag, 5th-order polynomial");
    gsl_test_rel (z[2], 2.0, 1e-9, "z1.real, 5th-order polynomial");
    gsl_test_rel (z[3], 0.0, 1e-9, "z1.imag, 5th-order polynomial");
    gsl_test_rel (z[4], 3.0, 1e-9, "z2.real, 5th-order polynomial");
    gsl_test_rel (z[5], 0.0, 1e-9, "z2.imag, 5th-order polynomial");
    gsl_test_rel (z[6], 4.0, 1e-9, "z3.real, 5th-order polynomial");
    gsl_test_rel (z[7], 0.0, 1e-9, "z3.imag, 5th-order polynomial");
    gsl_test_rel (z[8], 5.0, 1e-9, "z4.real, 5th-order polynomial");
    gsl_test_rel (z[9], 0.0, 1e-9, "z4.imag, 5th-order polynomial");
  }

  {
    /* : 8-th order polynomial y = x^8 + x^4 + 1 */

    double a[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 };
    double z[8*2];

    double C = 0.5;
    double S = sqrt (3.0) / 2.0;

    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (9);

    int status = gsl_poly_complex_solve (a, 9, w, z);

    gsl_poly_complex_workspace_free (w);

    gsl_test (status, "gsl_poly_complex_solve, 8th-order polynomial");

    gsl_test_rel (z[0], -S, 1e-9, "z0.real, 8th-order polynomial");
    gsl_test_rel (z[1], C, 1e-9, "z0.imag, 8th-order polynomial");
    gsl_test_rel (z[2], -S, 1e-9, "z1.real, 8th-order polynomial");
    gsl_test_rel (z[3], -C, 1e-9, "z1.imag, 8th-order polynomial");
    gsl_test_rel (z[4], -C, 1e-9, "z2.real, 8th-order polynomial");
    gsl_test_rel (z[5], S, 1e-9, "z2.imag, 8th-order polynomial");
    gsl_test_rel (z[6], -C, 1e-9, "z3.real, 8th-order polynomial");
    gsl_test_rel (z[7], -S, 1e-9, "z3.imag, 8th-order polynomial");
    gsl_test_rel (z[8], C, 1e-9, "z4.real, 8th-order polynomial");
    gsl_test_rel (z[9], S, 1e-9, "z4.imag, 8th-order polynomial");
    gsl_test_rel (z[10], C, 1e-9, "z5.real, 8th-order polynomial");
    gsl_test_rel (z[11], -S, 1e-9, "z5.imag, 8th-order polynomial");
    gsl_test_rel (z[12], S, 1e-9, "z6.real, 8th-order polynomial");
    gsl_test_rel (z[13], C, 1e-9, "z6.imag, 8th-order polynomial");
    gsl_test_rel (z[14], S, 1e-9, "z7.real, 8th-order polynomial");
    gsl_test_rel (z[15], -C, 1e-9, "z7.imag, 8th-order polynomial");

  }

  {
    int i;

    double xa[7] = {0.16, 0.97, 1.94, 2.74, 3.58, 3.73, 4.70 };
    double ya[7] = {0.73, 1.11, 1.49, 1.84, 2.30, 2.41, 3.07 };

    double dd_expected[7] = {  7.30000000000000e-01,
                               4.69135802469136e-01,
                              -4.34737219941284e-02,
                               2.68681098870099e-02,
                              -3.22937056934996e-03,
                               6.12763259971375e-03,
                              -6.45402453527083e-03 };

    double dd[7], coeff[7], work[7];
    
    gsl_poly_dd_init (dd, xa, ya, 7);

    for (i = 0; i < 7; i++)
      {
        gsl_test_rel (dd[i], dd_expected[i], 1e-10, "divided difference dd[%d]", i);
      }

    for (i = 0; i < 7; i++)
      {
        double y = gsl_poly_dd_eval(dd, xa, 7, xa[i]);
        gsl_test_rel (y, ya[i], 1e-10, "divided difference y[%d]", i);
      }

    gsl_poly_dd_taylor (coeff, 1.5, dd, xa, 7, work);
    
    for (i = 0; i < 7; i++)
      {
        double y = gsl_poly_eval(coeff, 7, xa[i] - 1.5);
        gsl_test_rel (y, ya[i], 1e-10, "taylor expansion about 1.5 y[%d]", i);
      }
  }


  /* now summarize the results */

  exit (gsl_test_summary ());
}