diff options
Diffstat (limited to 'support/dktools/dk4iter.c')
-rw-r--r-- | support/dktools/dk4iter.c | 1226 |
1 files changed, 0 insertions, 1226 deletions
diff --git a/support/dktools/dk4iter.c b/support/dktools/dk4iter.c deleted file mode 100644 index d958ac2b1e..0000000000 --- a/support/dktools/dk4iter.c +++ /dev/null @@ -1,1226 +0,0 @@ -/* -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 : */ - |