test.c

来自「math library from gnu」· C语言 代码 · 共 565 行 · 第 1/2 页

/* poly/test.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_poly.h>

int
main (void)
{
  const double eps = 100.0 * GSL_DBL_EPSILON;

  gsl_ieee_env_setup ();

  /* Polynomial evaluation */

  {
    double x, y;
    double c[3] = { 1.0, 0.5, 0.3 };
    x = 0.5;
    y = gsl_poly_eval (c, 3, x);
    gsl_test_rel (y, 1 + 0.5 * x + 0.3 * x * x, eps,
                  "gsl_poly_eval({1, 0.5, 0.3}, 0.5)");
  }

  {
    double x, y;
    double d[11] = { 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1 };
    x = 1.0;
    y = gsl_poly_eval (d, 11, x);
    gsl_test_rel (y, 1.0, eps,
                  "gsl_poly_eval({1,-1, 1, -1, 1, -1, 1, -1, 1, -1, 1}, 1.0)");

  }

  {
    gsl_complex x, y;
    double c[1] = {0.3};
    GSL_SET_REAL (&x, 0.75);
    GSL_SET_IMAG (&x, 1.2);
    y = gsl_poly_complex_eval (c, 1, x);

    gsl_test_rel (GSL_REAL (y), 0.3, eps, "y.real, gsl_poly_complex_eval ({0.3}, 0.75 + 1.2i)");
    gsl_test_rel (GSL_IMAG (y), 0.0, eps, "y.imag, gsl_poly_complex_eval ({0.3}, 0.75 + 1.2i)");
  }

  {
    gsl_complex x, y;
    double c[4] = {2.1, -1.34, 0.76, 0.45};
    GSL_SET_REAL (&x, 0.49);
    GSL_SET_IMAG (&x, 0.95);
    y = gsl_poly_complex_eval (c, 4, x);

    gsl_test_rel (GSL_REAL (y), 0.3959143, eps, "y.real, gsl_poly_complex_eval ({2.1, -1.34, 0.76, 0.45}, 0.49 + 0.95i)");
    gsl_test_rel (GSL_IMAG (y), -0.6433305, eps, "y.imag, gsl_poly_complex_eval ({2.1, -1.34, 0.76, 0.45}, 0.49 + 0.95i)");
  }

  {
    gsl_complex x, y;
    gsl_complex c[1];
    GSL_SET_REAL (&c[0], 0.674);
    GSL_SET_IMAG (&c[0], -1.423);
    GSL_SET_REAL (&x, -1.44);
    GSL_SET_IMAG (&x, 9.55);
    y = gsl_complex_poly_complex_eval (c, 1, x);

    gsl_test_rel (GSL_REAL (y), 0.674, eps, "y.real, gsl_complex_poly_complex_eval ({0.674 - 1.423i}, -1.44 + 9.55i)");
    gsl_test_rel (GSL_IMAG (y), -1.423, eps, "y.imag, gsl_complex_poly_complex_eval ({0.674 - 1.423i}, -1.44 + 9.55i)");
  }

  {
    gsl_complex x, y;
    gsl_complex c[4];
    GSL_SET_REAL (&c[0], -2.31);
    GSL_SET_IMAG (&c[0], 0.44);
    GSL_SET_REAL (&c[1], 4.21);
    GSL_SET_IMAG (&c[1], -3.19);
    GSL_SET_REAL (&c[2], 0.93);
    GSL_SET_IMAG (&c[2], 1.04);
    GSL_SET_REAL (&c[3], -0.42);
    GSL_SET_IMAG (&c[3], 0.68);
    GSL_SET_REAL (&x, 0.49);
    GSL_SET_IMAG (&x, 0.95);
    y = gsl_complex_poly_complex_eval (c, 4, x);

    gsl_test_rel (GSL_REAL (y), 1.82462012, eps, "y.real, gsl_complex_poly_complex_eval ({-2.31 + 0.44i, 4.21 - 3.19i, 0.93 + 1.04i, -0.42 + 0.68i}, 0.49 + 0.95i)");
    gsl_test_rel (GSL_IMAG (y), 2.30389412, eps, "y.imag, gsl_complex_poly_complex_eval ({-2.31 + 0.44i, 4.21 - 3.19i, 0.93 + 1.04i, -0.42 + 0.68i}, 0.49 + 0.95i)");
  }

  /* Quadratic */

  {
    double x0, x1;

    int n = gsl_poly_solve_quadratic (4.0, -20.0, 26.0, &x0, &x1);

    gsl_test (n != 0, "gsl_poly_solve_quadratic, no roots, (2x - 5)^2 = -1");
  }

  {
    double x0, x1;

    int n = gsl_poly_solve_quadratic (4.0, -20.0, 25.0, &x0, &x1);

    gsl_test (n != 2, "gsl_poly_solve_quadratic, one root, (2x - 5)^2 = 0");
    gsl_test_rel (x0, 2.5, 1e-9, "x0, (2x - 5)^2 = 0");
    gsl_test_rel (x1, 2.5, 1e-9, "x1, (2x - 5)^2 = 0");
    gsl_test (x0 != x1, "x0 == x1, (2x - 5)^2 = 0");
  }

  {
    double x0, x1;

    int n = gsl_poly_solve_quadratic (4.0, -20.0, 21.0, &x0, &x1);

    gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, (2x - 5)^2 = 4");
    gsl_test_rel (x0, 1.5, 1e-9, "x0, (2x - 5)^2 = 4");
    gsl_test_rel (x1, 3.5, 1e-9, "x1, (2x - 5)^2 = 4");
  }

  {
    double x0, x1;

    int n = gsl_poly_solve_quadratic (4.0, 7.0, 0.0, &x0, &x1);

    gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, x(4x + 7) = 0");
    gsl_test_rel (x0, -1.75, 1e-9, "x0, x(4x + 7) = 0");
    gsl_test_rel (x1, 0.0, 1e-9, "x1, x(4x + 7) = 0");
  }

  {
    double x0, x1;

    int n = gsl_poly_solve_quadratic (5.0, 0.0, -20.0, &x0, &x1);

    gsl_test (n != 2,
              "gsl_poly_solve_quadratic, two roots b = 0, 5 x^2 = 20");
    gsl_test_rel (x0, -2.0, 1e-9, "x0, 5 x^2 = 20");
    gsl_test_rel (x1, 2.0, 1e-9, "x1, 5 x^2 = 20");
  }


  {
    double x0, x1;

    int n = gsl_poly_solve_quadratic (0.0, 3.0, -21.0, &x0, &x1);

    gsl_test (n != 1,
              "gsl_poly_solve_quadratic, one root (linear) 3 x - 21 = 0");
    gsl_test_rel (x0, 7.0, 1e-9, "x0, 3x - 21 = 0");
  }


  {
    double x0, x1;
    int n = gsl_poly_solve_quadratic (0.0, 0.0, 1.0, &x0, &x1);

    gsl_test (n != 0,
              "gsl_poly_solve_quadratic, no roots 1 = 0");
  }


  /* Cubic */

  {
    double x0, x1, x2;

    int n = gsl_poly_solve_cubic (0.0, 0.0, -27.0, &x0, &x1, &x2);

    gsl_test (n != 1, "gsl_poly_solve_cubic, one root, x^3 = 27");
    gsl_test_rel (x0, 3.0, 1e-9, "x0, x^3 = 27");
  }

  {
    double x0, x1, x2;

    int n = gsl_poly_solve_cubic (-51.0, 867.0, -4913.0, &x0, &x1, &x2);

    gsl_test (n != 3, "gsl_poly_solve_cubic, three roots, (x-17)^3=0");
    gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)^3=0");
    gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)^3=0");
    gsl_test_rel (x2, 17.0, 1e-9, "x2, (x-17)^3=0");
  }

  {
    double x0, x1, x2;

    int n = gsl_poly_solve_cubic (-57.0, 1071.0, -6647.0, &x0, &x1, &x2);

    gsl_test (n != 3,
              "gsl_poly_solve_cubic, three roots, (x-17)(x-17)(x-23)=0");
    gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-17)(x-23)=0");
    gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)(x-17)(x-23)=0");
    gsl_test_rel (x2, 23.0, 1e-9, "x2, (x-17)(x-17)(x-23)=0");
  }

  {
    double x0, x1, x2;

    int n = gsl_poly_solve_cubic (-11.0, -493.0, +6647.0, &x0, &x1, &x2);

    gsl_test (n != 3,
              "gsl_poly_solve_cubic, three roots, (x+23)(x-17)(x-17)=0");
    gsl_test_rel (x0, -23.0, 1e-9, "x0, (x+23)(x-17)(x-17)=0");
    gsl_test_rel (x1, 17.0, 1e-9, "x1, (x+23)(x-17)(x-17)=0");
    gsl_test_rel (x2, 17.0, 1e-9, "x2, (x+23)(x-17)(x-17)=0");
  }

  {
    double x0, x1, x2;

    int n = gsl_poly_solve_cubic (-143.0, 5087.0, -50065.0, &x0, &x1, &x2);

    gsl_test (n != 3,
              "gsl_poly_solve_cubic, three roots, (x-17)(x-31)(x-95)=0");
    gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-31)(x-95)=0");
    gsl_test_rel (x1, 31.0, 1e-9, "x1, (x-17)(x-31)(x-95)=0");
    gsl_test_rel (x2, 95.0, 1e-9, "x2, (x-17)(x-31)(x-95)=0");
  }

  {
    double x0, x1, x2;

    int n = gsl_poly_solve_cubic (-109.0, 803.0, 50065.0, &x0, &x1, &x2);

    gsl_test (n != 3,
              "gsl_poly_solve_cubic, three roots, (x+17)(x-31)(x-95)=0");
    gsl_test_rel (x0, -17.0, 1e-9, "x0, (x+17)(x-31)(x-95)=0");
    gsl_test_rel (x1, 31.0, 1e-9, "x1, (x+17)(x-31)(x-95)=0");
    gsl_test_rel (x2, 95.0, 1e-9, "x2, (x+17)(x-31)(x-95)=0");
  }

  /* Quadratic with complex roots */

  {
    gsl_complex z0, z1;

    int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 26.0, &z0, &z1);

    gsl_test (n != 2,
              "gsl_poly_complex_solve_quadratic, 2 roots (2x - 5)^2 = -1");
    gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = -1");
    gsl_test_rel (GSL_IMAG (z0), -0.5, 1e-9, "z0.imag, (2x - 5)^2 = -1");

    gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = -1");
    gsl_test_rel (GSL_IMAG (z1), 0.5, 1e-9, "z1.imag, (2x - 5)^2 = -1");
  }

  {
    gsl_complex z0, z1;

    int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 25.0, &z0, &z1);

    gsl_test (n != 2,
              "gsl_poly_complex_solve_quadratic, one root, (2x - 5)^2 = 0");
    gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = 0");
    gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag (2x - 5)^2 = 0");
    gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = 0");
    gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag (2x - 5)^2 = 0");
    gsl_test (GSL_REAL (z0) != GSL_REAL (z1),
              "z0.real == z1.real, (2x - 5)^2 = 0");
    gsl_test (GSL_IMAG (z0) != GSL_IMAG (z1),
              "z0.imag == z1.imag, (2x - 5)^2 = 0");
  }