summaryrefslogtreecommitdiff
path: root/support/dktools/dk4iter.c
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2020-09-20 03:03:26 +0000
committerNorbert Preining <norbert@preining.info>2020-09-20 03:03:26 +0000
commit1f457376b478257b88d4a857f5ec1b6155442dd7 (patch)
tree2a06a60551dea362cf8cb0cb0ba66c78608717c4 /support/dktools/dk4iter.c
parentac690ca29ad5bf8a5203a65fd6252f7b564f4727 (diff)
CTAN sync 202009200303
Diffstat (limited to 'support/dktools/dk4iter.c')
-rw-r--r--support/dktools/dk4iter.c1226
1 files changed, 1226 insertions, 0 deletions
diff --git a/support/dktools/dk4iter.c b/support/dktools/dk4iter.c
new file mode 100644
index 0000000000..d958ac2b1e
--- /dev/null
+++ b/support/dktools/dk4iter.c
@@ -0,0 +1,1226 @@
+/*
+Copyright (C) 2018-2020, Dirk Krause
+SPDX-License-Identifier: BSD-3-Clause
+*/
+
+/*
+ WARNING: This file was generated by the dkct program (see
+ http://dktools.sourceforge.net/ for details).
+ Changes you make here will be lost if dkct is run again!
+ You should modify the original source and run dkct on it.
+ Original source: dk4iter.ctr
+*/
+
+/** @file dk4iter.c The dk4iter module.
+*/
+
+
+#line 401 "dk4iter.ctr"
+
+
+#include "dk4conf.h"
+
+#include <stdio.h>
+
+#if DK4_HAVE_ASSERT_H
+#ifndef ASSERT_H_INCLUDED
+#include <assert.h>
+#define ASSERT_H_INCLUDED 1
+#endif
+#endif
+
+#if DK4_HAVE_STDLIB_H
+#ifndef STDLIB_H_INCLUDED
+#include <stdlib.h>
+#define STDLIB_H_INCLUDED 1
+#endif
+#endif
+
+#if DK4_HAVE_LIMITS_H
+#ifndef LIMITS_H_INCLUDED
+#include <limits.h>
+#define LIMITS_H_INCLUDED 1
+#endif
+#endif
+
+#if DK4_HAVE_STDINT_H
+#ifndef STDINT_H_INCLUDED
+#include <stdint.h>
+#define STDINT_H_INCLUDED 1
+#endif
+#endif
+
+
+#include "dk4mem.h"
+#include "dk4iter.h"
+#include "dk4math.h"
+
+
+
+
+
+#line 443 "dk4iter.ctr"
+
+
+
+/** Initial function value condition.
+*/
+enum {
+ CONDITION_ILLEGAL = 0 ,
+ /** fa < 0, fb >= 0.
+ */
+ CONDITION_FA_LESS_ZERO ,
+
+ /** fa <= 0, fb > 0.
+ */
+ CONDITION_FA_LEQ_ZERO ,
+
+ /** fa > 0, fb <= 0.
+ */
+ CONDITION_FA_GREATER_ZERO ,
+
+ /** fa >= 0, fb < 0.
+ */
+ CONDITION_FA_GEQ_ZERO
+};
+
+
+
+/** Note which interval border was changed in previous step.
+*/
+enum {
+ /** No interval border change yet.
+ */
+ CHANGED_NONE = 0,
+
+ /** Previous step changed interval border a.
+ */
+ CHANGED_A = -1,
+
+ /** Previous step changed interval border b.
+ */
+ CHANGED_B = 1
+};
+
+
+
+/** Flags which restrictions apply to x values.
+*/
+enum {
+ /** Minimum specified.
+ */
+ MINMAX_MINIMUM = 0x0001 ,
+
+ /** Maximum specified.
+ */
+ MINMAX_MAXIMUM = 0x0002 ,
+
+ /** Specified minimum is exclusive.
+ */
+ MINMAX_MIN_EXCL = 0x0004 ,
+
+ /** Specified maximum is exclusive.
+ */
+ MINMAX_MAX_EXCL = 0x0008
+};
+
+
+
+void
+dk4iter_ctx_init(dk4_iter_ctx_t *ctx)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ DK4_MEMRES(ctx,sizeof(dk4_iter_ctx_t));
+ ctx->xmin = 0.0;
+ ctx->xmax = 0.0;
+ ctx->eps_x = 1.0e-8;
+ ctx->eps_y = 1.0e-8;
+ ctx->maxpass = (unsigned long)(DK4_ITER_PASSES_REGULAR);
+ ctx->exact = 0;
+ ctx->algo = DK4_ITER_ALG_RF_ANDERSON_BJOERCK;
+ ctx->minmax = 0;
+ }
+}
+
+
+
+dk4_iter_ctx_t *
+dk4iter_ctx_open(void)
+{
+ dk4_iter_ctx_t *back = NULL;
+ back = dk4mem_new(dk4_iter_ctx_t,1,NULL);
+ if (NULL != back) {
+ dk4iter_ctx_init(back);
+ }
+ return back;
+}
+
+
+
+void
+dk4iter_ctx_close(dk4_iter_ctx_t *ctx)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ dk4mem_free(ctx);
+ }
+}
+
+
+
+void
+dk4iter_ctx_set_eps_x(dk4_iter_ctx_t *ctx, double eps)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ ctx->eps_x = eps;
+ }
+}
+
+
+
+void
+dk4iter_ctx_set_eps_y(dk4_iter_ctx_t *ctx, double eps)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ ctx->eps_y = eps;
+ }
+}
+
+
+
+void
+dk4iter_ctx_set_maxpass(dk4_iter_ctx_t *ctx, unsigned long passes)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ ctx->maxpass = passes;
+ }
+}
+
+
+
+void
+dk4iter_ctx_set_exact(dk4_iter_ctx_t *ctx, int flag)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ ctx->exact = flag;
+ }
+}
+
+
+
+void
+dk4iter_ctx_set_algorithm(dk4_iter_ctx_t *ctx, int algorithm)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ ctx->algo = algorithm;
+ }
+}
+
+
+
+void
+dk4iter_ctx_set_min(dk4_iter_ctx_t *ctx, double xmin)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ ctx->xmin = xmin;
+ ctx->minmax |= MINMAX_MINIMUM;
+ }
+}
+
+
+
+void
+dk4iter_ctx_set_max(dk4_iter_ctx_t *ctx, double xmax)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ ctx->xmax = xmax;
+ ctx->minmax |= MINMAX_MAXIMUM;
+ }
+}
+
+
+
+void
+dk4iter_ctx_set_exclusive_min(dk4_iter_ctx_t *ctx, double xmin)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ ctx->xmin = xmin;
+ ctx->minmax |= (MINMAX_MINIMUM | MINMAX_MIN_EXCL);
+ }
+}
+
+
+
+void
+dk4iter_ctx_set_exclusive_max(dk4_iter_ctx_t *ctx, double xmax)
+{
+#if DK4_USE_ASSERT
+ assert(NULL != ctx);
+#endif
+ if (NULL != ctx) {
+ ctx->xmax = xmax;
+ ctx->minmax |= (MINMAX_MAXIMUM | MINMAX_MAX_EXCL);
+ }
+}
+
+
+
+static
+int
+find_condition(double fa, double fb)
+{
+ int back = CONDITION_ILLEGAL;
+
+ if ((0.0 > fa) && (0.0 <= fb)) {
+ back = CONDITION_FA_LESS_ZERO;
+ }
+ else {
+ if ((0.0 < fa) && (0.0 >= fb)) {
+ back = CONDITION_FA_GREATER_ZERO;
+ }
+ else {
+ if ((0.0 >= fa) && (0.0 < fb)) {
+ back = CONDITION_FA_LEQ_ZERO;
+ }
+ else {
+ if ((0.0 <= fa) && (0.0 > fb)) {
+ back = CONDITION_FA_GEQ_ZERO;
+ }
+ }
+ }
+ }
+ return back;
+}
+
+
+
+/** Calculate gamma value for PEGASUS variant.
+ @param fg Function value at border last set.
+ @param fx Function value at recent x position.
+ @return Gamma value.
+*/
+static
+double
+dk4iter_gamma_pegasus(double fg, double fx)
+{
+ double back;
+ back = fg / (fg + fx);
+ if (0 == dk4ma_is_finite(back)) {
+ back = 0.5;
+ }
+ else {
+ if (0.0 >= back) {
+ back = 0.5;
+ }
+#if 0
+ /* 2018-07-06
+ Can not happen, fg and fx have same sign.
+ */
+ else {
+ if (1.0 < back) {
+ back = 1.0;
+ }
+ }
+#endif
+ }
+ return back;
+}
+
+
+
+/** Calculate gamma value for ANDERSON-BJOERCK variant.
+ @param fg Function value at border last set.
+ @param fx Function value at recent x position.
+ @return Gamma value.
+*/
+static
+double
+dk4iter_gamma_anderson_bjoerck(double fg, double fx)
+{
+ double back;
+ back = 1.0 - fx / fg;
+ if (0 == dk4ma_is_finite(back)) {
+ back = 0.5;
+ }
+ else {
+ if (0.0 >= back) {
+ back = 0.5;
+ }
+#if 0
+ /* 2018-07-06
+ Can not happen as fx and fy have same sign.
+ */
+ else {
+ if (1.0 < back) {
+ back = 1.0;
+ }
+ }
+#endif
+ }
+ return back;
+}
+
+
+
+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
+)
+{
+ dk4_iter_ctx_t mctx; /* Copy of context */
+ double fa = 0.0; /* Function value for border a */
+ double fb = 0.0; /* Function value for border b */
+ double x = 0.0; /* Text x value */
+ double fx = 0.0; /* Function value for x */
+ double afx = 0.0; /* Absolute function value for x */
+ double xo = 0.0; /* Previous step x value */
+ double gamma = 0.5; /* Correction factor */
+ unsigned long passno = 0UL; /* Number of current pass */
+ int res = 0; /* Operation result */
+ int cond = 0; /* Condition of interval borders */
+ int cc = 0; /* 1=continue, 0=finished, -1=abort */
+ int pc = CHANGED_NONE; /* Previous border change */
+ int nc = CHANGED_NONE; /* Next change */
+ int back = DK4_ITER_RESULT_E_ARGS;
+
+
+#line 801 "dk4iter.ctr"
+ /* Check function call arguments
+ */
+#if DK4_USE_ASSERT
+ assert(NULL != d);
+ assert(NULL != fct);
+#endif
+ if ((NULL == d) || (NULL == fct)) {
+
+#line 808 "dk4iter.ctr"
+ goto finished;
+ }
+ if (a == b) {
+
+#line 811 "dk4iter.ctr"
+ goto finished;
+ }
+
+ /* Copy or initialize context
+ */
+ if (NULL != ctx) {
+
+#line 817 "dk4iter.ctr"
+ DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t));
+ }
+ else {
+
+#line 820 "dk4iter.ctr"
+ dk4iter_ctx_init(&mctx);
+ }
+
+ /* Check algorithm specification
+ */
+ if (DK4_ITER_ALG_BISECTION != mctx.algo) {
+ if (DK4_ITER_ALG_RF_PRIMITIVE != mctx.algo) {
+ if (DK4_ITER_ALG_RF_ILLINOIS != mctx.algo) {
+ if (DK4_ITER_ALG_RF_PEGASUS != mctx.algo) {
+ if (DK4_ITER_ALG_RF_ANDERSON_BJOERCK != mctx.algo) {
+ goto finished;
+
+#line 831 "dk4iter.ctr"
+ }
+ }
+ }
+ }
+ }
+
+ /* Calculate borders at beginning, check
+ */
+ res = (*fct)(&fa, a, ps);
+ if ((0 == res) || (!(dk4ma_is_finite(fa)))) {
+ back = DK4_ITER_RESULT_E_FCT;
+
+#line 842 "dk4iter.ctr"
+ goto finished;
+ }
+ res = (*fct)(&fb, b, ps);
+ if ((0 == res) || (!(dk4ma_is_finite(fb)))) {
+ back = DK4_ITER_RESULT_E_FCT;
+
+#line 847 "dk4iter.ctr"
+ goto finished;
+ }
+ cond = find_condition(fa, fb);
+ if (CONDITION_ILLEGAL == cond) {
+
+#line 851 "dk4iter.ctr"
+ goto finished;
+ }
+
+ /* Prepare iteration loop
+ */
+ passno = 0UL;
+ x = xo = 0.0;
+ cc = 1;
+ pc = CHANGED_NONE;
+
+ /* Run iteration loop
+ */
+ do {
+ /* Keep previous x position
+ */
+ xo = x;
+ /* Increase pass number, but avoid wrapping
+ */
+ if (ULONG_MAX > passno) { passno++; }
+
+#line 870 "dk4iter.ctr"
+ /* Calculate x position
+ */
+
+
+#line 873 "dk4iter.ctr"
+ switch (mctx.algo) {
+ case DK4_ITER_ALG_RF_PRIMITIVE :
+ case DK4_ITER_ALG_RF_ILLINOIS :
+ case DK4_ITER_ALG_RF_PEGASUS :
+ case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : { /* Regula falsi */
+ x = (a * fb - b * fa) / (fb - fa);
+ } break;
+ default : { /* Bisection */
+ x = (a + b) / 2.0;
+ } break;
+ }
+ if (0 == dk4ma_is_finite(x)) {
+ back = DK4_ITER_RESULT_E_INFINITE;
+ cc = -1;
+ }
+ else {
+
+#line 889 "dk4iter.ctr"
+ /* Calculate function for x, if x is usable
+ */
+ res = (*fct)(&fx, x, ps);
+ if (0 == res) {
+
+#line 893 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_FCT;
+ cc = -1;
+ }
+ if (0 == dk4ma_is_finite(fx)) {
+
+#line 897 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_INFINITE;
+ cc = -1;
+ }
+ afx = fabs(fx);
+ if (0 == dk4ma_is_finite(afx)) {
+
+#line 902 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_INFINITE;
+ cc = -1;
+ }
+ /* Check whether we are done, if we could continue
+ */
+ if (1 == cc) {
+
+#line 908 "dk4iter.ctr"
+
+
+#line 909 "dk4iter.ctr"
+
+
+#line 910 "dk4iter.ctr"
+
+
+#line 911 "dk4iter.ctr"
+ if (isgreater(mctx.eps_y,0.0)) {
+ /* Must check y */
+ if (afx < mctx.eps_y) {
+
+#line 914 "dk4iter.ctr"
+ if ((x == xo) && (1UL < passno)) {
+ cc = 0;
+
+#line 916 "dk4iter.ctr"
+ }
+ else {
+ if (0 == mctx.exact) {
+ if(isgreater(mctx.eps_x,0.0)) {
+ if (isless(fabs(x - xo),mctx.eps_x)) {
+ if (1UL < passno) {
+ cc = 0;
+
+#line 923 "dk4iter.ctr"
+ }
+ }
+ }
+ else {
+ cc = 0;
+
+#line 928 "dk4iter.ctr"
+ }
+ }
+ }
+ }
+ }
+ else {
+
+#line 934 "dk4iter.ctr"
+ /* No y check */
+ if ((x == xo) && (1UL < passno)) {
+ cc = 0;
+
+#line 937 "dk4iter.ctr"
+ }
+ else {
+ if (0 == mctx.exact) {
+ if(isgreater(mctx.eps_x,0.0)) {
+ if (isless(fabs(x - xo),mctx.eps_x)) {
+ if (1UL < passno) {
+ cc = 0;
+
+#line 944 "dk4iter.ctr"
+ }
+ }
+ }
+ else {
+ /* Neither eps_x nor eps_y,
+ so we wait for x=xo.
+ */
+ }
+ }
+ }
+ }
+ }
+ /* Continue for usable y values only
+ */
+ if (1 == cc) {
+ /* Find direction for next change
+ */
+ nc = CHANGED_NONE;
+ switch (cond) {
+ case CONDITION_FA_LESS_ZERO : {
+ if (0.0 > fx) {
+ nc = CHANGED_A;
+ }
+ else {
+ nc = CHANGED_B;
+ }
+ } break;
+ case CONDITION_FA_LEQ_ZERO : {
+ if (0.0 >= fx) {
+ nc = CHANGED_A;
+ }
+ else {
+ nc = CHANGED_B;
+ }
+ } break;
+ case CONDITION_FA_GREATER_ZERO : {
+ if (0.0 < fx) {
+ nc = CHANGED_A;
+ }
+ else {
+ nc = CHANGED_B;
+ }
+ } break;
+ case CONDITION_FA_GEQ_ZERO : {
+ if (0.0 <= fx) {
+ nc = CHANGED_A;
+ }
+ else {
+ nc = CHANGED_B;
+ }
+ } break;
+ }
+ /* Apply interval border change
+ */
+ switch (nc) {
+ case CHANGED_A: {
+
+#line 1000 "dk4iter.ctr"
+ if (CHANGED_A == pc) {
+ switch (mctx.algo) {
+ case DK4_ITER_ALG_RF_ILLINOIS : {
+ fb = 0.5 * fb;
+ if (0 == dk4ma_is_finite(fb)) {
+ back = DK4_ITER_RESULT_E_INFINITE;
+ cc = -1;
+
+#line 1007 "dk4iter.ctr"
+ }
+ } break;
+ case DK4_ITER_ALG_RF_PEGASUS : {
+ gamma = dk4iter_gamma_pegasus(fa, fx);
+ fb *= gamma;
+ if (0 == dk4ma_is_finite(fb)) {
+ back = DK4_ITER_RESULT_E_INFINITE;
+ cc = -1;
+
+#line 1015 "dk4iter.ctr"
+ }
+ } break;
+ case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : {
+ gamma = dk4iter_gamma_anderson_bjoerck(
+ fa, fx
+ );
+ fb *= gamma;
+ if (0 == dk4ma_is_finite(fb)) {
+ back = DK4_ITER_RESULT_E_INFINITE;
+ cc = -1;
+
+#line 1025 "dk4iter.ctr"
+ }
+ } break;
+ }
+ }
+ a = x;
+ fa = fx;
+ pc = CHANGED_A;
+ } break;
+ case CHANGED_B: {
+
+#line 1034 "dk4iter.ctr"
+ if (CHANGED_B == pc) {
+ switch (mctx.algo) {
+ case DK4_ITER_ALG_RF_ILLINOIS : {
+ fa = 0.5 * fa;
+ if (0 == dk4ma_is_finite(fa)) {
+ back = DK4_ITER_RESULT_E_INFINITE;
+ cc = -1;
+
+#line 1041 "dk4iter.ctr"
+ }
+ } break;
+ case DK4_ITER_ALG_RF_PEGASUS : {
+ gamma = dk4iter_gamma_pegasus(
+ fb, fx
+ );
+ fa *= gamma;
+ if (0 == dk4ma_is_finite(fa)) {
+ back = DK4_ITER_RESULT_E_INFINITE;
+ cc = -1;
+
+#line 1051 "dk4iter.ctr"
+ }
+ } break;
+ case DK4_ITER_ALG_RF_ANDERSON_BJOERCK : {
+ gamma = dk4iter_gamma_anderson_bjoerck(
+ fb, fx
+ );
+ fa *= gamma;
+ if (0 == dk4ma_is_finite(fa)) {
+ back = DK4_ITER_RESULT_E_INFINITE;
+ cc = -1;
+
+#line 1061 "dk4iter.ctr"
+ }
+ } break;
+ }
+ }
+ b = x;
+ fb = fx;
+ pc = CHANGED_B;
+ } break;
+ default: {
+ /* ERROR: Must not happen */
+ back = DK4_ITER_RESULT_E_CONV;
+
+#line 1072 "dk4iter.ctr"
+ cc = -1;
+ } break;
+ }
+ }
+ }
+ /* Check number of passes
+ */
+ if ((1 == cc) && (0UL < mctx.maxpass) && (passno >= mctx.maxpass)) {
+ back = DK4_ITER_RESULT_E_PASSES;
+
+#line 1081 "dk4iter.ctr"
+ cc = -1;
+ }
+ } while (1 == cc);
+
+ /* Success
+ */
+ if (0 == cc) {
+ *d = x;
+ if (NULL != pp) { *pp = passno; }
+ back = DK4_ITER_RESULT_SUCCESS;
+ }
+
+ finished:
+
+
+#line 1095 "dk4iter.ctr"
+ return back;
+}
+
+
+
+/** Run Newton iteration.
+ @param d Destination (address of result variable).
+ @param pp Address of variable to store number of passes on success.
+ @param fct Iteration function, returns two values into the array
+ at address d: the function value and the first derivative
+ value.
+ @param ps Parameter set, may be NULL if fct does not use it.
+ @param x0 Start point.
+ @param ctx Iteration context, may be NULL.
+ @return DK4_ITER_RESULT_SUCCESS on success, one from
+ DK4_ITER_RESULT_E_PASSES, DK4_ITER_RESULT_E_INFINITE,
+ DK4_ITER_RESULT_E_OOR, DK4_ITER_RESULT_E_CONV, DK4_ITER_RESULT_E_FCT,
+ or DK4_ITER_RESULT_E_ARGS on error.
+*/
+static
+int
+dk4iter_newton(
+ double *d,
+ unsigned long *pp,
+ dk4_iter_fct_t *fct,
+ void const *ps,
+ double x0,
+ dk4_iter_ctx_t const *ctx
+)
+{
+ double v[2];
+ double xn = 0.0;
+ unsigned long passno = 0UL;
+ int back = DK4_ITER_RESULT_E_ARGS;
+ int cc = 1;
+ int res = 0;
+
+
+#line 1132 "dk4iter.ctr"
+#if DK4_USE_ASSERT
+ assert(NULL != d);
+ assert(NULL != fct);
+#endif
+ res = (*fct)(v, x0, ps);
+ if (0 == res) {
+ cc = -1;
+
+#line 1139 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_FCT;
+ }
+ if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) {
+ cc = -1;
+
+#line 1143 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_INFINITE;
+ }
+ while (1 == cc) {
+ if (ULONG_MAX > passno) { passno++; }
+
+#line 1147 "dk4iter.ctr"
+ xn = x0 - v[0] / v[1];
+ if (dk4ma_is_finite(xn)) {
+
+#line 1149 "dk4iter.ctr"
+ if (0 != (MINMAX_MINIMUM & (ctx->minmax))) {
+ if (0 != (MINMAX_MIN_EXCL & (ctx->minmax))) {
+ if (xn <= ctx->xmin) {
+ cc = -1;
+
+#line 1153 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_OOR;
+ }
+ }
+ else {
+ if (xn < ctx->xmin) {
+ cc = -1;
+
+#line 1159 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_OOR;
+ }
+ }
+ }
+ if (0 != (MINMAX_MAXIMUM & (ctx->minmax))) {
+ if (0 != (MINMAX_MAX_EXCL & (ctx->minmax))) {
+ if (xn >= ctx->xmax) {
+ cc = -1;
+
+#line 1167 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_OOR;
+ }
+ }
+ else {
+ if (xn > ctx->xmax) {
+ cc = -1;
+
+#line 1173 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_OOR;
+ }
+ }
+ }
+ if (1 == cc) {
+ res = (*fct)(v, xn, ps);
+ if (0 == res) {
+ cc = -1;
+
+#line 1181 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_FCT;
+ }
+ if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) {
+ cc = -1;
+
+#line 1185 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_INFINITE;
+ }
+ if (1 == cc) {
+
+
+#line 1189 "dk4iter.ctr"
+
+
+#line 1190 "dk4iter.ctr"
+ if (isgreater(ctx->eps_y,0.0)) {
+
+#line 1191 "dk4iter.ctr"
+ if (fabs(v[0]) < ctx->eps_y) {
+
+#line 1192 "dk4iter.ctr"
+ if (xn == x0) {
+ cc = 0;
+
+#line 1194 "dk4iter.ctr"
+ }
+ else {
+ if (0 == ctx->exact) {
+
+#line 1197 "dk4iter.ctr"
+ if(isgreater(ctx->eps_x,0.0)) {
+ if (isless(fabs(xn - x0),ctx->eps_x)) {
+ cc = 0;
+
+#line 1200 "dk4iter.ctr"
+ }
+ }
+ else {
+ cc = 0;
+
+#line 1204 "dk4iter.ctr"
+ }
+ }
+ }
+ }
+ }
+ else {
+
+#line 1210 "dk4iter.ctr"
+ if (xn == x0) {
+ cc = 0;
+
+#line 1212 "dk4iter.ctr"
+ }
+ else {
+ if (0 == ctx->exact) {
+ if (isgreater(ctx->eps_x,0.0)) {
+ if (isless(fabs(xn - x0),ctx->eps_x)) {
+ cc = 0;
+
+#line 1218 "dk4iter.ctr"
+ }
+ }
+ else {
+ /* Neither eps_x nor eps_y,
+ so we wait for x=x0.
+ */
+ }
+ }
+ }
+ }
+ }
+ x0 = xn;
+ }
+ }
+ else {
+ cc = -1;
+
+#line 1234 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_INFINITE;
+ }
+ if ((1 == cc) && (0UL < ctx->maxpass) && (passno >= ctx->maxpass)) {
+ cc = -1;
+ back = DK4_ITER_RESULT_E_PASSES;
+ }
+ }
+ if (0 == cc) {
+ *d = xn;
+ if (NULL != pp) { *pp = passno; }
+ back = DK4_ITER_RESULT_SUCCESS;
+ }
+
+
+#line 1247 "dk4iter.ctr"
+ return back;
+}
+
+
+
+/** Run fix point iteration.
+ @param d Destination (address of result variable).
+ @param pp Address of variable to store number of passes on success.
+ @param fct Iteration function, the phi part of phi(x)=x.
+ @param ps Parameter set, may be NULL if fct does not use it.
+ @param x0 Start point.
+ @param ctx Iteration context, may be NULL.
+ @return DK4_ITER_RESULT_SUCCESS on success, one from
+ DK4_ITER_RESULT_E_PASSES, DK4_ITER_RESULT_E_INFINITE,
+ DK4_ITER_RESULT_E_OOR, DK4_ITER_RESULT_E_CONV, DK4_ITER_RESULT_E_FCT,
+ or DK4_ITER_RESULT_E_ARGS on error.
+*/
+static
+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
+)
+{
+#if 0
+ double a; /* Border from previous steps */
+ double b; /* Border from previous steps */
+#endif
+ double xn = 0.0; /* New x value */
+ unsigned long passno = 0UL; /* Iteration step number */
+ int res = 0; /* Function evaluation result */
+ int cc = 1; /* Flag: Can continue */
+ int back = DK4_ITER_RESULT_E_ARGS;
+
+
+#line 1285 "dk4iter.ctr"
+#if DK4_USE_ASSERT
+ assert(NULL != d);
+ assert(NULL != fct);
+#endif
+ while (1 == cc) {
+ if (ULONG_MAX > passno) { passno++; }
+
+#line 1291 "dk4iter.ctr"
+ /*
+ Calculate new x, check result
+ */
+ res = (*fct)(&xn, x0, ps);
+ if (0 == res) {
+ cc = -1;
+
+#line 1297 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_FCT;
+ }
+ if (!(dk4ma_is_finite(xn))) {
+ cc = -1;
+
+#line 1301 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_INFINITE;
+ }
+ if (1 == cc) {
+
+#line 1304 "dk4iter.ctr"
+ /*
+ Check whether specified interval is exceeded
+ */
+ if (0 != (MINMAX_MINIMUM & (ctx->minmax))) {
+ if (0 != (MINMAX_MIN_EXCL & (ctx->minmax))) {
+ if (xn <= ctx->xmin) {
+ cc = -1;
+
+#line 1311 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_OOR;
+ }
+ }
+ else {
+ if (xn < ctx->xmin) {
+ cc = -1;
+
+#line 1317 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_OOR;
+ }
+ }
+ }
+ if (0 != (MINMAX_MAXIMUM & (ctx->minmax))) {
+ if (0 != (MINMAX_MAX_EXCL & (ctx->minmax))) {
+ if (xn >= ctx->xmax) {
+ cc = -1;
+
+#line 1325 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_OOR;
+ }
+ }
+ else {
+ if (xn > ctx->xmax) {
+ cc = -1;
+
+#line 1331 "dk4iter.ctr"
+ back = DK4_ITER_RESULT_E_OOR;
+ }
+ }
+ }
+ if (1 == cc) {
+ if (1 == cc) {
+ /*
+ Check whether we are finished
+ */
+ if (xn == x0) {
+ cc = 0;
+
+#line 1342 "dk4iter.ctr"
+ }
+ else {
+ if ((0 == ctx->exact) && (0.0 < ctx->eps_x)) {
+ if (isless(fabs(xn-x0),ctx->eps_x)) {
+ cc = 0;
+
+#line 1347 "dk4iter.ctr"
+ }
+#if TRACE_DEBUG
+ else {
+
+#line 1350 "dk4iter.ctr"
+ }
+#endif
+ }
+ else {
+
+#line 1354 "dk4iter.ctr"
+ }
+ }
+ }
+ }
+ x0 = xn;
+ }
+ /* Stop if too many passes
+ */
+ if (1 == cc) {
+ if ((0UL < ctx->maxpass) && (passno >= ctx->maxpass)) {
+ cc = -1;
+ back = DK4_ITER_RESULT_E_PASSES;
+ }
+ }
+ }
+ if (0 == cc) {
+ *d = xn;
+ if (NULL != pp) { *pp = passno; }
+ back = DK4_ITER_RESULT_SUCCESS;
+ }
+
+
+#line 1375 "dk4iter.ctr"
+ return back;
+}
+
+
+
+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
+)
+{
+ dk4_iter_ctx_t mctx;
+ int back = DK4_ITER_RESULT_E_ARGS;
+
+ /* Check function call arguments
+ */
+#if DK4_USE_ASSERT
+ assert(NULL != d);
+ assert(NULL != fct);
+#endif
+ if ((NULL == d) || (NULL == fct)) {
+ goto finished;
+ }
+ /* Copy or initialize context
+ */
+ if (NULL != ctx) {
+ DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t));
+ }
+ else {
+ dk4iter_ctx_init(&mctx);
+ mctx.algo = DK4_ITER_ALG_NEWTON;
+ }
+ /* Check x0 in interval
+ */
+ if (0 != (MINMAX_MINIMUM & (mctx.minmax))) {
+ if (0 != (MINMAX_MIN_EXCL & (mctx.minmax))) {
+ if (x0 <= mctx.xmin) {
+ goto finished;
+ }
+ }
+ else {
+ if (x0 < mctx.xmin) {
+ goto finished;
+ }
+ }
+ }
+ if (0 != (MINMAX_MAXIMUM & (mctx.minmax))) {
+ if (0 != (MINMAX_MAX_EXCL & (mctx.minmax))) {
+ if (x0 >= mctx.xmax) {
+ goto finished;
+ }
+ }
+ else {
+ if (x0 > mctx.xmax) {
+ goto finished;
+ }
+ }
+ }
+ /* Call function for specified algorithm.
+ */
+ switch (mctx.algo) {
+ case DK4_ITER_ALG_NEWTON : {
+ back = dk4iter_newton(d, pp, fct, ps, x0, &mctx);
+ } break;
+ case DK4_ITER_ALG_FIX_POINT : {
+ back = dk4iter_fix_point(d, pp, fct, ps, x0, &mctx);
+ } break;
+ }
+
+ finished:
+ return back;
+}
+
+
+
+
+/* vim: set ai sw=4 ts=4 : */
+