htswanalysis/MACS/lib/gsl/gsl-1.11/poly/test.c

   1 /* poly/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_test.h>
  24 #include <gsl/gsl_ieee_utils.h>
  25 #include <gsl/gsl_poly.h>
  26
  27 int
  28 main (void)
  29 {
  30   const double eps = 100.0 * GSL_DBL_EPSILON;
  31
  32   gsl_ieee_env_setup ();
  33
  34   /* Polynomial evaluation */
  35
  36   {
  37     double x, y;
  38     double c[3] = { 1.0, 0.5, 0.3 };
  39     x = 0.5;
  40     y = gsl_poly_eval (c, 3, x);
  41     gsl_test_rel (y, 1 + 0.5 * x + 0.3 * x * x, eps,
  42                   "gsl_poly_eval({1, 0.5, 0.3}, 0.5)");
  43   }
  44
  45   {
  46     double x, y;
  47     double d[11] = { 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1 };
  48     x = 1.0;
  49     y = gsl_poly_eval (d, 11, x);
  50     gsl_test_rel (y, 1.0, eps,
  51                   "gsl_poly_eval({1,-1, 1, -1, 1, -1, 1, -1, 1, -1, 1}, 1.0)");
  52
  53   }
  54
  55   {
  56     gsl_complex x, y;
  57     double c[1] = {0.3};
  58     GSL_SET_REAL (&x, 0.75);
  59     GSL_SET_IMAG (&x, 1.2);
  60     y = gsl_poly_complex_eval (c, 1, x);
  61
  62     gsl_test_rel (GSL_REAL (y), 0.3, eps, "y.real, gsl_poly_complex_eval ({0.3}, 0.75 + 1.2i)");
  63     gsl_test_rel (GSL_IMAG (y), 0.0, eps, "y.imag, gsl_poly_complex_eval ({0.3}, 0.75 + 1.2i)");
  64   }
  65
  66   {
  67     gsl_complex x, y;
  68     double c[4] = {2.1, -1.34, 0.76, 0.45};
  69     GSL_SET_REAL (&x, 0.49);
  70     GSL_SET_IMAG (&x, 0.95);
  71     y = gsl_poly_complex_eval (c, 4, x);
  72
  73     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)");
  74     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)");
  75   }
  76
  77   {
  78     gsl_complex x, y;
  79     gsl_complex c[1];
  80     GSL_SET_REAL (&c[0], 0.674);
  81     GSL_SET_IMAG (&c[0], -1.423);
  82     GSL_SET_REAL (&x, -1.44);
  83     GSL_SET_IMAG (&x, 9.55);
  84     y = gsl_complex_poly_complex_eval (c, 1, x);
  85
  86     gsl_test_rel (GSL_REAL (y), 0.674, eps, "y.real, gsl_complex_poly_complex_eval ({0.674 - 1.423i}, -1.44 + 9.55i)");
  87     gsl_test_rel (GSL_IMAG (y), -1.423, eps, "y.imag, gsl_complex_poly_complex_eval ({0.674 - 1.423i}, -1.44 + 9.55i)");
  88   }
  89
  90   {
  91     gsl_complex x, y;
  92     gsl_complex c[4];
  93     GSL_SET_REAL (&c[0], -2.31);
  94     GSL_SET_IMAG (&c[0], 0.44);
  95     GSL_SET_REAL (&c[1], 4.21);
  96     GSL_SET_IMAG (&c[1], -3.19);
  97     GSL_SET_REAL (&c[2], 0.93);
  98     GSL_SET_IMAG (&c[2], 1.04);
  99     GSL_SET_REAL (&c[3], -0.42);
 100     GSL_SET_IMAG (&c[3], 0.68);
 101     GSL_SET_REAL (&x, 0.49);
 102     GSL_SET_IMAG (&x, 0.95);
 103     y = gsl_complex_poly_complex_eval (c, 4, x);
 104
 105     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)");
 106     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)");
 107   }
 108
 109   /* Quadratic */
 110
 111   {
 112     double x0, x1;
 113
 114     int n = gsl_poly_solve_quadratic (4.0, -20.0, 26.0, &x0, &x1);
 115
 116     gsl_test (n != 0, "gsl_poly_solve_quadratic, no roots, (2x - 5)^2 = -1");
 117   }
 118
 119   {
 120     double x0, x1;
 121
 122     int n = gsl_poly_solve_quadratic (4.0, -20.0, 25.0, &x0, &x1);
 123
 124     gsl_test (n != 2, "gsl_poly_solve_quadratic, one root, (2x - 5)^2 = 0");
 125     gsl_test_rel (x0, 2.5, 1e-9, "x0, (2x - 5)^2 = 0");
 126     gsl_test_rel (x1, 2.5, 1e-9, "x1, (2x - 5)^2 = 0");
 127     gsl_test (x0 != x1, "x0 == x1, (2x - 5)^2 = 0");
 128   }
 129
 130   {
 131     double x0, x1;
 132
 133     int n = gsl_poly_solve_quadratic (4.0, -20.0, 21.0, &x0, &x1);
 134
 135     gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, (2x - 5)^2 = 4");
 136     gsl_test_rel (x0, 1.5, 1e-9, "x0, (2x - 5)^2 = 4");
 137     gsl_test_rel (x1, 3.5, 1e-9, "x1, (2x - 5)^2 = 4");
 138   }
 139
 140   {
 141     double x0, x1;
 142
 143     int n = gsl_poly_solve_quadratic (4.0, 7.0, 0.0, &x0, &x1);
 144
 145     gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, x(4x + 7) = 0");
 146     gsl_test_rel (x0, -1.75, 1e-9, "x0, x(4x + 7) = 0");
 147     gsl_test_rel (x1, 0.0, 1e-9, "x1, x(4x + 7) = 0");
 148   }
 149
 150   {
 151     double x0, x1;
 152
 153     int n = gsl_poly_solve_quadratic (5.0, 0.0, -20.0, &x0, &x1);
 154
 155     gsl_test (n != 2,
 156               "gsl_poly_solve_quadratic, two roots b = 0, 5 x^2 = 20");
 157     gsl_test_rel (x0, -2.0, 1e-9, "x0, 5 x^2 = 20");
 158     gsl_test_rel (x1, 2.0, 1e-9, "x1, 5 x^2 = 20");
 159   }
 160
 161
 162   {
 163     double x0, x1;
 164
 165     int n = gsl_poly_solve_quadratic (0.0, 3.0, -21.0, &x0, &x1);
 166
 167     gsl_test (n != 1,
 168               "gsl_poly_solve_quadratic, one root (linear) 3 x - 21 = 0");
 169     gsl_test_rel (x0, 7.0, 1e-9, "x0, 3x - 21 = 0");
 170   }
 171
 172
 173   {
 174     double x0, x1;
 175     int n = gsl_poly_solve_quadratic (0.0, 0.0, 1.0, &x0, &x1);
 176
 177     gsl_test (n != 0,
 178               "gsl_poly_solve_quadratic, no roots 1 = 0");
 179   }
 180
 181
 182   /* Cubic */
 183
 184   {
 185     double x0, x1, x2;
 186
 187     int n = gsl_poly_solve_cubic (0.0, 0.0, -27.0, &x0, &x1, &x2);
 188
 189     gsl_test (n != 1, "gsl_poly_solve_cubic, one root, x^3 = 27");
 190     gsl_test_rel (x0, 3.0, 1e-9, "x0, x^3 = 27");
 191   }
 192
 193   {
 194     double x0, x1, x2;
 195
 196     int n = gsl_poly_solve_cubic (-51.0, 867.0, -4913.0, &x0, &x1, &x2);
 197
 198     gsl_test (n != 3, "gsl_poly_solve_cubic, three roots, (x-17)^3=0");
 199     gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)^3=0");
 200     gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)^3=0");
 201     gsl_test_rel (x2, 17.0, 1e-9, "x2, (x-17)^3=0");
 202   }
 203
 204   {
 205     double x0, x1, x2;
 206
 207     int n = gsl_poly_solve_cubic (-57.0, 1071.0, -6647.0, &x0, &x1, &x2);
 208
 209     gsl_test (n != 3,
 210               "gsl_poly_solve_cubic, three roots, (x-17)(x-17)(x-23)=0");
 211     gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-17)(x-23)=0");
 212     gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)(x-17)(x-23)=0");
 213     gsl_test_rel (x2, 23.0, 1e-9, "x2, (x-17)(x-17)(x-23)=0");
 214   }
 215
 216   {
 217     double x0, x1, x2;
 218
 219     int n = gsl_poly_solve_cubic (-11.0, -493.0, +6647.0, &x0, &x1, &x2);
 220
 221     gsl_test (n != 3,
 222               "gsl_poly_solve_cubic, three roots, (x+23)(x-17)(x-17)=0");
 223     gsl_test_rel (x0, -23.0, 1e-9, "x0, (x+23)(x-17)(x-17)=0");
 224     gsl_test_rel (x1, 17.0, 1e-9, "x1, (x+23)(x-17)(x-17)=0");
 225     gsl_test_rel (x2, 17.0, 1e-9, "x2, (x+23)(x-17)(x-17)=0");
 226   }
 227
 228   {
 229     double x0, x1, x2;
 230
 231     int n = gsl_poly_solve_cubic (-143.0, 5087.0, -50065.0, &x0, &x1, &x2);
 232
 233     gsl_test (n != 3,
 234               "gsl_poly_solve_cubic, three roots, (x-17)(x-31)(x-95)=0");
 235     gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-31)(x-95)=0");
 236     gsl_test_rel (x1, 31.0, 1e-9, "x1, (x-17)(x-31)(x-95)=0");
 237     gsl_test_rel (x2, 95.0, 1e-9, "x2, (x-17)(x-31)(x-95)=0");
 238   }
 239
 240   {
 241     double x0, x1, x2;
 242
 243     int n = gsl_poly_solve_cubic (-109.0, 803.0, 50065.0, &x0, &x1, &x2);
 244
 245     gsl_test (n != 3,
 246               "gsl_poly_solve_cubic, three roots, (x+17)(x-31)(x-95)=0");
 247     gsl_test_rel (x0, -17.0, 1e-9, "x0, (x+17)(x-31)(x-95)=0");
 248     gsl_test_rel (x1, 31.0, 1e-9, "x1, (x+17)(x-31)(x-95)=0");
 249     gsl_test_rel (x2, 95.0, 1e-9, "x2, (x+17)(x-31)(x-95)=0");
 250   }
 251
 252   /* Quadratic with complex roots */
 253
 254   {
 255     gsl_complex z0, z1;
 256
 257     int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 26.0, &z0, &z1);
 258
 259     gsl_test (n != 2,
 260               "gsl_poly_complex_solve_quadratic, 2 roots (2x - 5)^2 = -1");
 261     gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = -1");
 262     gsl_test_rel (GSL_IMAG (z0), -0.5, 1e-9, "z0.imag, (2x - 5)^2 = -1");
 263
 264     gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = -1");
 265     gsl_test_rel (GSL_IMAG (z1), 0.5, 1e-9, "z1.imag, (2x - 5)^2 = -1");
 266   }
 267
 268   {
 269     gsl_complex z0, z1;
 270
 271     int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 25.0, &z0, &z1);
 272
 273     gsl_test (n != 2,
 274               "gsl_poly_complex_solve_quadratic, one root, (2x - 5)^2 = 0");
 275     gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = 0");
 276     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag (2x - 5)^2 = 0");
 277     gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = 0");
 278     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag (2x - 5)^2 = 0");
 279     gsl_test (GSL_REAL (z0) != GSL_REAL (z1),
 280               "z0.real == z1.real, (2x - 5)^2 = 0");
 281     gsl_test (GSL_IMAG (z0) != GSL_IMAG (z1),
 282               "z0.imag == z1.imag, (2x - 5)^2 = 0");
 283   }
 284
 285   {
 286     gsl_complex z0, z1;
 287
 288     int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 21.0, &z0, &z1);
 289
 290     gsl_test (n != 2,
 291               "gsl_poly_complex_solve_quadratic, two roots, (2x - 5)^2 = 4");
 292     gsl_test_rel (GSL_REAL (z0), 1.5, 1e-9, "z0.real, (2x - 5)^2 = 4");
 293     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (2x - 5)^2 = 4");
 294     gsl_test_rel (GSL_REAL (z1), 3.5, 1e-9, "z1.real, (2x - 5)^2 = 4");
 295     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (2x - 5)^2 = 4");
 296   }
 297
 298   {
 299     gsl_complex z0, z1;
 300
 301     int n = gsl_poly_complex_solve_quadratic (4.0, 7.0, 0.0, &z0, &z1);
 302
 303     gsl_test (n != 2,
 304               "gsl_poly_complex_solve_quadratic, two roots, x(4x + 7) = 0");
 305     gsl_test_rel (GSL_REAL (z0), -1.75, 1e-9, "z0.real, x(4x + 7) = 0");
 306     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, x(4x + 7) = 0");
 307     gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, x(4x + 7) = 0");
 308     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, x(4x + 7) = 0");
 309   }
 310
 311   {
 312     gsl_complex z0, z1;
 313
 314     int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, -20.0, &z0, &z1);
 315
 316     gsl_test (n != 2,
 317               "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = 20");
 318     gsl_test_rel (GSL_REAL (z0), -2.0, 1e-9, "z0.real, 5 x^2 = 20");
 319     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, 5 x^2 = 20");
 320     gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, 5 x^2 = 20");
 321     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, 5 x^2 = 20");
 322   }
 323
 324   {
 325     gsl_complex z0, z1;
 326
 327     int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, 20.0, &z0, &z1);
 328
 329     gsl_test (n != 2,
 330               "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = -20");
 331     gsl_test_rel (GSL_REAL (z0), 0.0, 1e-9, "z0.real, 5 x^2 = -20");
 332     gsl_test_rel (GSL_IMAG (z0), -2.0, 1e-9, "z0.imag, 5 x^2 = -20");
 333     gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, 5 x^2 = -20");
 334     gsl_test_rel (GSL_IMAG (z1), 2.0, 1e-9, "z1.imag, 5 x^2 = -20");
 335   }
 336
 337
 338   {
 339     gsl_complex z0, z1;
 340
 341     int n = gsl_poly_complex_solve_quadratic (0.0, 3.0, -21.0, &z0, &z1);
 342
 343     gsl_test (n != 1,
 344               "gsl_poly_complex_solve_quadratic, one root (linear) 3 x - 21 = 0");
 345
 346     gsl_test_rel (GSL_REAL (z0), 7.0, 1e-9, "z0.real, 3x - 21 = 0");
 347     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, 3x - 21 = 0");
 348   }
 349
 350
 351   {
 352     gsl_complex z0, z1;
 353
 354     int n = gsl_poly_complex_solve_quadratic (0.0, 0.0, 1.0, &z0, &z1);
 355     gsl_test (n != 0,
 356               "gsl_poly_complex_solve_quadratic, no roots 1 = 0");
 357   }
 358
 359
 360
 361   /* Cubic with complex roots */
 362
 363   {
 364     gsl_complex z0, z1, z2;
 365
 366     int n = gsl_poly_complex_solve_cubic (0.0, 0.0, -27.0, &z0, &z1, &z2);
 367
 368     gsl_test (n != 3, "gsl_poly_complex_solve_cubic, three root, x^3 = 27");
 369     gsl_test_rel (GSL_REAL (z0), -1.5, 1e-9, "z0.real, x^3 = 27");
 370     gsl_test_rel (GSL_IMAG (z0), -1.5 * sqrt (3.0), 1e-9,
 371                   "z0.imag, x^3 = 27");
 372     gsl_test_rel (GSL_REAL (z1), -1.5, 1e-9, "z1.real, x^3 = 27");
 373     gsl_test_rel (GSL_IMAG (z1), 1.5 * sqrt (3.0), 1e-9, "z1.imag, x^3 = 27");
 374     gsl_test_rel (GSL_REAL (z2), 3.0, 1e-9, "z2.real, x^3 = 27");
 375     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, x^3 = 27");
 376   }
 377
 378   {
 379     gsl_complex z0, z1, z2;
 380
 381     int n = gsl_poly_complex_solve_cubic (-1.0, 1.0, 39.0, &z0, &z1, &z2);
 382
 383     gsl_test (n != 3,
 384               "gsl_poly_complex_solve_cubic, three root, (x+3)(x^2-4x+13) = 0");
 385     gsl_test_rel (GSL_REAL (z0), -3.0, 1e-9, "z0.real, (x+3)(x^2+1) = 0");
 386     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+3)(x^2+1) = 0");
 387     gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, (x+3)(x^2+1) = 0");
 388     gsl_test_rel (GSL_IMAG (z1), -3.0, 1e-9, "z1.imag, (x+3)(x^2+1) = 0");
 389     gsl_test_rel (GSL_REAL (z2), 2.0, 1e-9, "z2.real, (x+3)(x^2+1) = 0");
 390     gsl_test_rel (GSL_IMAG (z2), 3.0, 1e-9, "z2.imag, (x+3)(x^2+1) = 0");
 391   }
 392
 393   {
 394     gsl_complex z0, z1, z2;
 395
 396     int n =
 397       gsl_poly_complex_solve_cubic (-51.0, 867.0, -4913.0, &z0, &z1, &z2);
 398
 399     gsl_test (n != 3,
 400               "gsl_poly_complex_solve_cubic, three roots, (x-17)^3=0");
 401     gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)^3=0");
 402     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)^3=0");
 403     gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)^3=0");
 404     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)^3=0");
 405     gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x-17)^3=0");
 406     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)^3=0");
 407   }
 408
 409   {
 410     gsl_complex z0, z1, z2;
 411
 412     int n =
 413       gsl_poly_complex_solve_cubic (-57.0, 1071.0, -6647.0, &z0, &z1, &z2);
 414
 415     gsl_test (n != 3,
 416               "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-17)(x-23)=0");
 417     gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-17)(x-23)=0");
 418     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-17)(x-23)=0");
 419     gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)(x-17)(x-23)=0");
 420     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-17)(x-23)=0");
 421     gsl_test_rel (GSL_REAL (z2), 23.0, 1e-9, "z2.real, (x-17)(x-17)(x-23)=0");
 422     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-17)(x-23)=0");
 423   }
 424
 425   {
 426     gsl_complex z0, z1, z2;
 427
 428     int n =
 429       gsl_poly_complex_solve_cubic (-11.0, -493.0, +6647.0, &z0, &z1, &z2);
 430
 431     gsl_test (n != 3,
 432               "gsl_poly_complex_solve_cubic, three roots, (x+23)(x-17)(x-17)=0");
 433     gsl_test_rel (GSL_REAL (z0), -23.0, 1e-9,
 434                   "z0.real, (x+23)(x-17)(x-17)=0");
 435     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+23)(x-17)(x-17)=0");
 436     gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x+23)(x-17)(x-17)=0");
 437     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x+23)(x-17)(x-17)=0");
 438     gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x+23)(x-17)(x-17)=0");
 439     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x+23)(x-17)(x-17)=0");
 440   }
 441
 442
 443   {
 444     gsl_complex z0, z1, z2;
 445
 446     int n =
 447       gsl_poly_complex_solve_cubic (-143.0, 5087.0, -50065.0, &z0, &z1, &z2);
 448
 449     gsl_test (n != 3,
 450               "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-31)(x-95)=0");
 451     gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-31)(x-95)=0");
 452     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-31)(x-95)=0");
 453     gsl_test_rel (GSL_REAL (z1), 31.0, 1e-9, "z1.real, (x-17)(x-31)(x-95)=0");
 454     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-31)(x-95)=0");
 455     gsl_test_rel (GSL_REAL (z2), 95.0, 1e-9, "z2.real, (x-17)(x-31)(x-95)=0");
 456     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-31)(x-95)=0");
 457   }
 458
 459
 460   {
 461     /* Wilkinson polynomial: y = (x-1)(x-2)(x-3)(x-4)(x-5) */
 462
 463     double a[6] = { -120, 274, -225, 85, -15, 1 };
 464     double z[6*2];
 465
 466     gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (6);
 467
 468     int status = gsl_poly_complex_solve (a, 6, w, z);
 469
 470     gsl_poly_complex_workspace_free (w);
 471
 472     gsl_test (status,
 473               "gsl_poly_complex_solve, 5th-order Wilkinson polynomial");
 474     gsl_test_rel (z[0], 1.0, 1e-9, "z0.real, 5th-order polynomial");
 475     gsl_test_rel (z[1], 0.0, 1e-9, "z0.imag, 5th-order polynomial");
 476     gsl_test_rel (z[2], 2.0, 1e-9, "z1.real, 5th-order polynomial");
 477     gsl_test_rel (z[3], 0.0, 1e-9, "z1.imag, 5th-order polynomial");
 478     gsl_test_rel (z[4], 3.0, 1e-9, "z2.real, 5th-order polynomial");
 479     gsl_test_rel (z[5], 0.0, 1e-9, "z2.imag, 5th-order polynomial");
 480     gsl_test_rel (z[6], 4.0, 1e-9, "z3.real, 5th-order polynomial");
 481     gsl_test_rel (z[7], 0.0, 1e-9, "z3.imag, 5th-order polynomial");
 482     gsl_test_rel (z[8], 5.0, 1e-9, "z4.real, 5th-order polynomial");
 483     gsl_test_rel (z[9], 0.0, 1e-9, "z4.imag, 5th-order polynomial");
 484   }
 485
 486   {
 487     /* : 8-th order polynomial y = x^8 + x^4 + 1 */
 488
 489     double a[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 };
 490     double z[8*2];
 491
 492     double C = 0.5;
 493     double S = sqrt (3.0) / 2.0;
 494
 495     gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (9);
 496
 497     int status = gsl_poly_complex_solve (a, 9, w, z);
 498
 499     gsl_poly_complex_workspace_free (w);
 500
 501     gsl_test (status, "gsl_poly_complex_solve, 8th-order polynomial");
 502
 503     gsl_test_rel (z[0], -S, 1e-9, "z0.real, 8th-order polynomial");
 504     gsl_test_rel (z[1], C, 1e-9, "z0.imag, 8th-order polynomial");
 505     gsl_test_rel (z[2], -S, 1e-9, "z1.real, 8th-order polynomial");
 506     gsl_test_rel (z[3], -C, 1e-9, "z1.imag, 8th-order polynomial");
 507     gsl_test_rel (z[4], -C, 1e-9, "z2.real, 8th-order polynomial");
 508     gsl_test_rel (z[5], S, 1e-9, "z2.imag, 8th-order polynomial");
 509     gsl_test_rel (z[6], -C, 1e-9, "z3.real, 8th-order polynomial");
 510     gsl_test_rel (z[7], -S, 1e-9, "z3.imag, 8th-order polynomial");
 511     gsl_test_rel (z[8], C, 1e-9, "z4.real, 8th-order polynomial");
 512     gsl_test_rel (z[9], S, 1e-9, "z4.imag, 8th-order polynomial");
 513     gsl_test_rel (z[10], C, 1e-9, "z5.real, 8th-order polynomial");
 514     gsl_test_rel (z[11], -S, 1e-9, "z5.imag, 8th-order polynomial");
 515     gsl_test_rel (z[12], S, 1e-9, "z6.real, 8th-order polynomial");
 516     gsl_test_rel (z[13], C, 1e-9, "z6.imag, 8th-order polynomial");
 517     gsl_test_rel (z[14], S, 1e-9, "z7.real, 8th-order polynomial");
 518     gsl_test_rel (z[15], -C, 1e-9, "z7.imag, 8th-order polynomial");
 519
 520   }
 521
 522   {
 523     int i;
 524
 525     double xa[7] = {0.16, 0.97, 1.94, 2.74, 3.58, 3.73, 4.70 };
 526     double ya[7] = {0.73, 1.11, 1.49, 1.84, 2.30, 2.41, 3.07 };
 527
 528     double dd_expected[7] = {  7.30000000000000e-01,
 529                                4.69135802469136e-01,
 530                               -4.34737219941284e-02,
 531                                2.68681098870099e-02,
 532                               -3.22937056934996e-03,
 533                                6.12763259971375e-03,
 534                               -6.45402453527083e-03 };
 535
 536     double dd[7], coeff[7], work[7];
 537
 538     gsl_poly_dd_init (dd, xa, ya, 7);
 539
 540     for (i = 0; i < 7; i++)
 541       {
 542         gsl_test_rel (dd[i], dd_expected[i], 1e-10, "divided difference dd[%d]", i);
 543       }
 544
 545     for (i = 0; i < 7; i++)
 546       {
 547         double y = gsl_poly_dd_eval(dd, xa, 7, xa[i]);
 548         gsl_test_rel (y, ya[i], 1e-10, "divided difference y[%d]", i);
 549       }
 550
 551     gsl_poly_dd_taylor (coeff, 1.5, dd, xa, 7, work);
 552
 553     for (i = 0; i < 7; i++)
 554       {
 555         double y = gsl_poly_eval(coeff, 7, xa[i] - 1.5);
 556         gsl_test_rel (y, ya[i], 1e-10, "taylor expansion about 1.5 y[%d]", i);
 557       }
 558   }
 559
 560
 561   /* now summarize the results */
 562
 563   exit (gsl_test_summary ());
 564 }