diff options
author | Norbert Preining <norbert@preining.info> | 2020-09-20 03:03:26 +0000 |
---|---|---|
committer | Norbert Preining <norbert@preining.info> | 2020-09-20 03:03:26 +0000 |
commit | 1f457376b478257b88d4a857f5ec1b6155442dd7 (patch) | |
tree | 2a06a60551dea362cf8cb0cb0ba66c78608717c4 /support/dktools/dk4iter.c | |
parent | ac690ca29ad5bf8a5203a65fd6252f7b564f4727 (diff) |
CTAN sync 202009200303
Diffstat (limited to 'support/dktools/dk4iter.c')
-rw-r--r-- | support/dktools/dk4iter.c | 1226 |
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 : */ + |