/[gsl]/gsl/poly/test.c
ViewVC logotype

Diff of /gsl/poly/test.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.1.3 by bjg, Fri Oct 3 09:27:49 2003 UTC revision 1.1.1.3.2.1 by andybuckley, Thu Oct 9 14:27:30 2003 UTC
# Line 355  main (void) Line 355  main (void)
355      gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (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");
356    }    }
357    
358      /* Quartic with complex roots */
359      {
360        gsl_complex z0, z1, z2, z3;
361    
362        int n =
363          gsl_poly_complex_solve_quartic (0.0, 0.0, 0.0, -1.0, &z0, &z1, &z2,
364                                          &z3);
365    
366        gsl_test (n != 4, "gsl_poly_complex_solve_quartic, four roots, x^4=1");
367        gsl_test_rel (GSL_REAL (z0), -1.0, 1e-9, "z0.real, x^4=1");
368        gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, x^4=1");
369        gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, x^4=1");
370        gsl_test_rel (GSL_IMAG (z1), -1.0, 1e-9, "z1.imag, x^4=1");
371        gsl_test_rel (GSL_REAL (z2), 0.0, 1e-9, "z2.real, x^4=1");
372        gsl_test_rel (GSL_IMAG (z2), 1.0, 1e-9, "z2.imag, x^4=1");
373        gsl_test_rel (GSL_REAL (z3), 1.0, 1e-9, "z3.real, x^4=1");
374        gsl_test_rel (GSL_IMAG (z3), 0.0, 1e-9, "z3.imag, x^4=1");
375      }
376    
377      {
378        gsl_complex z0, z1, z2, z3;
379    
380        int n =
381          gsl_poly_complex_solve_quartic (0.0, 0.0, 0.0, 81.0, &z0, &z1, &z2,
382                                          &z3);
383    
384        gsl_test (n != 4, "gsl_poly_complex_solve_quartic, four roots, x^4=-81");
385        gsl_test_rel (GSL_REAL (z0), -3.0 / sqrt (2), 1e-9, "z0.real, x^4=-81");
386        gsl_test_rel (GSL_IMAG (z0), -3.0 / sqrt (2), 1e-9, "z0.imag, x^4=-81");
387        gsl_test_rel (GSL_REAL (z1), -3.0 / sqrt (2), 1e-9, "z1.real, x^4=-81");
388        gsl_test_rel (GSL_IMAG (z1), 3.0 / sqrt (2), 1e-9, "z1.imag, x^4=-81");
389        gsl_test_rel (GSL_REAL (z2), 3.0 / sqrt (2), 1e-9, "z2.real, x^4=-81");
390        gsl_test_rel (GSL_IMAG (z2), -3.0 / sqrt (2), 1e-9, "z2.imag, x^4=-81");
391        gsl_test_rel (GSL_REAL (z3), 3.0 / sqrt (2), 1e-9, "z3.real, x^4=-81");
392        gsl_test_rel (GSL_IMAG (z3), 3.0 / sqrt (2), 1e-9, "z3.imag, x^4=-81");
393      }
394    
395      {
396        gsl_complex z0, z1, z2, z3;
397    
398        int n =
399          gsl_poly_complex_solve_quartic (0.0, 5.0, 0.0, 4.0, &z0, &z1, &z2, &z3);
400    
401        gsl_test (n != 4,
402                  "gsl_poly_complex_solve_quartic, four roots, (x-i)(x+i)(x-2i)(x+2i)=0");
403        gsl_test_rel (GSL_REAL (z0), 0.0, 1e-9,
404                      "z0.real, (x-i)(x+i)(x-2i)(x+2i)=0");
405        gsl_test_rel (GSL_IMAG (z0), -2.0, 1e-9,
406                      "z0.imag, (x-i)(x+i)(x-2i)(x+2i)=0");
407        gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9,
408                      "z1.real, (x-i)(x+i)(x-2i)(x+2i)=0");
409        gsl_test_rel (GSL_IMAG (z1), -1.0, 1e-9,
410                      "z1.imag, (x-i)(x+i)(x-2i)(x+2i)=0");
411        gsl_test_rel (GSL_REAL (z2), 0.0, 1e-9,
412                      "z2.real, (x-i)(x+i)(x-2i)(x+2i)=0");
413        gsl_test_rel (GSL_IMAG (z2), 1.0, 1e-9,
414                      "z2.imag, (x-i)(x+i)(x-2i)(x+2i)=0");
415        gsl_test_rel (GSL_REAL (z3), 0.0, 1e-9,
416                      "z3.real, (x-i)(x+i)(x-2i)(x+2i)=0");
417        gsl_test_rel (GSL_IMAG (z3), 2.0, 1e-9,
418                      "z3.imag, (x-i)(x+i)(x-2i)(x+2i)=0");
419      }
420    
421      {
422        gsl_complex z0, z1, z2, z3;
423    
424        int n =
425          gsl_poly_complex_solve_quartic (-4.0, 11.0, -14.0, 10.0, &z0, &z1, &z2,
426                                          &z3);
427    
428        gsl_test (n != 4,
429                  "gsl_poly_complex_solve_quartic, four roots, (x-1-i)(x-1+i)(x-1-2i)(x-1+2i)=0");
430        gsl_test_rel (GSL_REAL (z0), 1.0, 1e-9,
431                      "z0.real, (x-1-i)(x-1+i)(x-1-2i)(x-1+2i)=0");
432        gsl_test_rel (GSL_IMAG (z0), -2.0, 1e-9,
433                      "z0.imag, (x-1-i)(x-1+i)(x-1-2i)(x-1+2i)=0");
434        gsl_test_rel (GSL_REAL (z1), 1.0, 1e-9,
435                      "z1.real, (x-1-i)(x-1+i)(x-1-2i)(x-1+2i)=0");
436        gsl_test_rel (GSL_IMAG (z1), -1.0, 1e-9,
437                      "z1.imag, (x-1-i)(x-1+i)(x-1-2i)(x-1+2i)=0");
438        gsl_test_rel (GSL_REAL (z2), 1.0, 1e-9,
439                      "z2.real, (x-1-i)(x-1+i)(x-1-2i)(x-1+2i)=0");
440        gsl_test_rel (GSL_IMAG (z2), 1.0, 1e-9,
441                      "z2.imag, (x-1-i)(x-1+i)(x-1-2i)(x-1+2i)=0");
442        gsl_test_rel (GSL_REAL (z3), 1.0, 1e-9,
443                      "z3.real, (x-1-i)(x-1+i)(x-1-2i)(x-1+2i)=0");
444        gsl_test_rel (GSL_IMAG (z3), 2.0, 1e-9,
445                      "z3.imag, (x-1-i)(x-1+i)(x-1-2i)(x-1+2i)=0");
446      }
447    
448    
449    {    {
450      /* Wilkinson polynomial: y = (x-1)(x-2)(x-3)(x-4)(x-5) */      /* Wilkinson polynomial: y = (x-1)(x-2)(x-3)(x-4)(x-5) */
451    
452      double a[6] = { -120, 274, -225, 85, -15, 1 };      double a[6] = { -120, 274, -225, 85, -15, 1 };
453      double z[6*2];      double z[6 * 2];
454    
455      gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (6);      gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (6);
456    
# Line 386  main (void) Line 476  main (void)
476      /* : 8-th order polynomial y = x^8 + x^4 + 1 */      /* : 8-th order polynomial y = x^8 + x^4 + 1 */
477    
478      double a[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 };      double a[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 };
479      double z[8*2];      double z[8 * 2];
480    
481      double C = 0.5;      double C = 0.5;
482      double S = sqrt (3.0) / 2.0;      double S = sqrt (3.0) / 2.0;
# Line 421  main (void) Line 511  main (void)
511    {    {
512      int i;      int i;
513    
514      double xa[7] = {0.16, 0.97, 1.94, 2.74, 3.58, 3.73, 4.70 };      double xa[7] = { 0.16, 0.97, 1.94, 2.74, 3.58, 3.73, 4.70 };
515      double ya[7] = {0.73, 1.11, 1.49, 1.84, 2.30, 2.41, 3.07 };      double ya[7] = { 0.73, 1.11, 1.49, 1.84, 2.30, 2.41, 3.07 };
516    
517      double dd_expected[7] = {  7.30000000000000e-01,      double dd_expected[7] = { 7.30000000000000e-01,
518                                 4.69135802469136e-01,        4.69135802469136e-01,
519                                -4.34737219941284e-02,        -4.34737219941284e-02,
520                                 2.68681098870099e-02,        2.68681098870099e-02,
521                                -3.22937056934996e-03,        -3.22937056934996e-03,
522                                 6.12763259971375e-03,        6.12763259971375e-03,
523                                -6.45402453527083e-03 };        -6.45402453527083e-03
524        };
525    
526      double dd[7], coeff[7], work[7];      double dd[7], coeff[7], work[7];
527        
528      gsl_poly_dd_init (dd, xa, ya, 7);      gsl_poly_dd_init (dd, xa, ya, 7);
529    
530      for (i = 0; i < 7; i++)      for (i = 0; i < 7; i++)
531        {        {
532          gsl_test_rel (dd[i], dd_expected[i], 1e-10, "divided difference dd[%d]", i);          gsl_test_rel (dd[i], dd_expected[i], 1e-10,
533                          "divided difference dd[%d]", i);
534        }        }
535    
536      for (i = 0; i < 7; i++)      for (i = 0; i < 7; i++)
537        {        {
538          double y = gsl_poly_dd_eval(dd, xa, 7, xa[i]);          double y = gsl_poly_dd_eval (dd, xa, 7, xa[i]);
539          gsl_test_rel (y, ya[i], 1e-10, "divided difference y[%d]", i);          gsl_test_rel (y, ya[i], 1e-10, "divided difference y[%d]", i);
540        }        }
541    
542      gsl_poly_dd_taylor (coeff, 1.5, dd, xa, 7, work);      gsl_poly_dd_taylor (coeff, 1.5, dd, xa, 7, work);
543        
544      for (i = 0; i < 7; i++)      for (i = 0; i < 7; i++)
545        {        {
546          double y = gsl_poly_eval(coeff, 7, xa[i] - 1.5);          double y = gsl_poly_eval (coeff, 7, xa[i] - 1.5);
547          gsl_test_rel (y, ya[i], 1e-10, "taylor expansion about 1.5 y[%d]", i);          gsl_test_rel (y, ya[i], 1e-10, "taylor expansion about 1.5 y[%d]", i);
548        }        }
549    }    }

Legend:
Removed from v.1.1.1.3  
changed lines
  Added in v.1.1.1.3.2.1

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26