CR2RE Pipeline Reference Manual 1.7.0
irplib_polynomial.c
1/* $Id: irplib_polynomial.c,v 1.35 2013-01-29 08:43:33 jtaylor Exp $
2 *
3 * This file is part of the ESO Common Pipeline Library
4 * Copyright (C) 2001-2004 European Southern Observatory
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA
19 */
20
21/*
22 * $Author: jtaylor $
23 * $Date: 2013-01-29 08:43:33 $
24 * $Revision: 1.35 $
25 * $Name: not supported by cvs2svn $
26 */
27
28/*-----------------------------------------------------------------------------
29 Includes
30 -----------------------------------------------------------------------------*/
31
32#ifdef HAVE_CONFIG_H
33#include <config.h>
34#endif
35
36#include "irplib_polynomial.h"
37
38/* IRPLIB_SWAP_DOUBLE: */
39#include "irplib_utils.h"
40
41#include <assert.h>
42#include <math.h>
43/* DBL_MAX: */
44#include <float.h>
45
46/*----------------------------------------------------------------------------*/
52/*----------------------------------------------------------------------------*/
54
55/*-----------------------------------------------------------------------------
56 Static functions
57 -----------------------------------------------------------------------------*/
58
59static double irplib_polynomial_eval_2_max(double,
60 double,
61 double,
62 cpl_boolean,
63 double,
64 double);
65
66static double irplib_polynomial_eval_3_max(double,
67 double,
68 double,
69 double,
70 cpl_boolean,
71 double,
72 double,
73 double);
74
75
76static cpl_boolean
77irplib_polynomial_solve_1d_2(double, double, double, double *, double *);
78static cpl_boolean irplib_polynomial_solve_1d_3(double,
79 double,
80 double,
81 double,
82 double *,
83 double *,
84 double *,
85 cpl_boolean *,
86 cpl_boolean *);
87
88static void irplib_polynomial_solve_1d_31(double,
89 double,
90 double *,
91 double *,
92 double *,
93 cpl_boolean *);
94
95static void irplib_polynomial_solve_1d_32(double,
96 double,
97 double,
98 double *,
99 double *,
100 double *,
101 cpl_boolean *);
102
103static void irplib_polynomial_solve_1d_3r(double,
104 double,
105 double,
106 double,
107 double *,
108 double *,
109 double *);
110
111static void irplib_polynomial_solve_1d_3c(double,
112 double,
113 double,
114 double,
115 double,
116 double,
117 double *,
118 double *,
119 double *,
120 cpl_boolean *,
121 cpl_boolean *);
122
123static cpl_error_code irplib_polynomial_solve_1d_4(double,
124 double,
125 double,
126 double,
127 double,
128 cpl_size *,
129 double *,
130 double *,
131 double *,
132 double *);
133
134static cpl_error_code
135irplib_polynomial_solve_1d_zero(cpl_polynomial *,
136 cpl_vector *,
137 cpl_size *) CPL_ATTR_NONNULL;
138
139static cpl_error_code
140irplib_polynomial_solve_1d_nonzero(cpl_polynomial *,
141 cpl_vector *,
142 cpl_size *) CPL_ATTR_NONNULL;
143
144static cpl_error_code
145irplib_polynomial_divide_1d_root(cpl_polynomial *, double, double *);
146
147static cpl_error_code
148irplib_polynomial_solve_1d_guess(const cpl_polynomial *,
149 double *) CPL_ATTR_NONNULL;
150
151#ifdef IRPLIB_POLYNOMIAL_GUESS_ANASOL
152static double irplib_polynomial_depress_1d(cpl_polynomial *);
153#endif
154
155/*-----------------------------------------------------------------------------
156 Function codes
157 -----------------------------------------------------------------------------*/
158
159/*----------------------------------------------------------------------------*/
185/*----------------------------------------------------------------------------*/
186cpl_error_code
187irplib_polynomial_solve_1d_all(const cpl_polynomial *self,
188 cpl_vector *roots,
189 cpl_size *preal)
190{
191 cpl_error_code error;
192 const cpl_size degree = cpl_polynomial_get_degree(self);
193 cpl_polynomial *p;
194
195 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
196 cpl_ensure_code(roots != NULL, CPL_ERROR_NULL_INPUT);
197 cpl_ensure_code(preal != NULL, CPL_ERROR_NULL_INPUT);
198 cpl_ensure_code(cpl_polynomial_get_dimension(self) == 1,
199 CPL_ERROR_INVALID_TYPE);
200 cpl_ensure_code(degree > 0, CPL_ERROR_DATA_NOT_FOUND);
201 cpl_ensure_code(degree == cpl_vector_get_size(roots),
202 CPL_ERROR_INCOMPATIBLE_INPUT);
203
204 *preal = 0;
205
206 p = cpl_polynomial_duplicate(self);
207
208 error = irplib_polynomial_solve_1d_zero(p, roots, preal);
209
210 if (!error && *preal < degree) {
211 /* There are non-zero roots */
212
213 /* Whether roots need sorting (no need w. up to 4 non-zero roots) */
214 const cpl_boolean dosort = *preal > 0 || degree - *preal > 4;
215
216 assert(cpl_polynomial_get_degree(p) + *preal == degree);
217
218 error = irplib_polynomial_solve_1d_nonzero(p, roots, preal);
219
220 if (!error && dosort) {
221 cpl_vector *reals =
222 cpl_vector_wrap(*preal, cpl_vector_get_data(roots));
223 cpl_vector_sort(reals, CPL_SORT_ASCENDING);
224 (void)cpl_vector_unwrap(reals);
225 }
226 }
227
228 cpl_polynomial_delete(p);
229
230 return error ? cpl_error_set_where(cpl_func) : CPL_ERROR_NONE;
231}
232
233/*----------------------------------------------------------------------------*/
246/*----------------------------------------------------------------------------*/
247static cpl_error_code
248irplib_polynomial_solve_1d_zero(cpl_polynomial *self,
249 cpl_vector *roots,
250 cpl_size *preal)
251{
252 cpl_size nzero;
253 const cpl_size degree = cpl_polynomial_get_degree(self);
254
255 /* Count number of zero-value roots */
256 for (nzero = 0; nzero < degree; nzero++) {
257 if (cpl_polynomial_get_coeff(self, &nzero) != 0.0)
258 break;
259 }
260
261 if (nzero > 0) {
262 cpl_size i = 0;
263 for (; i <= degree - nzero; i++) {
264 const cpl_size icopy = i + nzero;
265 const double value = cpl_polynomial_get_coeff(self, &icopy);
266
267 if (cpl_polynomial_set_coeff(self, &i, value))
268 return cpl_error_set_where(cpl_func);
269 }
270 for (; i <= degree; i++) {
271 if (cpl_polynomial_set_coeff(self, &i, 0.0))
272 return cpl_error_set_where(cpl_func);
273 cpl_vector_set(roots, (*preal)++, 0.0);
274 }
275 }
276
277 return CPL_ERROR_NONE;
278}
279
281
282/*----------------------------------------------------------------------------*/
295/*----------------------------------------------------------------------------*/
296static cpl_error_code
297irplib_polynomial_solve_1d_nonzero(cpl_polynomial *self,
298 cpl_vector *roots,
299 cpl_size *preal)
300{
301 cpl_error_code error = CPL_ERROR_NONE;
302 const cpl_size ncoeffs = 1 + cpl_polynomial_get_degree(self);
303
304
305 cpl_ensure_code(ncoeffs > 1, CPL_ERROR_DATA_NOT_FOUND);
306 cpl_ensure_code(*preal >= 0, CPL_ERROR_ILLEGAL_INPUT);
307 cpl_ensure_code(ncoeffs + *preal == 1 + cpl_vector_get_size(roots),
308 CPL_ERROR_INCOMPATIBLE_INPUT);
309
310 switch (ncoeffs) {
311 case 2: {
312 const cpl_size i1 = 1;
313 const cpl_size i0 = 0;
314 const double p1 = cpl_polynomial_get_coeff(self, &i1);
315 const double p0 = cpl_polynomial_get_coeff(self, &i0);
316
317 assert(p1 != 0.0);
318
319 cpl_vector_set(roots, (*preal)++, -p0 / p1);
320 break;
321 }
322 case 3: {
323 const cpl_size i2 = 2;
324 const cpl_size i1 = 1;
325 const cpl_size i0 = 0;
326 const double p2 = cpl_polynomial_get_coeff(self, &i2);
327 const double p1 = cpl_polynomial_get_coeff(self, &i1);
328 const double p0 = cpl_polynomial_get_coeff(self, &i0);
329 double x1, x2;
330
331 assert(p2 != 0.0);
332
333 if (irplib_polynomial_solve_1d_2(p2, p1, p0, &x1, &x2)) {
334 /* This is the complex root in the upper imaginary half-plane */
335 cpl_vector_set(roots, (*preal), x1);
336 cpl_vector_set(roots, (*preal) + 1, x2);
337 }
338 else {
339 cpl_vector_set(roots, (*preal)++, x1);
340 cpl_vector_set(roots, (*preal)++, x2);
341 }
342 break;
343 }
344 case 4: {
345 const cpl_size i3 = 3;
346 const cpl_size i2 = 2;
347 const cpl_size i1 = 1;
348 const cpl_size i0 = 0;
349 const double p3 = cpl_polynomial_get_coeff(self, &i3);
350 const double p2 = cpl_polynomial_get_coeff(self, &i2);
351 const double p1 = cpl_polynomial_get_coeff(self, &i1);
352 const double p0 = cpl_polynomial_get_coeff(self, &i0);
353 double x1, x2, x3;
354
355 assert(p3 != 0.0);
356
357 if (irplib_polynomial_solve_1d_3(p3, p2, p1, p0, &x1, &x2, &x3,
358 NULL, NULL)) {
359 cpl_vector_set(roots, (*preal)++, x1);
360 /* This is the complex root in the upper imaginary half-plane */
361 cpl_vector_set(roots, (*preal), x2);
362 cpl_vector_set(roots, (*preal) + 1, x3);
363 }
364 else {
365 cpl_vector_set(roots, (*preal)++, x1);
366 cpl_vector_set(roots, (*preal)++, x2);
367 cpl_vector_set(roots, (*preal)++, x3);
368 }
369 break;
370 }
371 case 5: {
372 const cpl_size i4 = 4;
373 const cpl_size i3 = 3;
374 const cpl_size i2 = 2;
375 const cpl_size i1 = 1;
376 const cpl_size i0 = 0;
377 const double p4 = cpl_polynomial_get_coeff(self, &i4);
378 const double p3 = cpl_polynomial_get_coeff(self, &i3);
379 const double p2 = cpl_polynomial_get_coeff(self, &i2);
380 const double p1 = cpl_polynomial_get_coeff(self, &i1);
381 const double p0 = cpl_polynomial_get_coeff(self, &i0);
382 double x1, x2, x3, x4;
383 cpl_size nreal;
384
385 assert(p4 != 0.0);
386
387 error = irplib_polynomial_solve_1d_4(p4, p3, p2, p1, p0, &nreal,
388 &x1, &x2, &x3, &x4);
389 if (!error) {
390 cpl_vector_set(roots, (*preal), x1);
391 cpl_vector_set(roots, (*preal) + 1, x2);
392 cpl_vector_set(roots, (*preal) + 2, x3);
393 cpl_vector_set(roots, (*preal) + 3, x4);
394
395 *preal += nreal;
396 }
397 break;
398 }
399
400 default: {
401 /* Try to reduce the problem by finding a single root */
402 double root = 0.0;
403
404 error = irplib_polynomial_solve_1d_guess(self, &root);
405
406 if (!error) {
407 cpl_vector_set(roots, (*preal)++, root);
408
409 irplib_polynomial_divide_1d_root(self, root, NULL);
410
411 error = irplib_polynomial_solve_1d_nonzero(self, roots, preal);
412 }
413
414 break;
415 }
416 }
417
418 return error ? cpl_error_set_where(cpl_func) : CPL_ERROR_NONE;
419}
420
421/*----------------------------------------------------------------------------*/
431/*----------------------------------------------------------------------------*/
432static cpl_error_code
433irplib_polynomial_solve_1d_guess(const cpl_polynomial *self, double *proot)
434{
435 cpl_errorstate prestate = cpl_errorstate_get();
436 cpl_error_code error = CPL_ERROR_NONE;
437 const cpl_size degree = cpl_polynomial_get_degree(self);
438 const cpl_size ncand = 5; /* 2 is enough for current guessing strategy */
439 double rcand[ncand]; /* Stack with first guesses to try */
440 cpl_size icand = 0;
441 size_t ipos = 0, ineg = 0;
442 double rpos[1], rneg[1];
443 cpl_boolean do_bisect = CPL_FALSE;
444 cpl_size itry;
445
446
447 /* If the derivative at the first guess happens to be zero, then
448 the first guess is no good, so try a few different ones. */
449
450 for (itry = 0;; itry++) {
451 switch (itry) {
452 case 0: {
453 /* Try the arithmetic mean of the roots */
454 const double pn0 = cpl_polynomial_get_coeff(self, &degree);
455 const cpl_size n1 = degree - 1;
456 const double pn1 = cpl_polynomial_get_coeff(self, &n1);
457 double rmean;
458
459 assert(pn0 != 0.0);
460
461 rmean = -pn1 / (pn0 * (double)degree);
462
463 rcand[icand++] = rmean;
464
465 break;
466 }
467
468 case 1: {
469 /* Try the geometric mean of the roots */
470 const cpl_size i0 = 0;
471 const double c0 = cpl_polynomial_get_coeff(self, &i0);
472 double rmean;
473
474 assert(c0 != 0.0);
475
476 rmean = pow(fabs(c0), 1.0 / (double)degree);
477
478 rcand[icand++] = rmean;
479
480 break;
481 }
482
483 case 2: {
484 /* Try to get starting guesses with opposite sign residuals */
485 const cpl_size i0 = 0;
486 const double c0 = cpl_polynomial_get_coeff(self, &i0);
487 const double pn0 = cpl_polynomial_get_coeff(self, &degree);
488 const cpl_size n1 = degree - 1;
489 const double pn1 = cpl_polynomial_get_coeff(self, &n1);
490 double rmean;
491
492 assert(pn0 != 0.0);
493
494 rmean = -pn1 / (pn0 * (double)degree);
495
496 rcand[icand++] = rmean + c0;
497 rcand[icand++] = rmean - c0;
498
499 break;
500 }
501
502#ifdef IRPLIB_POLYNOMIAL_GUESS_ANASOL
503 case 3: {
504 /* Try an analytical solution to a (shifted) monomial */
505 cpl_polynomial *copy = cpl_polynomial_duplicate(self);
506 const cpl_size i0 = 0;
507 const double rmean = irplib_polynomial_depress_1d(copy);
508 const double c0 = cpl_polynomial_get_coeff(copy, &i0);
509 const double radius = pow(fabs(c0), 1.0 / (double)degree);
510
511 rcand[icand++] = rmean + radius;
512 if (radius != 0.0) /* Should always be true */
513 rcand[icand++] = rmean - radius;
514
515 cpl_polynomial_delete(copy);
516
517 break;
518 }
519#endif
520
521 default:
522 /* From here on only first guesses increasingly refined via
523 bisection are tried */
524 if (ipos > 0 && ineg > 0) {
525 rcand[icand++] = 0.5 * (rpos[0] + rneg[0]);
526 do_bisect = CPL_TRUE;
527 }
528 break;
529 }
530
531 if (icand > 0) {
532 double grad;
533 double root = rcand[--icand];
534 const double resid = cpl_polynomial_eval_1d(self, root, &grad);
535
536#ifdef IRPLIB_POLYNOMIAL_DEBUG
537 if (itry > 0)
538 cpl_msg_warning(cpl_func, "RETRY(%d)=%g, degree=%d, r=%g, d=%g",
539 (int)itry, root, (int)degree, resid, grad);
540#endif
541
542 error = cpl_polynomial_solve_1d(self, root, proot, 1);
543 if (!error) {
544 cpl_errorstate_set(prestate);
545 break;
546 }
547
548 if (do_bisect) {
549 *(resid > 0.0 ? rpos : rneg) = root;
550 }
551 else {
552 /* Try to collect first guess with residuals with opposite signs */
553 if (resid > 0.0) {
554 if (ipos == 0)
555 rpos[ipos++] = root;
556 }
557 else if (ineg == 0) {
558 rneg[ineg++] = root;
559 }
560
561 if (ipos == 0 || ineg == 0) {
562 const double resid2 =
563 cpl_polynomial_eval_1d(self, *proot, &grad);
564 if (resid2 > 0.0) {
565 if (ipos == 0)
566 rpos[ipos++] = *proot;
567 }
568 else if (ineg == 0) {
569 rneg[ineg++] = *proot;
570 }
571 }
572 }
573 }
574 else {
575 break;
576 }
577 }
578
579 return error ? cpl_error_set_where(cpl_func) : CPL_ERROR_NONE;
580}
581
582/*----------------------------------------------------------------------------*/
596/*----------------------------------------------------------------------------*/
597static cpl_boolean
598irplib_polynomial_solve_1d_2(double p2,
599 double p1,
600 double p0,
601 double *px1,
602 double *px2)
603{
604 const double sqrtD =
605 sqrt(p1 * p1 < 4.0 * p2 * p0 ? 4.0 * p2 * p0 - p1 * p1
606 : p1 * p1 - 4.0 * p2 * p0);
607 cpl_boolean is_complex = CPL_FALSE;
608 double x1 = -0.5 * p1 / p2; /* Double root */
609 double x2;
610
611 /* Compute residual, assuming D == 0 */
612 double res0 = irplib_polynomial_eval_2_max(p2, p1, p0, CPL_FALSE, x1, x1);
613 double res;
614
615 assert(px1 != NULL);
616 assert(px2 != NULL);
617
618 *px2 = *px1 = x1;
619
620 /* Compute residual, assuming D > 0 */
621
622 /* x1 is the root with largest absolute value */
623 if (p1 > 0.0) {
624 x1 = -0.5 * (p1 + sqrtD);
625 }
626 else {
627 x1 = -0.5 * (p1 - sqrtD);
628 }
629 /* Compute smaller root via division to avoid
630 loss of precision due to cancellation */
631 x2 = p0 / x1;
632 x1 /= p2; /* Scale x1 with leading coefficient */
633
634 res = irplib_polynomial_eval_2_max(p2, p1, p0, CPL_FALSE, x1, x2);
635
636 if (res < res0) {
637 res0 = res;
638 if (x2 > x1) {
639 *px1 = x1;
640 *px2 = x2;
641 }
642 else {
643 *px1 = x2;
644 *px2 = x1;
645 }
646 }
647
648 /* Compute residual, assuming D < 0 */
649
650 x1 = -0.5 * p1 / p2; /* Real part of complex root */
651 x2 = 0.5 * sqrtD / fabs(p2); /* Positive, imaginary part of root */
652
653 res = irplib_polynomial_eval_2_max(p2, p1, p0, CPL_TRUE, x1, x2);
654
655 if (res < res0) {
656 *px1 = x1;
657 *px2 = x2;
658 is_complex = CPL_TRUE;
659 }
660
661 return is_complex;
662}
663
664
665/*----------------------------------------------------------------------------*/
678/*----------------------------------------------------------------------------*/
679static double
680irplib_polynomial_eval_2_max(double p2,
681 double p1,
682 double p0,
683 cpl_boolean is_c,
684 double x1,
685 double x2)
686{
687 double res;
688
689 if (is_c) {
690 res = fabs(p0 + x1 * (p1 + x1 * p2) - p2 * x2 * x2);
691 }
692 else {
693 const double r1 = fabs(p0 + x1 * (p1 + x1 * p2));
694 const double r2 = fabs(p0 + x2 * (p1 + x2 * p2));
695
696 res = r1 > r2 ? r1 : r2;
697 }
698
699 return res;
700}
701
702
703/*----------------------------------------------------------------------------*/
718/*----------------------------------------------------------------------------*/
719static double
720irplib_polynomial_eval_3_max(double p3,
721 double p2,
722 double p1,
723 double p0,
724 cpl_boolean is_c,
725 double x1,
726 double x2,
727 double x3)
728{
729 const double r1 = fabs(p0 + x1 * (p1 + x1 * (p2 + x1 * p3)));
730 double res;
731
732 if (is_c) {
733 const double r2 = fabs(p0 + x2 * (p1 + x2 * (p2 + x2 * p3)) -
734 x3 * x3 * (3.0 * p3 * x2 + p2));
735
736 res = r1 > r2 ? r1 : r2;
737 }
738 else {
739 const double r2 = fabs(p0 + x2 * (p1 + x2 * (p2 + x2 * p3)));
740 const double r3 = fabs(p0 + x3 * (p1 + x3 * (p2 + x3 * p3)));
741 res = r1 > r2 ? (r1 > r3 ? r1 : r3) : (r2 > r3 ? r2 : r3);
742 }
743
744 /* cpl_msg_info(cpl_func, "%d: %g (%g)", __LINE__, res, r1); */
745
746 return res;
747}
748
749
750/*----------------------------------------------------------------------------*/
769/*----------------------------------------------------------------------------*/
770static cpl_boolean
771irplib_polynomial_solve_1d_3(double p3,
772 double p2,
773 double p1,
774 double p0,
775 double *px1,
776 double *px2,
777 double *px3,
778 cpl_boolean *pdbl1,
779 cpl_boolean *pdbl2)
780{
781 cpl_boolean is_complex = CPL_FALSE;
782 const double a = p2 / p3;
783 const double b = p1 / p3;
784 const double c = p0 / p3;
785
786 const double q = (a * a - 3.0 * b);
787 const double r = (a * (2.0 * a * a - 9.0 * b) + 27.0 * c);
788
789 const double Q = q / 9.0;
790 const double R = r / 54.0;
791
792 const double Q3 = Q * Q * Q;
793 const double R2 = R * R;
794
795 double x1 = DBL_MAX; /* Fix (false) uninit warning */
796 double x2 = DBL_MAX; /* Fix (false) uninit warning */
797 double x3 = DBL_MAX; /* Fix (false) uninit warning */
798 double xx1 = DBL_MAX; /* Fix (false) uninit warning */
799 double xx2 = DBL_MAX; /* Fix (false) uninit warning */
800 double xx3 = DBL_MAX; /* Fix (false) uninit warning */
801
802 double res = DBL_MAX;
803 cpl_boolean is_first = CPL_TRUE;
804
805 cpl_boolean dbl2;
806
807
808 assert(px1 != NULL);
809
810 if (pdbl1 != NULL)
811 *pdbl1 = CPL_FALSE;
812 if (pdbl2 != NULL)
813 *pdbl2 = CPL_FALSE;
814
815 dbl2 = CPL_FALSE;
816
817 /*
818 All branches (for which the roots are defined) are evaluated, and
819 the branch with the smallest maximum-residual is chosen.
820 When two maximum-residual are identical, preference is given to
821 the purely real solution and if necessary to the solution with a
822 double root.
823 */
824
825 if ((R2 >= Q3 && R != 0.0) || R2 > Q3) {
826 cpl_boolean is_c = CPL_FALSE;
827 double resx;
828
829 irplib_polynomial_solve_1d_3c(a, c, Q, Q3, R, R2, &x1, &x2, &x3, &is_c,
830 &dbl2);
831
832
833 res = resx =
834 irplib_polynomial_eval_3_max(p3, p2, p1, p0, is_c, x1, x2, x3);
835
836 is_first = CPL_FALSE;
837
838 if (pdbl1 != NULL)
839 *pdbl1 = CPL_FALSE;
840 if (!is_c && pdbl2 != NULL)
841 *pdbl2 = dbl2;
842 is_complex = is_c;
843 }
844
845 if (Q > 0.0 && fabs(R / (Q * sqrt(Q))) <= 1.0) {
846 double resx;
847
848 /* this test is actually R2 < Q3, written in a form suitable
849 for exact computation with integers */
850
851 /* assert( Q > 0.0 ); */
852
853 irplib_polynomial_solve_1d_3r(a, c, Q, R, &xx1, &xx2, &xx3);
854
855 resx = irplib_polynomial_eval_3_max(p3, p2, p1, p0, CPL_FALSE, xx1, xx2,
856 xx3);
857
858 if (is_first || (dbl2 ? resx < res : resx <= res)) {
859 is_first = CPL_FALSE;
860 res = resx;
861 x1 = xx1;
862 x2 = xx2;
863 x3 = xx3;
864 if (pdbl1 != NULL)
865 *pdbl1 = CPL_FALSE;
866 if (pdbl2 != NULL)
867 *pdbl2 = CPL_FALSE;
868 is_complex = CPL_FALSE;
869 }
870 }
871
872 if (Q >= 0) {
873 cpl_boolean dbl1 = CPL_FALSE;
874 double resx;
875
876
877 irplib_polynomial_solve_1d_32(a, c, Q, &xx1, &xx2, &xx3, &dbl2);
878
879 resx = irplib_polynomial_eval_3_max(p3, p2, p1, p0, CPL_FALSE, xx1, xx2,
880 xx3);
881 /*
882 cpl_msg_info(cpl_func, "%d: %g = %g - %g (%u)", __LINE__,
883 res - resx, res, resx, is_complex);
884 */
885
886 if (is_first || resx <= res) {
887 /* is_first = CPL_FALSE; */
888 res = resx;
889 x1 = xx1;
890 x2 = xx2;
891 x3 = xx3;
892 if (pdbl1 != NULL)
893 *pdbl1 = CPL_FALSE;
894 if (pdbl2 != NULL)
895 *pdbl2 = dbl2;
896 is_complex = CPL_FALSE;
897 }
898
899
900 /* This branch also covers the case where the depressed cubic
901 polynomial has zero as triple root (i.e. Q == R == 0) */
902
903 irplib_polynomial_solve_1d_31(a, Q, &xx1, &xx2, &xx3, &dbl1);
904
905 resx = irplib_polynomial_eval_3_max(p3, p2, p1, p0, CPL_FALSE, xx1, xx2,
906 xx3);
907
908 if (resx <= res) {
909 /* is_first = CPL_FALSE; */
910 /*res = resx;*/
911 x1 = xx1;
912 x2 = xx2;
913 x3 = xx3;
914 if (pdbl1 != NULL)
915 *pdbl1 = dbl1;
916 if (pdbl2 != NULL)
917 *pdbl2 = CPL_FALSE;
918 is_complex = CPL_FALSE;
919 }
920 }
921
922 if (px2 != NULL && px3 != NULL) {
923 *px1 = x1;
924 *px2 = x2;
925 *px3 = x3;
926 }
927 else if (is_complex) {
928 *px1 = x1;
929 }
930 else {
931 *px1 = x3;
932 }
933
934 return is_complex;
935}
936
937/*----------------------------------------------------------------------------*/
951/*----------------------------------------------------------------------------*/
952static void
953irplib_polynomial_solve_1d_31(double a,
954 double Q,
955 double *px1,
956 double *px2,
957 double *px3,
958 cpl_boolean *pdbl1)
959{
960 const double sqrtQ = sqrt(Q);
961
962 double x1, x2, x3;
963
964 x2 = x1 = -sqrtQ - a / 3.0;
965 x3 = 2.0 * sqrtQ - a / 3.0;
966 if (pdbl1 != NULL)
967 *pdbl1 = CPL_TRUE;
968
969 *px1 = x1;
970 *px2 = x2;
971 *px3 = x3;
972
973 return;
974}
975
976/*----------------------------------------------------------------------------*/
991/*----------------------------------------------------------------------------*/
992static void
993irplib_polynomial_solve_1d_32(double a,
994 double c,
995 double Q,
996 double *px1,
997 double *px2,
998 double *px3,
999 cpl_boolean *pdbl2)
1000{
1001 const double sqrtQ = sqrt(Q);
1002
1003 double x1 = DBL_MAX;
1004 double x2 = DBL_MAX;
1005 double x3 = DBL_MAX;
1006
1007 if (a > 0.0) {
1008 /* a and sqrt(Q) have same sign - or Q is zero */
1009 x1 = -2.0 * sqrtQ - a / 3.0;
1010 /* FIXME: Two small roots with opposite signs may
1011 end up here, with the sign lost for one of them */
1012 x3 = x2 = -a < x1 ? -sqrt(fabs(c / x1)) : sqrt(fabs(c / x1));
1013 if (pdbl2 != NULL)
1014 *pdbl2 = CPL_TRUE;
1015 }
1016 else if (a < 0.0) {
1017 /* a and sqrt(Q) have opposite signs - or Q is zero */
1018 x3 = x2 = sqrtQ - a / 3.0;
1019 x1 = -c / (x2 * x2);
1020 if (pdbl2 != NULL)
1021 *pdbl2 = CPL_TRUE;
1022 }
1023 else {
1024 x1 = -2.0 * sqrtQ;
1025 x3 = x2 = sqrtQ;
1026 if (pdbl2 != NULL)
1027 *pdbl2 = CPL_TRUE;
1028 }
1029
1030 *px1 = x1;
1031 *px2 = x2;
1032 *px3 = x3;
1033
1034 return;
1035}
1036
1037/*----------------------------------------------------------------------------*/
1057/*----------------------------------------------------------------------------*/
1058static void
1059irplib_polynomial_solve_1d_3c(double a,
1060 double c,
1061 double Q,
1062 double Q3,
1063 double R,
1064 double R2,
1065 double *px1,
1066 double *px2,
1067 double *px3,
1068 cpl_boolean *pis_c,
1069 cpl_boolean *pdbl2)
1070{
1071 /* Due to finite precision some double roots may be missed, and
1072 will be considered to be a pair of complex roots z = x +/-
1073 epsilon i close to the real axis. */
1074
1075 /* Another case: A double root, which is small relative to the
1076 last root, may cause this branch to be taken - with the
1077 imaginary part eventually being truncated to zero. */
1078
1079 const double sgnR = (R >= 0 ? 1.0 : -1.0);
1080 const double A = -sgnR * pow(fabs(R) + sqrt(R2 - Q3), 1.0 / 3.0);
1081 const double B = Q / A;
1082
1083 double x1 = DBL_MAX;
1084 double x2 = DBL_MAX;
1085 double x3 = DBL_MAX;
1086 cpl_boolean is_complex = CPL_FALSE;
1087
1088 if ((A > -B && a > 0.0) || (A < -B && a < 0.0)) {
1089 /* A+B has same sign as a */
1090
1091 /* Real part of complex conjugate */
1092 x2 = -0.5 * (A + B) - a / 3.0; /* No cancellation */
1093 /* Positive, imaginary part of complex conjugate */
1094 x3 = 0.5 * CPL_MATH_SQRT3 * fabs(A - B);
1095
1096 x1 = -c / (x2 * x2 + x3 * x3);
1097 }
1098 else {
1099 /* A+B and a have opposite signs - or exactly one is zero */
1100 x1 = A + B - a / 3.0;
1101 /* Positive, imaginary part of complex conjugate */
1102 x3 = 0.5 * CPL_MATH_SQRT3 * fabs(A - B);
1103
1104 if (x3 > 0.0) {
1105 /* Real part of complex conjugate */
1106 x2 = -0.5 * (A + B) - a / 3.0; /* FIXME: Cancellation */
1107 }
1108 else {
1109 x2 = -a < x1 ? -sqrt(fabs(c / x1)) : sqrt(fabs(c / x1));
1110 x3 = 0.0;
1111 }
1112 }
1113
1114 if (x3 > 0.0) {
1115 is_complex = CPL_TRUE;
1116 }
1117 else {
1118 /* Whoaa, the imaginary part was truncated to zero
1119 - return a real, double root */
1120 x3 = x2;
1121 if (pdbl2 != NULL)
1122 *pdbl2 = CPL_TRUE;
1123 }
1124
1125 *px1 = x1;
1126 *px2 = x2;
1127 *px3 = x3;
1128 *pis_c = is_complex;
1129
1130 return;
1131}
1132
1133/*----------------------------------------------------------------------------*/
1148/*----------------------------------------------------------------------------*/
1149static void
1150irplib_polynomial_solve_1d_3r(double a,
1151 double c,
1152 double Q,
1153 double R,
1154 double *px1,
1155 double *px2,
1156 double *px3)
1157{
1158 const double sqrtQ = sqrt(Q);
1159 const double theta = acos(R / (Q * sqrtQ)); /* theta in range [0; pi] */
1160
1161 /* -1.0 <= cos((theta + CPL_MATH_2PI) / 3.0) <= -0.5
1162 -0.5 <= cos((theta - CPL_MATH_2PI) / 3.0) <= 0.5
1163 0.5 <= cos((theta ) / 3.0) <= 1.0 */
1164
1165#define TR1 (-2.0 * sqrtQ * cos(theta / 3.0))
1166#define TR2 (-2.0 * sqrtQ * cos((theta - CPL_MATH_2PI) / 3.0))
1167#define TR3 (-2.0 * sqrtQ * cos((theta + CPL_MATH_2PI) / 3.0))
1168
1169 /* TR1 < TR2 < TR3, except when theta == 0, then TR2 == TR3 */
1170
1171 /* The three roots must be transformed back via subtraction with a/3.
1172 To prevent loss of precision due to cancellation, the root which
1173 is closest to a/3 is computed using the relation
1174 p3 * x1 * x2 * x3 = -p0 */
1175
1176 double x1 = DBL_MAX;
1177 double x2 = DBL_MAX;
1178 double x3 = DBL_MAX;
1179
1180 if (a > 0.0) {
1181 x1 = TR1 - a / 3.0;
1182 if (TR2 > 0.0 && (TR2 + TR3) > 2.0 * a) {
1183 /* FIXME: Cancellation may still effect x3 ? */
1184 x3 = TR3 - a / 3.0;
1185 x2 = -c / (x1 * x3);
1186 }
1187 else {
1188 /* FIXME: Cancellation may still effect x2, especially
1189 if x2, x3 is (almost) a double root, i.e.
1190 if theta is close to zero. */
1191 x2 = TR2 - a / 3.0;
1192
1193 x3 = -c / (x1 * x2);
1194 }
1195 }
1196 else if (a < 0.0) {
1197 x3 = TR3 - a / 3.0;
1198 if (TR2 < 0.0 && (TR1 + TR2) > 2.0 * a) {
1199 x1 = TR1 - a / 3.0;
1200 x2 = -c / (x1 * x3);
1201 }
1202 else {
1203 x2 = TR2 - a / 3.0;
1204 x1 = -c / (x2 * x3);
1205 }
1206 }
1207 else {
1208 x1 = TR1;
1209 x2 = TR2;
1210 x3 = TR3;
1211 }
1212
1213 assert(x1 < x3);
1214
1215 if (x1 > x2) {
1216 /* In absence of round-off:
1217 theta == PI: x1 == x2,
1218 theta < PI: x1 < x2,
1219
1220 The only way x1 could exceed x2 would be due to round-off when
1221 theta is close to PI */
1222
1223 x1 = x2 = 0.5 * (x1 + x2);
1224 }
1225 else if (x2 > x3) {
1226 /* In absence of round-off:
1227 theta == 0: x2 == x3,
1228 theta > 0: x2 < x3,
1229
1230 For small theta:
1231 Round-off can cause x2 to become greater than x3 */
1232
1233 x3 = x2 = 0.5 * (x2 + x3);
1234 }
1235
1236 *px1 = x1;
1237 *px2 = x2;
1238 *px3 = x3;
1239
1240 return;
1241}
1242
1243/*----------------------------------------------------------------------------*/
1261/*----------------------------------------------------------------------------*/
1262static cpl_error_code
1263irplib_polynomial_solve_1d_4(double p4,
1264 double p3,
1265 double p2,
1266 double p1,
1267 double p0,
1268 cpl_size *preal,
1269 double *px1,
1270 double *px2,
1271 double *px3,
1272 double *px4)
1273{
1274 /* Construct the monic, depressed quartic using Horners scheme on 1 / p4 */
1275 const double a = (p2 - 0.375 * p3 * p3 / p4) / p4;
1276 const double b = (p1 - 0.5 * (p2 - 0.25 * p3 * p3 / p4) * p3 / p4) / p4;
1277 const double c =
1278 (p0 - 0.25 * (p1 - 0.25 * (p2 - 0.1875 * p3 * p3 / p4) * p3 / p4) * p3 /
1279 p4) /
1280 p4;
1281
1282 double x1 = DBL_MAX; /* Fix (false) uninit warning */
1283 double x2 = DBL_MAX; /* Fix (false) uninit warning */
1284 double x3 = DBL_MAX; /* Fix (false) uninit warning */
1285 double x4 = DBL_MAX; /* Fix (false) uninit warning */
1286
1287 assert(preal != NULL);
1288 assert(px1 != NULL);
1289 assert(px2 != NULL);
1290 assert(px3 != NULL);
1291 assert(px4 != NULL);
1292
1293 *preal = 4;
1294
1295 if (c == 0.0) {
1296 /* The depressed quartic has zero as root */
1297 /* Since the sum of the roots is zero, at least one is negative
1298 and at least one is positive - unless they are all zero */
1299 cpl_boolean dbl1, dbl2;
1300 const cpl_boolean is_real =
1301 !irplib_polynomial_solve_1d_3(1.0, 0.0, a, b, &x1, &x3, &x4, &dbl1,
1302 &dbl2);
1303
1304 x1 -= 0.25 * p3 / p4;
1305 x2 = -0.25 * p3 / p4;
1306 x3 -= 0.25 * p3 / p4;
1307 if (is_real) {
1308 if (dbl2) {
1309 x4 = x3;
1310 assert(x1 <= x2);
1311 assert(x2 <= x3);
1312 }
1313 else {
1314 x4 -= 0.25 * p3 / p4;
1315 /* Need (only) a guarded swap of x2, x3 */
1316 if (x2 > x3) {
1317 IRPLIB_SWAP_DOUBLE(x2, x3);
1318 }
1319 if (dbl1) {
1320 assert(x1 <= x2); /* The cubic may have 0 as triple root */
1321 assert(x2 <= x3);
1322 assert(x2 <= x4);
1323 }
1324 else {
1325 assert(x1 < x2);
1326 assert(x2 < x4);
1327 }
1328 }
1329 }
1330 else {
1331 *preal = 2;
1332
1333 if (x1 > x2) {
1334 assert(x3 <= x2); /* Don't swap a complex root */
1335
1336 IRPLIB_SWAP_DOUBLE(x1, x2);
1337 }
1338 else {
1339 assert(x3 >= x2);
1340 }
1341 }
1342 }
1343 else if (b == 0.0) {
1344 /* The monic, depressed quartic is a monic, biquadratic equation */
1345 double u1, u2;
1346 const cpl_boolean is_complex =
1347 irplib_polynomial_solve_1d_2(1.0, a, c, &u1, &u2);
1348
1349 if (is_complex) {
1350 /* All four roots are conjugate, complex */
1351 const double norm = sqrt(u1 * u1 + u2 * u2);
1352 const double v1 = sqrt(0.5 * (norm + u1));
1353 const double v2 = u2 / sqrt(2.0 * (norm + u1));
1354
1355
1356 x1 = -0.25 * p3 / p4 - v1;
1357 x3 = -0.25 * p3 / p4 + v1;
1358
1359 x4 = x2 = v2;
1360
1361 *preal = 0;
1362 }
1363 else if (u1 >= 0.0) {
1364 /* All four roots are real */
1365 const double sv1 = sqrt(u1);
1366 const double sv2 = sqrt(u2);
1367
1368
1369 *preal = 4;
1370
1371 x1 = -0.25 * p3 / p4 - sv2;
1372 x2 = -0.25 * p3 / p4 - sv1;
1373 x3 = -0.25 * p3 / p4 + sv1;
1374 x4 = -0.25 * p3 / p4 + sv2;
1375 }
1376 else if (u2 < 0.0) {
1377 /* All four roots are conjugate, complex */
1378 const double sv1 = sqrt(-u2);
1379 const double sv2 = sqrt(-u1);
1380
1381
1382 *preal = 0;
1383
1384 x1 = x3 = -0.25 * p3 / p4;
1385
1386 x2 = sv1;
1387 x4 = sv2;
1388 }
1389 else {
1390 /* Two roots are real, two roots are conjugate, complex */
1391 const double sv1 = sqrt(-u1);
1392 const double sv2 = sqrt(u2);
1393
1394
1395 *preal = 2;
1396
1397 x1 = -0.25 * p3 / p4 - sv2;
1398 x2 = -0.25 * p3 / p4 + sv2;
1399
1400 x3 = -0.25 * p3 / p4;
1401 x4 = sv1;
1402 }
1403 }
1404 else {
1405 /* Need a root from the nested, monic cubic */
1406 const double q2 = -a;
1407 const double q1 = -4.0 * c;
1408 const double q0 = 4.0 * a * c - b * b;
1409 double u1, sqrtd, sqrtrd;
1410 double z1, z2, z3, z4;
1411
1412 cpl_boolean is_complex1, is_complex2;
1413
1414 /* Largest cubic root ensures real square roots when solving the
1415 quartic equation */
1416 (void)irplib_polynomial_solve_1d_3(1.0, q2, q1, q0, &u1, NULL, NULL,
1417 NULL, NULL);
1418
1419
1420 assert(u1 > a);
1421
1422 sqrtd = sqrt(u1 - a);
1423
1424 sqrtrd = 0.5 * b / sqrtd;
1425
1426 is_complex1 = irplib_polynomial_solve_1d_2(1.0, sqrtd,
1427 0.5 * u1 - sqrtrd, &z1, &z2);
1428
1429 is_complex2 = irplib_polynomial_solve_1d_2(1.0, -sqrtd,
1430 0.5 * u1 + sqrtrd, &z3, &z4);
1431
1432 z1 -= 0.25 * p3 / p4;
1433 z3 -= 0.25 * p3 / p4;
1434 if (!is_complex1)
1435 z2 -= 0.25 * p3 / p4;
1436 if (!is_complex2)
1437 z4 -= 0.25 * p3 / p4;
1438
1439 if (!is_complex1 && is_complex2) {
1440 *preal = 2;
1441 x1 = z1;
1442 x2 = z2;
1443 x3 = z3;
1444 x4 = z4;
1445 }
1446 else if (is_complex1 && !is_complex2) {
1447 *preal = 2;
1448 x1 = z3;
1449 x2 = z4;
1450 x3 = z1;
1451 x4 = z2;
1452 }
1453 else if (is_complex1 && is_complex2) {
1454 *preal = 0;
1455
1456 if (z1 < z3 || (z1 == z3 && z2 <= z4)) {
1457 x1 = z1;
1458 x2 = z2;
1459 x3 = z3;
1460 x4 = z4;
1461 }
1462 else {
1463 x1 = z3;
1464 x2 = z4;
1465 x3 = z1;
1466 x4 = z2;
1467 }
1468 }
1469 else {
1470 *preal = 4;
1471
1472 if (z3 >= z2) {
1473 x1 = z1;
1474 x2 = z2;
1475 x3 = z3;
1476 x4 = z4;
1477 }
1478 else if (z4 <= z1) {
1479 x1 = z3;
1480 x2 = z4;
1481 x3 = z1;
1482 x4 = z2;
1483 }
1484 else if (z2 > z4) {
1485 x1 = z3;
1486 x2 = z1;
1487 x3 = z4;
1488 x4 = z2;
1489 }
1490 else {
1491 x1 = z1;
1492 x2 = z3;
1493 x3 = z2;
1494 x4 = z4;
1495 }
1496 }
1497 }
1498
1499 *px1 = x1;
1500 *px2 = x2;
1501 *px3 = x3;
1502 *px4 = x4;
1503
1504 return CPL_ERROR_NONE;
1505}
1506
1507#ifdef IRPLIB_POLYNOMIAL_GUESS_ANASOL
1508/*----------------------------------------------------------------------------*/
1516/*----------------------------------------------------------------------------*/
1517static double
1518irplib_polynomial_depress_1d(cpl_polynomial *self)
1519{
1520 const cpl_size degree = cpl_polynomial_get_degree(self);
1521 const cpl_size nc1 = degree - 1;
1522 const double an = cpl_polynomial_get_coeff(self, &degree);
1523 const double an1 = cpl_polynomial_get_coeff(self, &nc1);
1524 const double rmean = an != 0.0 ? -an1 / (an * (double)degree) : 0.0;
1525
1526
1527 cpl_ensure(degree > 0, CPL_ERROR_DATA_NOT_FOUND, 0.0);
1528
1529 assert(an != 0.0);
1530
1531 if (rmean != 0.0) {
1532 cpl_polynomial_shift_1d(self, 0, rmean);
1533 }
1534
1535 /* Divide polynomial by leading coefficient */
1536 for (cpl_size i = 0; i < nc1; i++) {
1537 const double ai = cpl_polynomial_get_coeff(self, &i) / an;
1538 cpl_polynomial_set_coeff(self, &i, ai);
1539 }
1540
1541 cpl_polynomial_set_coeff(self, &nc1, 0.0); /* Ensure exact values */
1542 cpl_polynomial_set_coeff(self, &degree, 1.0); /* Ensure exact values */
1543
1544 return rmean;
1545}
1546#endif
1547
1548/*----------------------------------------------------------------------------*/
1563/*----------------------------------------------------------------------------*/
1564static cpl_error_code
1565irplib_polynomial_divide_1d_root(cpl_polynomial *p, double r, double *pres)
1566{
1567 const cpl_size n = cpl_polynomial_get_degree(p);
1568 double sum;
1569 cpl_size i;
1570
1571
1572 cpl_ensure_code(p != NULL, CPL_ERROR_NULL_INPUT);
1573 cpl_ensure_code(n > 0, CPL_ERROR_DATA_NOT_FOUND);
1574
1575 sum = cpl_polynomial_get_coeff(p, &n);
1576 cpl_polynomial_set_coeff(p, &n, 0.0);
1577
1578 for (i = n - 1; i >= 0; i--) {
1579 const double coeff = cpl_polynomial_get_coeff(p, &i);
1580
1581 cpl_polynomial_set_coeff(p, &i, sum);
1582
1583 sum = coeff + r * sum;
1584 }
1585
1586 if (pres != NULL)
1587 *pres = sum;
1588
1589 return CPL_ERROR_NONE;
1590}
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2