%% 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
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
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
#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"
$!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 : */