summaryrefslogtreecommitdiff
path: root/support/dktools/dk4iter.ctr
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2020-10-12 03:04:00 +0000
committerNorbert Preining <norbert@preining.info>2020-10-12 03:04:00 +0000
commit0ce40abb18ec02ec6fd6bcc5f21612c88daa7578 (patch)
tree416289fe1448873fd8ca33051f50ad85bffa8aaa /support/dktools/dk4iter.ctr
parentfdb18507cd80dc17f5a5256153d34668b4f4e61c (diff)
CTAN sync 202010120303
Diffstat (limited to 'support/dktools/dk4iter.ctr')
-rw-r--r--support/dktools/dk4iter.ctr1456
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 : */
-