finds complex roots of a polynomial using Durand-Kerner method: http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method */
| 2483 | /* finds complex roots of a polynomial using Durand-Kerner method: |
| 2484 | http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method */ |
| 2485 | double cv::solvePoly( InputArray _coeffs0, OutputArray _roots0, int maxIters ) |
| 2486 | { |
| 2487 | typedef Complex<double> C; |
| 2488 | |
| 2489 | double maxDiff = 0; |
| 2490 | int iter, i, j; |
| 2491 | Mat coeffs0 = _coeffs0.getMat(); |
| 2492 | int ctype = _coeffs0.type(); |
| 2493 | int cdepth = CV_MAT_DEPTH(ctype); |
| 2494 | |
| 2495 | CV_Assert( CV_MAT_DEPTH(ctype) >= CV_32F && CV_MAT_CN(ctype) <= 2 ); |
| 2496 | CV_Assert( coeffs0.rows == 1 || coeffs0.cols == 1 ); |
| 2497 | |
| 2498 | int n = coeffs0.cols + coeffs0.rows - 2; |
| 2499 | |
| 2500 | _roots0.create(n, 1, CV_MAKETYPE(cdepth, 2), -1, true, DEPTH_MASK_FLT); |
| 2501 | Mat roots0 = _roots0.getMat(); |
| 2502 | |
| 2503 | AutoBuffer<C> buf(n*2+2); |
| 2504 | C *coeffs = buf, *roots = coeffs + n + 1; |
| 2505 | Mat coeffs1(coeffs0.size(), CV_MAKETYPE(CV_64F, coeffs0.channels()), coeffs0.channels() == 2 ? coeffs : roots); |
| 2506 | coeffs0.convertTo(coeffs1, coeffs1.type()); |
| 2507 | if( coeffs0.channels() == 1 ) |
| 2508 | { |
| 2509 | const double* rcoeffs = (const double*)roots; |
| 2510 | for( i = 0; i <= n; i++ ) |
| 2511 | coeffs[i] = C(rcoeffs[i], 0); |
| 2512 | } |
| 2513 | |
| 2514 | C p(1, 0), r(1, 1); |
| 2515 | for( i = 0; i < n; i++ ) |
| 2516 | { |
| 2517 | roots[i] = p; |
| 2518 | p = p * r; |
| 2519 | } |
| 2520 | |
| 2521 | maxIters = maxIters <= 0 ? 1000 : maxIters; |
| 2522 | for( iter = 0; iter < maxIters; iter++ ) |
| 2523 | { |
| 2524 | maxDiff = 0; |
| 2525 | for( i = 0; i < n; i++ ) |
| 2526 | { |
| 2527 | p = roots[i]; |
| 2528 | C num = coeffs[n], denom = coeffs[n]; |
| 2529 | int num_same_root = 1; |
| 2530 | for( j = 0; j < n; j++ ) |
| 2531 | { |
| 2532 | num = num*p + coeffs[n-j-1]; |
| 2533 | if( j != i ) |
| 2534 | { |
| 2535 | if ( (p - roots[j]).re != 0 || (p - roots[j]).im != 0 ) |
| 2536 | denom = denom * (p - roots[j]); |
| 2537 | else |
| 2538 | num_same_root++; |
| 2539 | } |
| 2540 | } |
| 2541 | num /= denom; |
| 2542 | if( num_same_root > 1) |