/* 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 #if DK4_HAVE_ASSERT_H #ifndef ASSERT_H_INCLUDED #include #define ASSERT_H_INCLUDED 1 #endif #endif #if DK4_HAVE_STDLIB_H #ifndef STDLIB_H_INCLUDED #include #define STDLIB_H_INCLUDED 1 #endif #endif #if DK4_HAVE_LIMITS_H #ifndef LIMITS_H_INCLUDED #include #define LIMITS_H_INCLUDED 1 #endif #endif #if DK4_HAVE_STDINT_H #ifndef STDINT_H_INCLUDED #include #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 : */