diff options
author | Norbert Preining <norbert@preining.info> | 2020-10-12 03:04:00 +0000 |
---|---|---|
committer | Norbert Preining <norbert@preining.info> | 2020-10-12 03:04:00 +0000 |
commit | 0ce40abb18ec02ec6fd6bcc5f21612c88daa7578 (patch) | |
tree | 416289fe1448873fd8ca33051f50ad85bffa8aaa /support/dktools/dk4iter.ctr | |
parent | fdb18507cd80dc17f5a5256153d34668b4f4e61c (diff) |
CTAN sync 202010120303
Diffstat (limited to 'support/dktools/dk4iter.ctr')
-rw-r--r-- | support/dktools/dk4iter.ctr | 1456 |
1 files changed, 0 insertions, 1456 deletions
diff --git a/support/dktools/dk4iter.ctr b/support/dktools/dk4iter.ctr deleted file mode 100644 index 1cd125b366..0000000000 --- a/support/dktools/dk4iter.ctr +++ /dev/null @@ -1,1456 +0,0 @@ -%% options - -copyright owner = Dirk Krause -copyright year = 2018-xxxx -SPDX-License-Identifier: BSD-3-Clause - -%% header - -/** @file dk4iter.h Iteration algorithms for root finding. - -This module implements the following root finding iteration algorithms: -- Bisection -- Regula falsi (primitive form, Illinois, Pegasus, Anderson-Bjoerck) -- Newton -- fix point - -The function to iterate must be implemented as dk4_iter_fct_t. -This function type returns an integer value (non-zero to indicate a -successful calculation, 0 to indicate an error). -The function expects the following arguments: -- Result address<br> - Address of a variable or an array to store the calculation result. - Functions for Newton iteration algorithm store 2 values: Function value - and value of the first derivative. -- X position<br> - The x value you want a function value for. -- Address of parameter set. - Additional parameters probably required by the function. I.e. for a - polynomial calculation function you may specify the coefficients here. - This parameter is optional. - -Details for an iteration may be specified in an iteration context. -The dk4iter_ctx_open() creates such a context and returns a pointer. -Use dk4iter_ctx_close() to release the context when done with it. -Alternatively use a automatic/static variable of the dk4_iter_ctx_t type -and initialize it using the dk4iter_ctx_init() function. - -The dk4iter_ctx_set_algorithm() function chooses the iteration method. - -For Newton and fixpoint the context may specify an x interval the -iteration must not leave. Leaving the interval results in abort. -The dk4iter_ctx_set_min(), dk4iter_ctx_set_exclusive_min(), -dk4iter_ctx_set_max() and dk4iter_ctx_set_exclusive_max() functions can -be used to set closed and open interval borders. - -Tolerance values (epsilon) may be specified for -- y direction unless fixpoint is used and/or -- x direction. -Use dk4iter_ctx_set_eps_y() and dk4iter_ctx_set_eps_x() to set the -tolerances. -Use 0.0 or negative tolerances to skip one check (not both!). - -Once the iteration reached allowed y difference you can decide -to stop when the x difference is in allowed tolerance range too or -to continue until full machine precision is reached (new x is exactly -the same as previous x). -I do not recommend attempts to iterate to full machine precision as this -(a) may fail for some functions resulting in an oscillation and -(b) uses a larger number of iteration passes and. -Use dk4iter_ctx_set_exact() to control this. - -The dk4iter_ctx_set_maxpass() function sets the maximum number of -iteration passes (iteration steps). The iteration is aborted if there -is no success within this maximum number of passes. -Although you can use 0 to set an unlimited number of passes I do not recommend -to do so. - -The dk4iter_interval() function runs an iteration for a specified interval. - -The dk4iter_start_point() function runs an iteration if just one -starting x value is specified. - -If no context is specified the following defaults are used: - -Option | Default -:----: | :------ -Algorithm | DK4_ITER_ALG_RF_ANDERSON_BJOERCK for dk4iter_interval(), DK4_ITER_ALG_NEWTON for dk4iter_start_point(). -Max passes | 256. -Y tolerance | 1.0e-8 -X tolerance | 1.0e-8 -Full machine precision | No. -Restricted x interval | None. - -*/ - - - -/** Border values for maximum number of iteration passes. -*/ -enum { - /** Number of passes for average cases. - */ - DK4_ITER_PASSES_REGULAR = 256 , - - /** Do not set a limit on the number of passes. - Warning: Might result in an endless loop! - */ - DK4_ITER_PASSES_UNLIMITED = 0 -}; - - - -/** Iteration algorithms. -*/ -enum { - /** Interval bisection. - */ - DK4_ITER_ALG_BISECTION = 0 , - - /** Primitive from of regula falsi. - */ - DK4_ITER_ALG_RF_PRIMITIVE , - - /** Regula falsi, Illinois variant. - */ - DK4_ITER_ALG_RF_ILLINOIS , - - /** Regula falsi, Pegasus variant. - */ - DK4_ITER_ALG_RF_PEGASUS , - - /** Regula falsi, Anderson-Bjoerck variant. - */ - DK4_ITER_ALG_RF_ANDERSON_BJOERCK , - - /** Newton algorithm. - The function must place 2 elements - in the array specified by pointer: - Function value and derivative value. - */ - DK4_ITER_ALG_NEWTON , - - /** Fix point algorithm. - The function must be the phi(x) - function from the - phi(x) = x equation. - */ - DK4_ITER_ALG_FIX_POINT -}; - - - -/** Iteration result. -*/ -enum { - /** Iteration succeeded. - */ - DK4_ITER_RESULT_SUCCESS = 1 , - - /** Error: Too many passes without success. - */ - DK4_ITER_RESULT_E_PASSES = 0 , - - /** Error: Infinite value or NaN in calculation. - */ - DK4_ITER_RESULT_E_INFINITE = -1 , - - /** Error: X left initial interval. - Can happen with Newton and fix point - algorithm. - */ - DK4_ITER_RESULT_E_OOR = -2 , - - /** Error: Not converging. - Fix point algorithm: Interval was enlarging. - */ - DK4_ITER_RESULT_E_CONV = -3 , - - /** Function calculation failed. - */ - DK4_ITER_RESULT_E_FCT = -4 , - - /** Invalid arguments passed to function call. - */ - DK4_ITER_RESULT_E_ARGS = -5 -}; - - - -/** Function to iterate. - @param d Destination address. One element for all algorithms - except DK4_ITER_ALG_NEWTON which saves function value and derivative - value in a 2 elements array. - @param x X position to calculate the function value for. - @param ps Parameter set, may be NULL. - @return Non-zero value on success, 0 on error. -*/ -typedef int dk4_iter_fct_t(double *d, double x, const void *ps); - - -/** Iteration context. -*/ -typedef struct { - double xmin; /**< Minimum of x interval. */ - double xmax; /**< Maximum of x interval. */ - double eps_x; /**< Epsilon for x change. */ - double eps_y; /**< Epsilon for absolute y value. */ - unsigned long maxpass; /**< Maximum number of passes. */ - int exact; /**< Flag: End only if there is no x change. */ - int algo; /**< Iteration algorithm to use. */ - int minmax; /**< Flags: min(1), max(2) set. */ -} dk4_iter_ctx_t; - - - -#ifdef __cplusplus -extern "C" { -#endif - -/** Open new iteration context. - Allocate memory for new iteration context, set up default value. - A context allocated by this function must be released after - use by the dk4iter_ctx_close() function. - @return Valid pointer to new context on success, NULL on error. -*/ - -dk4_iter_ctx_t * -dk4iter_ctx_open(void); - - -/** Initialize a context to default values. - @param ctx Context to initialize. -*/ - -void -dk4iter_ctx_init(dk4_iter_ctx_t *ctx); - - -/** Close iteration context. - Release memory assigned to the context. - @param ctx Context created by dk4iter_ctx_open(). -*/ - -void -dk4iter_ctx_close(dk4_iter_ctx_t *ctx); - - -/** Set epsilon value for x change. - @param ctx Context to modify. - @param eps Epsilon value for x change. -*/ - -void -dk4iter_ctx_set_eps_x(dk4_iter_ctx_t *ctx, double eps); - - -/** Set maximum absolute y value. - @param ctx Context to modify. - @param eps Maximum absolute y value. -*/ - -void -dk4iter_ctx_set_eps_y(dk4_iter_ctx_t *ctx, double eps); - - -/** Set maximum number of passes (iteration steps). - @param ctx Context to modify. - @param passes Maximum number of passes, 0 or negative values indicate - an unlimited number of steps. -*/ - -void -dk4iter_ctx_set_maxpass(dk4_iter_ctx_t *ctx, unsigned long passes); - - -/** Set flag for exact iteration. - If the flag is activated, iteration is continued until there is - no longer any change in the x value. For some functions you can - retrieve a result in machine precision using this flag, for other - functions the iteration may fail. - @param ctx Context to modify. - @param flag New flag value, 0=inactive, other=active. - The recommended value is 0. -*/ - -void -dk4iter_ctx_set_exact(dk4_iter_ctx_t *ctx, int flag); - - -/** Set up iteration algorithm. - @param ctx Context to modify. - @param algorithm Algorithm to use. -*/ - -void -dk4iter_ctx_set_algorithm(dk4_iter_ctx_t *ctx, int algorithm); - - -/** Set interval minimum. - The allowed interval is only used for Newton and fixed point - algorithm. - @param ctx Context to modify. - @param xmin Minimum x value allowed. -*/ - -void -dk4iter_ctx_set_min(dk4_iter_ctx_t *ctx, double xmin); - - -/** Set exlusive interval minimum (open interval border). - The allowed interval is only used for Newton and fixed point - algorithm. - @param ctx Context to modify. - @param xmin Minimum x value allowed. -*/ - -void -dk4iter_ctx_set_exclusive_min(dk4_iter_ctx_t *ctx, double xmin); - - -/** Set interval maximum. - The allowed interval is only used for Newton and fixed point - algorithm. - @param ctx Context to modify. - @param xmax Maximum x value allowed. -*/ - -void -dk4iter_ctx_set_max(dk4_iter_ctx_t *ctx, double xmax); - - -/** Set exclusive interval maximum (open interval border). - The allowed interval is only used for Newton and fixed point - algorithm. - @param ctx Context to modify. - @param xmax Maximum x value allowed. -*/ - -void -dk4iter_ctx_set_exclusive_max(dk4_iter_ctx_t *ctx, double xmax); - - -/** Run iteration. - The algorithm in the context must be one from: - DK4_ITER_ALG_BISECTION, DK4_ITER_ALG_RF_PRIMITIVE, DK4_ITER_ALG_RF_ILLINOIS, - DK4_ITER_ALG_RF_PEGASUS DK4_ITER_ALG_RF_ANDERSON_BJOERCK. - Without a context DK4_ITER_ALG_RF_ANDERSON_BJOERCK is used as default. - @param d Address of result variable. - @param pp Address of variable to store number of passes, may be NULL. - @param fct Function to find root for. - @param ps Parameter set for function, may be NULL. - @param a One interval border. - @param b Other interval border. - @param ctx Iteration context. - @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 or - DK4_ITER_RESULT_E_FCT on error. -*/ - -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 -); - - -/** Run iteration. - The algorithm in the context must be one from: - DK4_ITER_ALG_NEWTON, DK4_ITER_ALG_FIX_POINT. - Without a context, DK4_ITER_ALG_NEWTON is used as default. - @param d Address of result variable. - @param pp Address of variable to store number of passes, may be NULL. - @param fct Function to find root for. - For Newton iteration this function must set two values - in the destination array: function value and derivative value. - For fixpoint iteration the function calculates the - phi(x) part of phi(x)=x. - @param ps Parameter set for function, may be NULL. - @param x0 Start point. - @param ctx Iteration context. - @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 or - DK4_ITER_RESULT_E_FCT on error. -*/ - -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 -); - - -#ifdef __cplusplus -} -#endif - -/* vim: set ai sw=4 ts=4 : */ - -%% module - - -#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" - - - -$!trace-include - - - -/** 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; - $? "+ dk4iter_interval" - /* Check function call arguments - */ -#if DK4_USE_ASSERT - assert(NULL != d); - assert(NULL != fct); -#endif - if ((NULL == d) || (NULL == fct)) { $? "! d or fct" - goto finished; - } - if (a == b) { $? "! zero length interval" - goto finished; - } - - /* Copy or initialize context - */ - if (NULL != ctx) { $? ". use context" - DK4_MEMCPY(&mctx,ctx,sizeof(dk4_iter_ctx_t)); - } - else { $? ". no context, use defaults" - 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; $? "! illegal algorithm" - } - } - } - } - } - - /* Calculate borders at beginning, check - */ - res = (*fct)(&fa, a, ps); - if ((0 == res) || (!(dk4ma_is_finite(fa)))) { - back = DK4_ITER_RESULT_E_FCT; $? "! function calculation for a" - goto finished; - } - res = (*fct)(&fb, b, ps); - if ((0 == res) || (!(dk4ma_is_finite(fb)))) { - back = DK4_ITER_RESULT_E_FCT; $? "! function calculation for b" - goto finished; - } - cond = find_condition(fa, fb); - if (CONDITION_ILLEGAL == cond) { $? "! invalid interval" - 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++; } $? ". pass %lu", passno - /* Calculate x position - */ - $? ". a = %lg fa = %lg b = %lg fb = %lg", a, fa, b, fb - 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 { $? ". x = %lg", x - /* Calculate function for x, if x is usable - */ - res = (*fct)(&fx, x, ps); - if (0 == res) { $? "! function calculation for x" - back = DK4_ITER_RESULT_E_FCT; - cc = -1; - } - if (0 == dk4ma_is_finite(fx)) { $? "! function calculation for x" - back = DK4_ITER_RESULT_E_INFINITE; - cc = -1; - } - afx = fabs(fx); - if (0 == dk4ma_is_finite(afx)) { $? "! absolute value for fx" - back = DK4_ITER_RESULT_E_INFINITE; - cc = -1; - } - /* Check whether we are done, if we could continue - */ - if (1 == cc) { $? ". y = %lg", fx - $!trace-code if (1UL < passno) { - $? ". x-delta = %lg", fabs(x - xo) - $!trace-code } - if (isgreater(mctx.eps_y,0.0)) { - /* Must check y */ - if (afx < mctx.eps_y) { $? ". y small enough" - if ((x == xo) && (1UL < passno)) { - cc = 0; $? ". OK exact match" - } - 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; $? ". OK small step" - } - } - } - else { - cc = 0; $? ". OK no x step restriction" - } - } - } - } - } - else { $? ". no y check required" - /* No y check */ - if ((x == xo) && (1UL < passno)) { - cc = 0; $? ". OK exact match" - } - 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; $? ". OK x step small enough" - } - } - } - 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: { $? ". use x as a" - 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; $? "! fb" - } - } 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; $? "! fb" - } - } 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; $? "! fb" - } - } break; - } - } - a = x; - fa = fx; - pc = CHANGED_A; - } break; - case CHANGED_B: { $? ". use x as b" - 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; $? "! fa" - } - } 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; $? "! fa" - } - } 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; $? "! fa" - } - } break; - } - } - b = x; - fb = fx; - pc = CHANGED_B; - } break; - default: { - /* ERROR: Must not happen */ - back = DK4_ITER_RESULT_E_CONV; $? "! bug" - cc = -1; - } break; - } - } - } - /* Check number of passes - */ - if ((1 == cc) && (0UL < mctx.maxpass) && (passno >= mctx.maxpass)) { - back = DK4_ITER_RESULT_E_PASSES; $? "! too many passes" - cc = -1; - } - } while (1 == cc); - - /* Success - */ - if (0 == cc) { - *d = x; - if (NULL != pp) { *pp = passno; } - back = DK4_ITER_RESULT_SUCCESS; - } - - finished: - $? "- dk4iter_interval %d", back - 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; - $? "+ dk4iter_newton" -#if DK4_USE_ASSERT - assert(NULL != d); - assert(NULL != fct); -#endif - res = (*fct)(v, x0, ps); - if (0 == res) { - cc = -1; $? "! function calculation" - back = DK4_ITER_RESULT_E_FCT; - } - if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) { - cc = -1; $? "! infinite values" - back = DK4_ITER_RESULT_E_INFINITE; - } - while (1 == cc) { - if (ULONG_MAX > passno) { passno++; } $? ". begin pass %lu", passno - xn = x0 - v[0] / v[1]; - if (dk4ma_is_finite(xn)) { $? ". x = %lg", xn - if (0 != (MINMAX_MINIMUM & (ctx->minmax))) { - if (0 != (MINMAX_MIN_EXCL & (ctx->minmax))) { - if (xn <= ctx->xmin) { - cc = -1; $? "! out of range, too small" - back = DK4_ITER_RESULT_E_OOR; - } - } - else { - if (xn < ctx->xmin) { - cc = -1; $? "! out of range, too small" - 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; $? "! out of range, too large" - back = DK4_ITER_RESULT_E_OOR; - } - } - else { - if (xn > ctx->xmax) { - cc = -1; $? "! out of range, too large" - back = DK4_ITER_RESULT_E_OOR; - } - } - } - if (1 == cc) { - res = (*fct)(v, xn, ps); - if (0 == res) { - cc = -1; $? "! function calculation" - back = DK4_ITER_RESULT_E_FCT; - } - if (!((dk4ma_is_finite(v[0])) && (dk4ma_is_finite(v[1])))) { - cc = -1; $? "! infinite values" - back = DK4_ITER_RESULT_E_INFINITE; - } - if (1 == cc) { - $? ". y = %lg", v[0] - $? ". dy/dx = %lg", v[1] - if (isgreater(ctx->eps_y,0.0)) { $? ". check y" - if (fabs(v[0]) < ctx->eps_y) { $? ". y in range" - if (xn == x0) { - cc = 0; $? ". OK exact value" - } - else { - if (0 == ctx->exact) { $? ". check step" - if(isgreater(ctx->eps_x,0.0)) { - if (isless(fabs(xn - x0),ctx->eps_x)) { - cc = 0; $? ". OK small step" - } - } - else { - cc = 0; $? ". OK no x step constraint" - } - } - } - } - } - else { $? ". no y check" - if (xn == x0) { - cc = 0; $? ". OK exact match" - } - else { - if (0 == ctx->exact) { - if (isgreater(ctx->eps_x,0.0)) { - if (isless(fabs(xn - x0),ctx->eps_x)) { - cc = 0; $? ". OK x step small enough" - } - } - else { - /* Neither eps_x nor eps_y, - so we wait for x=x0. - */ - } - } - } - } - } - x0 = xn; - } - } - else { - cc = -1; $? "! infinite value" - 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; - } - $? "- dk4iter_newton %d", back - 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; - $? "+ dk4iter_fix_point" -#if DK4_USE_ASSERT - assert(NULL != d); - assert(NULL != fct); -#endif - while (1 == cc) { - if (ULONG_MAX > passno) { passno++; } $? ". begin pass %lu", passno - /* - Calculate new x, check result - */ - res = (*fct)(&xn, x0, ps); - if (0 == res) { - cc = -1; $? "! function evaluation" - back = DK4_ITER_RESULT_E_FCT; - } - if (!(dk4ma_is_finite(xn))) { - cc = -1; $? "! infinite result" - back = DK4_ITER_RESULT_E_INFINITE; - } - if (1 == cc) { $? ". x = %lg", xn - /* - 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; $? "! out of range, too small" - back = DK4_ITER_RESULT_E_OOR; - } - } - else { - if (xn < ctx->xmin) { - cc = -1; $? "! out of range, too small" - 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; $? "! out of range, too large" - back = DK4_ITER_RESULT_E_OOR; - } - } - else { - if (xn > ctx->xmax) { - cc = -1; $? "! out of range, too large" - back = DK4_ITER_RESULT_E_OOR; - } - } - } - if (1 == cc) { - if (1 == cc) { - /* - Check whether we are finished - */ - if (xn == x0) { - cc = 0; $? ". OK exact match" - } - else { - if ((0 == ctx->exact) && (0.0 < ctx->eps_x)) { - if (isless(fabs(xn-x0),ctx->eps_x)) { - cc = 0; $? ". OK change small enough" - } -#if TRACE_DEBUG - else { $? ". change %lg", fabs(xn-x0) - } -#endif - } - else { $? "! No condition to check %lg", (xn-x0) - } - } - } - } - 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; - } - $? "- dk4iter_fix_point %d", back - 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 : */ - |