00001 // This may look like C code, but it is really -*- C++ -*- 00002 /* 00003 ************************************************************************ 00004 * 00005 * Numerical Math Package 00006 * 00007 * The present package implements various algorithms of Numerical Math 00008 * 00009 * $Id: math_num.h,v 1.1 2004/05/21 21:02:52 maxx Exp $ 00010 * 00011 ************************************************************************ 00012 */ 00013 00014 #ifndef __GNUC__ 00015 #pragma once 00016 #endif 00017 #ifndef _math_num_h 00018 #define _math_num_h 1 00019 00020 #if defined(__GNUC__) 00021 #pragma interface 00022 #endif 00023 00024 #include "myenv.h" 00025 #include <float.h> 00026 #include <math.h> 00027 #include "builtin.h" 00028 #include "std.h" 00029 00030 namespace linalg 00031 { 00032 00033 /* 00034 *------------------------------------------------------------------------ 00035 * Some constants 00036 * Compile and run the program epsilon.c to determine the values below for 00037 * your computer 00038 */ 00039 00040 //#define EPSILON 2.22045e-16 // DBL_EPSILON 00041 //#define SQRT_EPSILON 1.49012e-08 00042 00043 /* 00044 *------------------------------------------------------------------------ 00045 * Brent's minimum and zero finders for 00046 * a function of a single argument 00047 */ 00048 00049 // A function of a single argument 00050 struct UnivariateFunctor 00051 { 00052 virtual double operator() (const double x) = 0; 00053 }; 00054 00055 // Obtain a zero of function f 00056 // over the interval [ax,bx] with the 00057 // accuracy tol. 00058 double zeroin(const double ax, const double bx, 00059 UnivariateFunctor& f, const double tol=DBL_EPSILON); 00060 00061 // Find a minimum of function f 00062 // over the interval [a,b] with the 00063 // accuracy tol. 00064 // Returns an approx. to the min location 00065 double fminbr(const double a, const double b, 00066 UnivariateFunctor& f, const double tol=DBL_EPSILON); 00067 00068 class Vector; // Opaque class used by the routines below 00069 00070 /* 00071 *------------------------------------------------------------------------ 00072 * Interpolation of the function 00073 * specified in the tabular form 00074 */ 00075 00076 00077 // Aitken-Lagrange interpolation to the 00078 // point q over the table of function values 00079 // y[i] = y(x[i]), i = y.lwb..y.upb 00080 00081 // Uniform mesh x[i] = x0 + s*(i-y.lwb) 00082 double ali(const double q, const double x0, const double s, const Vector& y); 00083 // Nonuniform grid with nodes in x[i] 00084 double ali(const double q, const Vector& x, const Vector& y); 00085 00086 /* 00087 *------------------------------------------------------------------------ 00088 * Multi-dimensional minimization 00089 */ 00090 00091 // A function of a _vector_ argument 00092 struct MultivariateFunctor 00093 { 00094 virtual double operator() (const Vector& x) = 0; 00095 }; 00096 00097 // Find a local minimum of a given 00098 // function by the Hooke-Jeeves method 00099 double hjmin( // Return the function value at min 00100 Vector& b, // Input: initial guess to min loc 00101 // Output: loc for the min found 00102 Vector& h, // Input: initial values for the 00103 // steps along each b(i) 00104 // Output: final steps right before 00105 // the termination 00106 MultivariateFunctor& f // A function being optimized 00107 ); 00108 00109 00110 // The same as above with the only difference 00111 // initial steps are given to be the same 00112 // along every direction. The final steps 00113 // aren't reported back though 00114 double hjmin(Vector& b, const double h0, 00115 MultivariateFunctor& f); 00116 00117 #endif 00118 }