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 |
|
|
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; |
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 |
} |
} |