diff options
Diffstat (limited to 'src/root_base.c')
| -rw-r--r-- | src/root_base.c | 134 |
1 files changed, 134 insertions, 0 deletions
diff --git a/src/root_base.c b/src/root_base.c new file mode 100644 index 0000000..b13d9db --- /dev/null +++ b/src/root_base.c @@ -0,0 +1,134 @@ +/* +Copyright 2016 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* + Base functions for root finding + + see integral_*.h for the values taken by ROOT_FUNC, etc... +*/ + + +// compute the root of a real function of 1 variable using the Newton method +int ROOT_FUNC(root_newton_inplace) (ROOT_FLOAT_TYPE* out, int (*func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), int (*deriv_func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), ROOT_FLOAT_TYPE tolerance, unsigned int maxiter, void* extra_args){ + unsigned int count=0; + ROOT_FLOAT_TYPE valf, valdf, delta, tmp; + int ret; + + #ifdef ROOT_FLOAT_INIT + ROOT_FLOAT_INIT(valf); + ROOT_FLOAT_INIT(valdf); + ROOT_FLOAT_INIT(delta); + ROOT_FLOAT_INIT(tmp); + #endif + + // evaluate at guess + ret=(*func)(&valf, *out, extra_args); + if(ret<0){ + #ifdef ROOT_FLOAT_FREE + ROOT_FLOAT_FREE(valf); + ROOT_FLOAT_FREE(valdf); + ROOT_FLOAT_FREE(delta); + ROOT_FLOAT_FREE(tmp); + #endif + return(ret); + } + // check that valf is a number + if(! ROOT_FLOAT_ISNUMBER(valf)){ + #ifdef ROOT_FLOAT_FREE + ROOT_FLOAT_FREE(valf); + ROOT_FLOAT_FREE(valdf); + ROOT_FLOAT_FREE(delta); + ROOT_FLOAT_FREE(tmp); + #endif + return(LIBINUM_ERROR_NAN); + } + + ROOT_FLOAT_ABS(delta, valf); + + // loop until tolerance is reached + while(ROOT_FLOAT_CMP(delta, tolerance) > 0 && count < maxiter){ + // new guess + (*deriv_func)(&valdf, *out, extra_args); + if(ret<0){ + #ifdef ROOT_FLOAT_FREE + ROOT_FLOAT_FREE(valf); + ROOT_FLOAT_FREE(valdf); + ROOT_FLOAT_FREE(delta); + ROOT_FLOAT_FREE(tmp); + #endif + return(ret); + } + // check that valdf is a number that is not 0 + if(! ROOT_FLOAT_ISNUMBER(valdf) || ROOT_FLOAT_ISZERO(valdf)){ + #ifdef ROOT_FLOAT_FREE + ROOT_FLOAT_FREE(valf); + ROOT_FLOAT_FREE(valdf); + ROOT_FLOAT_FREE(delta); + ROOT_FLOAT_FREE(tmp); + #endif + return(LIBINUM_ERROR_NAN); + } + + ROOT_FLOAT_DIV(tmp, valf, valdf); + ROOT_FLOAT_SUB(*out, *out, tmp); + + // evaluate at new guess + (*func)(&valf, *out, extra_args); + if(ret<0){ + #ifdef ROOT_FLOAT_FREE + ROOT_FLOAT_FREE(valf); + ROOT_FLOAT_FREE(valdf); + ROOT_FLOAT_FREE(delta); + ROOT_FLOAT_FREE(tmp); + #endif + return(ret); + } + // check that valf is a number + if(! ROOT_FLOAT_ISNUMBER(valf)){ + #ifdef ROOT_FLOAT_FREE + ROOT_FLOAT_FREE(valf); + ROOT_FLOAT_FREE(valdf); + ROOT_FLOAT_FREE(delta); + ROOT_FLOAT_FREE(tmp); + #endif + return(LIBINUM_ERROR_NAN); + } + + ROOT_FLOAT_ABS(delta, valf); + + // increase counter + count++; + } + + // free variables + #ifdef ROOT_FLOAT_FREE + ROOT_FLOAT_FREE(valf); + ROOT_FLOAT_FREE(valdf); + ROOT_FLOAT_FREE(delta); + ROOT_FLOAT_FREE(tmp); + #endif + + // fail if maxiter was reached + if(count==maxiter){ + return(LIBINUM_ERROR_MAXITER); + } + return(0); +} +int ROOT_FUNC(root_newton) (ROOT_FLOAT_TYPE* out, int (*func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), int (*deriv_func)(ROOT_FLOAT_TYPE*, ROOT_FLOAT_TYPE, void*), ROOT_FLOAT_TYPE init, ROOT_FLOAT_TYPE tolerance, unsigned int maxiter, void* extra_args){ + ROOT_FLOAT_SET(*out, init); + ROOT_FUNC(root_newton_inplace) (out, func, deriv_func, tolerance, maxiter, extra_args); + return(0); +} |