135 lines
3.6 KiB
C
135 lines
3.6 KiB
C
/*
|
|
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);
|
|
}
|