1 /*
2 Copyright (C) 2018-2021, Dirk Krause
3 SPDX-License-Identifier: BSD-3-Clause
4 */
5
6 /*
7 WARNING: This file was generated by the dkct program (see
8 http://dktools.sourceforge.net/ for details).
9 Changes you make here will be lost if dkct is run again!
10 You should modify the original source and run dkct on it.
11 Original source: dk4iter.ctr
12 */
13
14 /** @file dk4iter.c The dk4iter module.
15 */
16
17
18
19 #include "dk4conf.h"
20
21 #include <stdio.h>
22
23 #if DK4_HAVE_ASSERT_H
24 #ifndef ASSERT_H_INCLUDED
25 #include <assert.h>
26 #define ASSERT_H_INCLUDED 1
27 #endif
28 #endif
29
30 #if DK4_HAVE_STDLIB_H
31 #ifndef STDLIB_H_INCLUDED
32 #include <stdlib.h>
33 #define STDLIB_H_INCLUDED 1
34 #endif
35 #endif
36
37 #if DK4_HAVE_LIMITS_H
38 #ifndef LIMITS_H_INCLUDED
39 #include <limits.h>
40 #define LIMITS_H_INCLUDED 1
41 #endif
42 #endif
43
44 #if DK4_HAVE_STDINT_H
45 #ifndef STDINT_H_INCLUDED
46 #include <stdint.h>
47 #define STDINT_H_INCLUDED 1
48 #endif
49 #endif
50
51
52 #include <libdk4base/dk4mem.h>
53 #include <libdk4c/dk4iter.h>
54 #include <libdk4c/dk4math.h>
55
56
57
58
59
60
61
62 /** Initial function value condition.
63 */
64 enum {
65 CONDITION_ILLEGAL = 0 ,
66 /** fa < 0, fb >= 0.
67 */
68 CONDITION_FA_LESS_ZERO ,
69
70 /** fa <= 0, fb > 0.
71 */
72 CONDITION_FA_LEQ_ZERO ,
73
74 /** fa > 0, fb <= 0.
75 */
76 CONDITION_FA_GREATER_ZERO ,
77
78 /** fa >= 0, fb < 0.
79 */
80 CONDITION_FA_GEQ_ZERO
81 };
82
83
84
85 /** Note which interval border was changed in previous step.
86 */
87 enum {
88 /** No interval border change yet.
89 */
90 CHANGED_NONE = 0,
91
92 /** Previous step changed interval border a.
93 */
94 CHANGED_A = -1,
95
96 /** Previous step changed interval border b.
97 */
98 CHANGED_B = 1
99 };
100
101
102
103 /** Flags which restrictions apply to x values.
104 */
105 enum {
106 /** Minimum specified.
107 */
108 MINMAX_MINIMUM = 0x0001 ,
109
110 /** Maximum specified.
111 */
112 MINMAX_MAXIMUM = 0x0002 ,
113
114 /** Specified minimum is exclusive.
115 */
116 MINMAX_MIN_EXCL = 0x0004 ,
117
118 /** Specified maximum is exclusive.
119 */
120 MINMAX_MAX_EXCL = 0x0008
121 };
122
123
124
125 void
dk4iter_ctx_init(dk4_iter_ctx_t * ctx)126 dk4iter_ctx_init(dk4_iter_ctx_t *ctx)
127 {
128 #if DK4_USE_ASSERT
129 assert(NULL != ctx);
130 #endif
131 if (NULL != ctx) {
132 DK4_MEMRES(ctx,sizeof(dk4_iter_ctx_t));
133 ctx->xmin = 0.0;
134 ctx->xmax = 0.0;
135 ctx->eps_x = 1.0e-8;
136 ctx->eps_y = 1.0e-8;
137 ctx->maxpass = (unsigned long)(DK4_ITER_PASSES_REGULAR);
138 ctx->exact = 0;
139 ctx->algo = DK4_ITER_ALG_RF_ANDERSON_BJOERCK;
140 ctx->minmax = 0;
141 }
142 }
143
144
145
146 dk4_iter_ctx_t *
dk4iter_ctx_open(void)147 dk4iter_ctx_open(void)
148 {
149 dk4_iter_ctx_t *back = NULL;
150 back = dk4mem_new(dk4_iter_ctx_t,1,NULL);
151 if (NULL != back) {
152 dk4iter_ctx_init(back);
153 }
154 return back;
155 }
156
157
158
159 void
dk4iter_ctx_close(dk4_iter_ctx_t * ctx)160 dk4iter_ctx_close(dk4_iter_ctx_t *ctx)
161 {
162 #if DK4_USE_ASSERT
163 assert(NULL != ctx);
164 #endif
165 if (NULL != ctx) {
166 dk4mem_free(ctx);
167 }
168 }
169
170
171
172 void
dk4iter_ctx_set_eps_x(dk4_iter_ctx_t * ctx,double eps)173 dk4iter_ctx_set_eps_x(dk4_iter_ctx_t *ctx, double eps)
174 {
175 #if DK4_USE_ASSERT
176 assert(NULL != ctx);
177 #endif
178 if (NULL != ctx) {
179 ctx->eps_x = eps;
180 }
181 }
182
183
184
185 void
dk4iter_ctx_set_eps_y(dk4_iter_ctx_t * ctx,double eps)186 dk4iter_ctx_set_eps_y(dk4_iter_ctx_t *ctx, double eps)
187 {
188 #if DK4_USE_ASSERT
189 assert(NULL != ctx);
190 #endif
191 if (NULL != ctx) {
192 ctx->eps_y = eps;
193 }
194 }
195
196
197
198 void
dk4iter_ctx_set_maxpass(dk4_iter_ctx_t * ctx,unsigned long passes)199 dk4iter_ctx_set_maxpass(dk4_iter_ctx_t *ctx, unsigned long passes)
200 {
201 #if DK4_USE_ASSERT
202 assert(NULL != ctx);
203 #endif
204 if (NULL != ctx) {
205 ctx->maxpass = passes;
206 }
207 }
208
209
210
211 void
dk4iter_ctx_set_exact(dk4_iter_ctx_t * ctx,int flag)212 dk4iter_ctx_set_exact(dk4_iter_ctx_t *ctx, int flag)
213 {
214 #if DK4_USE_ASSERT
215 assert(NULL != ctx);
216 #endif
217 if (NULL != ctx) {
218 ctx->exact = flag;
219 }
220 }
221
222
223
224 void
dk4iter_ctx_set_algorithm(dk4_iter_ctx_t * ctx,int algorithm)225 dk4iter_ctx_set_algorithm(dk4_iter_ctx_t *ctx, int algorithm)
226 {
227 #if DK4_USE_ASSERT
228 assert(NULL != ctx);
229 #endif
230 if (NULL != ctx) {
231 ctx->algo = algorithm;
232 }
233 }
234
235
236
237 void
dk4iter_ctx_set_min(dk4_iter_ctx_t * ctx,double xmin)238 dk4iter_ctx_set_min(dk4_iter_ctx_t *ctx, double xmin)
239 {
240 #if DK4_USE_ASSERT
241 assert(NULL != ctx);
242 #endif
243 if (NULL != ctx) {
244 ctx->xmin = xmin;
245 ctx->minmax |= MINMAX_MINIMUM;
246 }
247 }
248
249
250
251 void
dk4iter_ctx_set_max(dk4_iter_ctx_t * ctx,double xmax)252 dk4iter_ctx_set_max(dk4_iter_ctx_t *ctx, double xmax)
253 {
254 #if DK4_USE_ASSERT
255 assert(NULL != ctx);
256 #endif
257 if (NULL != ctx) {
258 ctx->xmax = xmax;
259 ctx->minmax |= MINMAX_MAXIMUM;
260 }
261 }
262
263
264
265 void
dk4iter_ctx_set_exclusive_min(dk4_iter_ctx_t * ctx,double xmin)266 dk4iter_ctx_set_exclusive_min(dk4_iter_ctx_t *ctx, double xmin)
267 {
268 #if DK4_USE_ASSERT
269 assert(NULL != ctx);
270 #endif
271 if (NULL != ctx) {
272 ctx->xmin = xmin;
273 ctx->minmax |= (MINMAX_MINIMUM | MINMAX_MIN_EXCL);
274 }
275 }
276
277
278
279 void
dk4iter_ctx_set_exclusive_max(dk4_iter_ctx_t * ctx,double xmax)280 dk4iter_ctx_set_exclusive_max(dk4_iter_ctx_t *ctx, double xmax)
281 {
282 #if DK4_USE_ASSERT
283 assert(NULL != ctx);
284 #endif
285 if (NULL != ctx) {
286 ctx->xmax = xmax;
287 ctx->minmax |= (MINMAX_MAXIMUM | MINMAX_MAX_EXCL);
288 }
289 }
290
291
292
293 static
294 int
find_condition(double fa,double fb)295 find_condition(double fa, double fb)
296 {
297 int back = CONDITION_ILLEGAL;
298
299 if ((0.0 > fa) && (0.0 <= fb)) {
300 back = CONDITION_FA_LESS_ZERO;
301 }
302 else {
303 if ((0.0 < fa) && (0.0 >= fb)) {
304 back = CONDITION_FA_GREATER_ZERO;
305 }
306 else {
307 if ((0.0 >= fa) && (0.0 < fb)) {
308 back = CONDITION_FA_LEQ_ZERO;
309 }
310 else {
311 if ((0.0 <= fa) && (0.0 > fb)) {
312 back = CONDITION_FA_GEQ_ZERO;
313 }
314 }
315 }
316 }
317 return back;
318 }
319
320
321
322 /** Calculate gamma value for PEGASUS variant.
323 @param fg Function value at border last set.
324 @param fx Function value at recent x position.
325 @return Gamma value.
326 */
327 static
328 double
dk4iter_gamma_pegasus(double fg,double fx)329 dk4iter_gamma_pegasus(double fg, double fx)
330 {
331 double back;
332 back = fg / (fg + fx);
333 if (0 == dk4ma_is_finite(back)) {
334 back = 0.5;
335 }
336 else {
337 if (0.0 >= back) {
338 back = 0.5;
339 }
340 #if 0
341 /* 2018-07-06
342 Can not happen, fg and fx have same sign.
343 */
344 else {
345 if (1.0 < back) {
346 back = 1.0;
347 }
348 }
349 #endif
350 }
351 return back;
352 }
353
354
355
356 /** Calculate gamma value for ANDERSON-BJOERCK variant.
357 @param fg Function value at border last set.
358 @param fx Function value at recent x position.
359 @return Gamma value.
360 */
361 static
362 double
dk4iter_gamma_anderson_bjoerck(double fg,double fx)363 dk4iter_gamma_anderson_bjoerck(double fg, double fx)
364 {
365 double back;
366 back = 1.0 - fx / fg;
367 if (0 == dk4ma_is_finite(back)) {
368 back = 0.5;
369 }
370 else {
371 if (0.0 >= back) {
372 back = 0.5;
373 }
374 #if 0
375 /* 2018-07-06
376 Can not happen as fx and fy have same sign.
377 */
378 else {
379 if (1.0 < back) {
380 back = 1.0;
381 }
382 }
383 #endif
384 }
385 return back;
386 }
387
388
389
390 int
dk4iter_interval(double * d,unsigned long * pp,dk4_iter_fct_t * fct,void const * ps,double a,double b,dk4_iter_ctx_t const * ctx)391 dk4iter_interval(
392 double *d,
393 unsigned long *pp,
394 dk4_iter_fct_t *fct,
395 void const *ps,
396 double a,
397 double b,
398 dk4_iter_ctx_t const *ctx
399 )
400 {
401 dk4_iter_ctx_t mctx; /* Copy of context */
402 double fa = 0.0; /* Function value for border a */
403 double fb = 0.0; /* Function value for border b */
404 double x = 0.0; /* Text x value */
405 double fx = 0.0; /* Function value for x */
406 double afx = 0.0; /* Absolute function value for x */
407 double xo = 0.0; /* Previous step x value */
408 double gamma = 0.5; /* Correction factor */
409 unsigned long passno = 0UL; /* Number of current pass */
410 int res = 0; /* Operation result */
411 int cond = 0; /* Condition of interval borders */
412 int cc = 0; /* 1=continue, 0=finished, -1=abort */
413 int pc = CHANGED_NONE; /* Previous border change */
414 int nc = CHANGED_NONE; /* Next change */
415 int back = DK4_ITER_RESULT_E_ARGS;
416
417 /* Check function call arguments
418 */
419 #if DK4_USE_ASSERT
420 assert(NULL != d);
421 assert(NULL != fct);
422 #endif
423 if ((NULL == d) || (NULL == fct)) {
424 goto finished;
425 }
426 if (a == b) {
427 goto finished;
428 }
429
430 /* Copy or initialize context
431 */
432 if (NULL != ctx) {
433 DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t));
434 }
435 else {
436 dk4iter_ctx_init(&mctx);
437 }
438
439 /* Check algorithm specification
440 */
441 if (DK4_ITER_ALG_BISECTION != mctx.algo) {
442 if (DK4_ITER_ALG_RF_PRIMITIVE != mctx.algo) {
443 if (DK4_ITER_ALG_RF_ILLINOIS != mctx.algo) {
444 if (DK4_ITER_ALG_RF_PEGASUS != mctx.algo) {
445 if (DK4_ITER_ALG_RF_ANDERSON_BJOERCK != mctx.algo) {
446 goto finished;
447 }
448 }
449 }
450 }
451 }
452
453 /* Calculate borders at beginning, check
454 */
455 res = (*fct)(&fa, a, ps);
456 if ((0 == res) || (!(dk4ma_is_finite(fa)))) {
457 back = DK4_ITER_RESULT_E_FCT;
458 goto finished;
459 }
460 res = (*fct)(&fb, b, ps);
461 if ((0 == res) || (!(dk4ma_is_finite(fb)))) {
462 back = DK4_ITER_RESULT_E_FCT;
463 goto finished;
464 }
465 cond = find_condition(fa, fb);
466 if (CONDITION_ILLEGAL == cond) {
467 goto finished;
468 }
469
470 /* Prepare iteration loop
471 */
472 passno = 0UL;
473 x = xo = 0.0;
474 cc = 1;
475 pc = CHANGED_NONE;
476
477 /* Run iteration loop
478 */
479 do {
480 /* Keep previous x position
481 */
482 xo = x;
483 /* Increase pass number, but avoid wrapping
484 */
485 if (ULONG_MAX > passno) { passno++; }
486 /* Calculate x position
487 */
488
489 switch (mctx.algo) {
490 case DK4_ITER_ALG_RF_PRIMITIVE :
491 case DK4_ITER_ALG_RF_ILLINOIS :
492 case DK4_ITER_ALG_RF_PEGASUS :
493 case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : { /* Regula falsi */
494 x = (a * fb - b * fa) / (fb - fa);
495 } break;
496 default : { /* Bisection */
497 x = (a + b) / 2.0;
498 } break;
499 }
500 if (0 == dk4ma_is_finite(x)) {
501 back = DK4_ITER_RESULT_E_INFINITE;
502 cc = -1;
503 }
504 else {
505 /* Calculate function for x, if x is usable
506 */
507 res = (*fct)(&fx, x, ps);
508 if (0 == res) {
509 back = DK4_ITER_RESULT_E_FCT;
510 cc = -1;
511 }
512 if (0 == dk4ma_is_finite(fx)) {
513 back = DK4_ITER_RESULT_E_INFINITE;
514 cc = -1;
515 }
516 afx = fabs(fx);
517 if (0 == dk4ma_is_finite(afx)) {
518 back = DK4_ITER_RESULT_E_INFINITE;
519 cc = -1;
520 }
521 /* Check whether we are done, if we could continue
522 */
523 if (1 == cc) {
524
525
526
527 if (isgreater(mctx.eps_y,0.0)) {
528 /* Must check y */
529 if (afx < mctx.eps_y) {
530 if ((x == xo) && (1UL < passno)) {
531 cc = 0;
532 }
533 else {
534 if (0 == mctx.exact) {
535 if(isgreater(mctx.eps_x,0.0)) {
536 if (isless(fabs(x - xo),mctx.eps_x)) {
537 if (1UL < passno) {
538 cc = 0;
539 }
540 }
541 }
542 else {
543 cc = 0;
544 }
545 }
546 }
547 }
548 }
549 else {
550 /* No y check */
551 if ((x == xo) && (1UL < passno)) {
552 cc = 0;
553 }
554 else {
555 if (0 == mctx.exact) {
556 if(isgreater(mctx.eps_x,0.0)) {
557 if (isless(fabs(x - xo),mctx.eps_x)) {
558 if (1UL < passno) {
559 cc = 0;
560 }
561 }
562 }
563 else {
564 /* Neither eps_x nor eps_y,
565 so we wait for x=xo.
566 */
567 }
568 }
569 }
570 }
571 }
572 /* Continue for usable y values only
573 */
574 if (1 == cc) {
575 /* Find direction for next change
576 */
577 nc = CHANGED_NONE;
578 switch (cond) {
579 case CONDITION_FA_LESS_ZERO : {
580 if (0.0 > fx) {
581 nc = CHANGED_A;
582 }
583 else {
584 nc = CHANGED_B;
585 }
586 } break;
587 case CONDITION_FA_LEQ_ZERO : {
588 if (0.0 >= fx) {
589 nc = CHANGED_A;
590 }
591 else {
592 nc = CHANGED_B;
593 }
594 } break;
595 case CONDITION_FA_GREATER_ZERO : {
596 if (0.0 < fx) {
597 nc = CHANGED_A;
598 }
599 else {
600 nc = CHANGED_B;
601 }
602 } break;
603 case CONDITION_FA_GEQ_ZERO : {
604 if (0.0 <= fx) {
605 nc = CHANGED_A;
606 }
607 else {
608 nc = CHANGED_B;
609 }
610 } break;
611 }
612 /* Apply interval border change
613 */
614 switch (nc) {
615 case CHANGED_A: {
616 if (CHANGED_A == pc) {
617 switch (mctx.algo) {
618 case DK4_ITER_ALG_RF_ILLINOIS : {
619 fb = 0.5 * fb;
620 if (0 == dk4ma_is_finite(fb)) {
621 back = DK4_ITER_RESULT_E_INFINITE;
622 cc = -1;
623 }
624 } break;
625 case DK4_ITER_ALG_RF_PEGASUS : {
626 gamma = dk4iter_gamma_pegasus(fa, fx);
627 fb *= gamma;
628 if (0 == dk4ma_is_finite(fb)) {
629 back = DK4_ITER_RESULT_E_INFINITE;
630 cc = -1;
631 }
632 } break;
633 case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : {
634 gamma = dk4iter_gamma_anderson_bjoerck(
635 fa, fx
636 );
637 fb *= gamma;
638 if (0 == dk4ma_is_finite(fb)) {
639 back = DK4_ITER_RESULT_E_INFINITE;
640 cc = -1;
641 }
642 } break;
643 }
644 }
645 a = x;
646 fa = fx;
647 pc = CHANGED_A;
648 } break;
649 case CHANGED_B: {
650 if (CHANGED_B == pc) {
651 switch (mctx.algo) {
652 case DK4_ITER_ALG_RF_ILLINOIS : {
653 fa = 0.5 * fa;
654 if (0 == dk4ma_is_finite(fa)) {
655 back = DK4_ITER_RESULT_E_INFINITE;
656 cc = -1;
657 }
658 } break;
659 case DK4_ITER_ALG_RF_PEGASUS : {
660 gamma = dk4iter_gamma_pegasus(
661 fb, fx
662 );
663 fa *= gamma;
664 if (0 == dk4ma_is_finite(fa)) {
665 back = DK4_ITER_RESULT_E_INFINITE;
666 cc = -1;
667 }
668 } break;
669 case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : {
670 gamma = dk4iter_gamma_anderson_bjoerck(
671 fb, fx
672 );
673 fa *= gamma;
674 if (0 == dk4ma_is_finite(fa)) {
675 back = DK4_ITER_RESULT_E_INFINITE;
676 cc = -1;
677 }
678 } break;
679 }
680 }
681 b = x;
682 fb = fx;
683 pc = CHANGED_B;
684 } break;
685 default: {
686 /* ERROR: Must not happen */
687 back = DK4_ITER_RESULT_E_CONV;
688 cc = -1;
689 } break;
690 }
691 }
692 }
693 /* Check number of passes
694 */
695 if ((1 == cc) && (0UL < mctx.maxpass) && (passno >= mctx.maxpass)) {
696 back = DK4_ITER_RESULT_E_PASSES;
697 cc = -1;
698 }
699 } while (1 == cc);
700
701 /* Success
702 */
703 if (0 == cc) {
704 *d = x;
705 if (NULL != pp) { *pp = passno; }
706 back = DK4_ITER_RESULT_SUCCESS;
707 }
708
709 finished:
710
711 return back;
712 }
713
714
715
716 /** Run Newton iteration.
717 @param d Destination (address of result variable).
718 @param pp Address of variable to store number of passes on success.
719 @param fct Iteration function, returns two values into the array
720 at address d: the function value and the first derivative
721 value.
722 @param ps Parameter set, may be NULL if fct does not use it.
723 @param x0 Start point.
724 @param ctx Iteration context, may be NULL.
725 @return DK4_ITER_RESULT_SUCCESS on success, one from
726 DK4_ITER_RESULT_E_PASSES, DK4_ITER_RESULT_E_INFINITE,
727 DK4_ITER_RESULT_E_OOR, DK4_ITER_RESULT_E_CONV, DK4_ITER_RESULT_E_FCT,
728 or DK4_ITER_RESULT_E_ARGS on error.
729 */
730 static
731 int
dk4iter_newton(double * d,unsigned long * pp,dk4_iter_fct_t * fct,void const * ps,double x0,dk4_iter_ctx_t const * ctx)732 dk4iter_newton(
733 double *d,
734 unsigned long *pp,
735 dk4_iter_fct_t *fct,
736 void const *ps,
737 double x0,
738 dk4_iter_ctx_t const *ctx
739 )
740 {
741 double v[2];
742 double xn = 0.0;
743 unsigned long passno = 0UL;
744 int back = DK4_ITER_RESULT_E_ARGS;
745 int cc = 1;
746 int res = 0;
747
748 #if DK4_USE_ASSERT
749 assert(NULL != d);
750 assert(NULL != fct);
751 #endif
752 res = (*fct)(v, x0, ps);
753 if (0 == res) {
754 cc = -1;
755 back = DK4_ITER_RESULT_E_FCT;
756 }
757 if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) {
758 cc = -1;
759 back = DK4_ITER_RESULT_E_INFINITE;
760 }
761 while (1 == cc) {
762 if (ULONG_MAX > passno) { passno++; }
763 xn = x0 - v[0] / v[1];
764 if (dk4ma_is_finite(xn)) {
765 if (0 != (MINMAX_MINIMUM & (ctx->minmax))) {
766 if (0 != (MINMAX_MIN_EXCL & (ctx->minmax))) {
767 if (xn <= ctx->xmin) {
768 cc = -1;
769 back = DK4_ITER_RESULT_E_OOR;
770 }
771 }
772 else {
773 if (xn < ctx->xmin) {
774 cc = -1;
775 back = DK4_ITER_RESULT_E_OOR;
776 }
777 }
778 }
779 if (0 != (MINMAX_MAXIMUM & (ctx->minmax))) {
780 if (0 != (MINMAX_MAX_EXCL & (ctx->minmax))) {
781 if (xn >= ctx->xmax) {
782 cc = -1;
783 back = DK4_ITER_RESULT_E_OOR;
784 }
785 }
786 else {
787 if (xn > ctx->xmax) {
788 cc = -1;
789 back = DK4_ITER_RESULT_E_OOR;
790 }
791 }
792 }
793 if (1 == cc) {
794 res = (*fct)(v, xn, ps);
795 if (0 == res) {
796 cc = -1;
797 back = DK4_ITER_RESULT_E_FCT;
798 }
799 if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) {
800 cc = -1;
801 back = DK4_ITER_RESULT_E_INFINITE;
802 }
803 if (1 == cc) {
804
805
806 if (isgreater(ctx->eps_y,0.0)) {
807 if (fabs(v[0]) < ctx->eps_y) {
808 if (xn == x0) {
809 cc = 0;
810 }
811 else {
812 if (0 == ctx->exact) {
813 if(isgreater(ctx->eps_x,0.0)) {
814 if (isless(fabs(xn - x0),ctx->eps_x)) {
815 cc = 0;
816 }
817 }
818 else {
819 cc = 0;
820 }
821 }
822 }
823 }
824 }
825 else {
826 if (xn == x0) {
827 cc = 0;
828 }
829 else {
830 if (0 == ctx->exact) {
831 if (isgreater(ctx->eps_x,0.0)) {
832 if (isless(fabs(xn - x0),ctx->eps_x)) {
833 cc = 0;
834 }
835 }
836 else {
837 /* Neither eps_x nor eps_y,
838 so we wait for x=x0.
839 */
840 }
841 }
842 }
843 }
844 }
845 x0 = xn;
846 }
847 }
848 else {
849 cc = -1;
850 back = DK4_ITER_RESULT_E_INFINITE;
851 }
852 if ((1 == cc) && (0UL < ctx->maxpass) && (passno >= ctx->maxpass)) {
853 cc = -1;
854 back = DK4_ITER_RESULT_E_PASSES;
855 }
856 }
857 if (0 == cc) {
858 *d = xn;
859 if (NULL != pp) { *pp = passno; }
860 back = DK4_ITER_RESULT_SUCCESS;
861 }
862
863 return back;
864 }
865
866
867
868 /** Run fix point iteration.
869 @param d Destination (address of result variable).
870 @param pp Address of variable to store number of passes on success.
871 @param fct Iteration function, the phi part of phi(x)=x.
872 @param ps Parameter set, may be NULL if fct does not use it.
873 @param x0 Start point.
874 @param ctx Iteration context, may be NULL.
875 @return DK4_ITER_RESULT_SUCCESS on success, one from
876 DK4_ITER_RESULT_E_PASSES, DK4_ITER_RESULT_E_INFINITE,
877 DK4_ITER_RESULT_E_OOR, DK4_ITER_RESULT_E_CONV, DK4_ITER_RESULT_E_FCT,
878 or DK4_ITER_RESULT_E_ARGS on error.
879 */
880 static
881 int
dk4iter_fix_point(double * d,unsigned long * pp,dk4_iter_fct_t * fct,void const * ps,double x0,dk4_iter_ctx_t const * ctx)882 dk4iter_fix_point(
883 double *d,
884 unsigned long *pp,
885 dk4_iter_fct_t *fct,
886 void const *ps,
887 double x0,
888 dk4_iter_ctx_t const *ctx
889 )
890 {
891 #if 0
892 double a; /* Border from previous steps */
893 double b; /* Border from previous steps */
894 #endif
895 double xn = 0.0; /* New x value */
896 unsigned long passno = 0UL; /* Iteration step number */
897 int res = 0; /* Function evaluation result */
898 int cc = 1; /* Flag: Can continue */
899 int back = DK4_ITER_RESULT_E_ARGS;
900
901 #if DK4_USE_ASSERT
902 assert(NULL != d);
903 assert(NULL != fct);
904 #endif
905 while (1 == cc) {
906 if (ULONG_MAX > passno) { passno++; }
907 /*
908 Calculate new x, check result
909 */
910 res = (*fct)(&xn, x0, ps);
911 if (0 == res) {
912 cc = -1;
913 back = DK4_ITER_RESULT_E_FCT;
914 }
915 if (!(dk4ma_is_finite(xn))) {
916 cc = -1;
917 back = DK4_ITER_RESULT_E_INFINITE;
918 }
919 if (1 == cc) {
920 /*
921 Check whether specified interval is exceeded
922 */
923 if (0 != (MINMAX_MINIMUM & (ctx->minmax))) {
924 if (0 != (MINMAX_MIN_EXCL & (ctx->minmax))) {
925 if (xn <= ctx->xmin) {
926 cc = -1;
927 back = DK4_ITER_RESULT_E_OOR;
928 }
929 }
930 else {
931 if (xn < ctx->xmin) {
932 cc = -1;
933 back = DK4_ITER_RESULT_E_OOR;
934 }
935 }
936 }
937 if (0 != (MINMAX_MAXIMUM & (ctx->minmax))) {
938 if (0 != (MINMAX_MAX_EXCL & (ctx->minmax))) {
939 if (xn >= ctx->xmax) {
940 cc = -1;
941 back = DK4_ITER_RESULT_E_OOR;
942 }
943 }
944 else {
945 if (xn > ctx->xmax) {
946 cc = -1;
947 back = DK4_ITER_RESULT_E_OOR;
948 }
949 }
950 }
951 if (1 == cc) {
952 if (1 == cc) {
953 /*
954 Check whether we are finished
955 */
956 if (xn == x0) {
957 cc = 0;
958 }
959 else {
960 if ((0 == ctx->exact) && (0.0 < ctx->eps_x)) {
961 if (isless(fabs(xn-x0),ctx->eps_x)) {
962 cc = 0;
963 }
964 #if TRACE_DEBUG
965 else {
966 }
967 #endif
968 }
969 else {
970 }
971 }
972 }
973 }
974 x0 = xn;
975 }
976 /* Stop if too many passes
977 */
978 if (1 == cc) {
979 if ((0UL < ctx->maxpass) && (passno >= ctx->maxpass)) {
980 cc = -1;
981 back = DK4_ITER_RESULT_E_PASSES;
982 }
983 }
984 }
985 if (0 == cc) {
986 *d = xn;
987 if (NULL != pp) { *pp = passno; }
988 back = DK4_ITER_RESULT_SUCCESS;
989 }
990
991 return back;
992 }
993
994
995
996 int
dk4iter_start_point(double * d,unsigned long * pp,dk4_iter_fct_t * fct,void const * ps,double x0,dk4_iter_ctx_t const * ctx)997 dk4iter_start_point(
998 double *d,
999 unsigned long *pp,
1000 dk4_iter_fct_t *fct,
1001 void const *ps,
1002 double x0,
1003 dk4_iter_ctx_t const *ctx
1004 )
1005 {
1006 dk4_iter_ctx_t mctx;
1007 int back = DK4_ITER_RESULT_E_ARGS;
1008
1009 /* Check function call arguments
1010 */
1011 #if DK4_USE_ASSERT
1012 assert(NULL != d);
1013 assert(NULL != fct);
1014 #endif
1015 if ((NULL == d) || (NULL == fct)) {
1016 goto finished;
1017 }
1018 /* Copy or initialize context
1019 */
1020 if (NULL != ctx) {
1021 DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t));
1022 }
1023 else {
1024 dk4iter_ctx_init(&mctx);
1025 mctx.algo = DK4_ITER_ALG_NEWTON;
1026 }
1027 /* Check x0 in interval
1028 */
1029 if (0 != (MINMAX_MINIMUM & (mctx.minmax))) {
1030 if (0 != (MINMAX_MIN_EXCL & (mctx.minmax))) {
1031 if (x0 <= mctx.xmin) {
1032 goto finished;
1033 }
1034 }
1035 else {
1036 if (x0 < mctx.xmin) {
1037 goto finished;
1038 }
1039 }
1040 }
1041 if (0 != (MINMAX_MAXIMUM & (mctx.minmax))) {
1042 if (0 != (MINMAX_MAX_EXCL & (mctx.minmax))) {
1043 if (x0 >= mctx.xmax) {
1044 goto finished;
1045 }
1046 }
1047 else {
1048 if (x0 > mctx.xmax) {
1049 goto finished;
1050 }
1051 }
1052 }
1053 /* Call function for specified algorithm.
1054 */
1055 switch (mctx.algo) {
1056 case DK4_ITER_ALG_NEWTON : {
1057 back = dk4iter_newton(d, pp, fct, ps, x0, &mctx);
1058 } break;
1059 case DK4_ITER_ALG_FIX_POINT : {
1060 back = dk4iter_fix_point(d, pp, fct, ps, x0, &mctx);
1061 } break;
1062 }
1063
1064 finished:
1065 return back;
1066 }
1067
1068
1069
1070
1071 /* vim: set ai sw=4 ts=4 : */
1072
1073