fmprb_calc.h – calculus with real-valued functions

This module provides functions for operations of calculus over the real numbers (intended to include root-finding, optimization, integration, and so on). It is planned that the module will include two types of algorithms:

• Interval algorithms that give provably correct results. An example would be numerical integration on an interval by dividing the interval into small balls and evaluating the function on each ball, giving rigorous upper and lower bounds.
• Conventional numerical algorithms that use heuristics to estimate the accuracy of a result, without guaranteeing that it is correct. An example would be numerical integration based on pointwise evaluation, where the error is estimated by comparing the results with two different sets of evaluation points. Ball arithmetic then still tracks the accuracy of the function evaluations.

Any algorithms of the second kind will be clearly marked as such.

Types, macros and constants

fmprb_calc_func_t

Typedef for a pointer to a function with signature:

int func(fmprb_ptr out, const fmprb_t inp, void * param, long order, long prec)

implementing a univariate real function \(f(x)\). When called, func should write to out the first order coefficients in the Taylor series expansion of \(f(x)\) at the point inp, evaluated at a precision of prec bits. The param argument may be used to pass through additional parameters to the function. The return value is reserved for future use as an error code. It can be assumed that out and inp are not aliased and that order is positive.

FMPRB_CALC_SUCCESS

Return value indicating that an operation is successful.

FMPRB_CALC_IMPRECISE_INPUT

Return value indicating that the input to a function probably needs to be computed more accurately.

FMPRB_CALC_NO_CONVERGENCE

Return value indicating that an algorithm has failed to convergence, possibly due to the problem not having a solution, the algorithm not being applicable, or the precision being insufficient

Debugging

int fmprb_calc_verbose

If set, enables printing information about the calculation to standard output.

Root-finding

long fmprb_calc_isolate_roots(fmprb_ptr * found, int ** flags, fmprb_calc_func_t func, void * param, const fmprb_t interval, long maxdepth, long maxeval, long maxfound, long prec)

Rigorously isolates single roots of a real analytic function on the interior of an interval.

This routine writes an array of n interesting subintervals of interval to found and corresponding flags to flags, returning the integer n. The output has the following properties:

• The function has no roots on interval outside of the output subintervals.
• Subintervals are sorted in increasing order (with no overlap except possibly starting and ending with the same point).
• Subintervals with a flag of 1 contain exactly one (single) root.
• Subintervals with any other flag may or may not contain roots.

If no flags other than 1 occur, all roots of the function on interval have been isolated. If there are output subintervals on which the existence or nonexistence of roots could not be determined, the user may attempt further searches on those subintervals (possibly with increased precision and/or increased bounds for the breaking criteria). Note that roots of multiplicity higher than one and roots located exactly at endpoints cannot be isolated by the algorithm.

The following breaking criteria are implemented:

• At most maxdepth recursive subdivisions are attempted. The smallest details that can be distinguished are therefore about \(2^{-\text{maxdepth}}\) times the width of interval. A typical, reasonable value might be between 20 and 50.
• If the total number of tested subintervals exceeds maxeval, the algorithm is terminated and any untested subintervals are added to the output. The total number of calls to func is thereby restricted to a small multiple of maxeval (the actual count can be slightly higher depending on implementation details). A typical, reasonable value might be between 100 and 100000.
• The algorithm terminates if maxfound roots have been isolated. In particular, setting maxfound to 1 can be used to locate just one root of the function even if there are numerous roots. To try to find all roots, LONG_MAX may be passed.

The argument prec denotes the precision used to evaluate the function. It is possibly also used for some other arithmetic operations performed internally by the algorithm. Note that it probably does not make sense for maxdepth to exceed prec.

Warning: it is assumed that subdivision points of interval can be represented exactly as floating-point numbers in memory. Do not pass \(1 \pm 2^{-10^{100}}\) as input.

int fmprb_calc_refine_root_bisect(fmprb_t r, fmprb_calc_func_t func, void * param, const fmprb_t start, long iter, long prec)

Given an interval start known to contain a single root of func, refines it using iter bisection steps. The algorithm can return a failure code if the sign of the function at an evaluation point is ambiguous. The output r is set to a valid isolating interval (possibly just start) even if the algorithm fails.

void fmprb_calc_newton_conv_factor(fmpr_t conv_factor, fmprb_calc_func_t func, void * param, const fmprb_t conv_region, long prec)

Given an interval \(I\) specified by conv_region, evaluates a bound for \(C = \sup_{t,u \in I} \frac{1}{2} |f''(t)| / |f'(u)|\), where \(f\) is the function specified by func and param. The bound is obtained by evaluating \(f'(I)\) and \(f''(I)\) directly. If \(f\) is ill-conditioned, \(I\) may need to be extremely precise in order to get an effective, finite bound for C.

int fmprb_calc_newton_step(fmprb_t xnew, fmprb_calc_func_t func, void * param, const fmprb_t x, const fmprb_t conv_region, const fmpr_t conv_factor, long prec)

Performs a single step with an interval version of Newton’s method. The input consists of the function \(f\) specified by func and param, a ball \(x = [m-r, m+r]\) known to contain a single root of \(f\), a ball \(I\) (conv_region) containing \(x\) with an associated bound (conv_factor) for \(C = \sup_{t,u \in I} \frac{1}{2} |f''(t)| / |f'(u)|\), and a working precision prec.

The Newton update consists of setting \(x' = [m'-r', m'+r']\) where \(m' = m - f(m) / f'(m)\) and \(r' = C r^2\). The expression \(m - f(m) / f'(m)\) is evaluated using ball arithmetic at a working precision of prec bits, and the rounding error during this evaluation is accounted for in the output. We now check that \(x' \in I\) and \(r' < r\). If both conditions are satisfied, we set xnew to \(x'\) and return FMPRB_CALC_SUCCESS. If either condition fails, we set xnew to \(x\) and return FMPRB_CALC_NO_CONVERGENCE, indicating that no progress is made.

int fmprb_calc_refine_root_newton(fmprb_t r, fmprb_calc_func_t func, void * param, const fmprb_t start, const fmprb_t conv_region, const fmpr_t conv_factor, long eval_extra_prec, long prec)

Refines a precise estimate of a single root of a function to high precision by performing several Newton steps, using nearly optimally chosen doubling precision steps.

The inputs are defined as for fmprb_calc_newton_step, except for the precision parameters: prec is the target accuracy and eval_extra_prec is the estimated number of guard bits that need to be added to evaluate the function accurately close to the root (for example, if the function is a polynomial with large coefficients of alternating signs and Horner’s rule is used to evaluate it, the extra precision should typically be approximately the bit size of the coefficients).

This function returns FMPRB_CALC_SUCCESS if all attempted Newton steps are successful (note that this does not guarantee that the computed root is accurate to prec bits, which has to be verified by the user), only that it is more accurate than the starting ball.

On failure, FMPRB_CALC_IMPRECISE_INPUT or FMPRB_CALC_NO_CONVERGENCE may be returned. In this case, r is set to a ball for the root which is valid but likely does have full accuracy (it can possibly just be equal to the starting ball).