00001 /* EGlib "Efficient General Library" provides some basic structures and 00002 * algorithms commons in many optimization algorithms. 00003 * 00004 * Copyright (C) 2005 Daniel Espinoza and Marcos Goycoolea. 00005 * 00006 * This library is free software; you can redistribute it and/or modify it 00007 * under the terms of the GNU Lesser General Public License as published by the 00008 * Free Software Foundation; either version 2.1 of the License, or (at your 00009 * option) any later version. 00010 * 00011 * This library is distributed in the hope that it will be useful, but 00012 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 00013 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 00014 * License for more details. 00015 * 00016 * You should have received a copy of the GNU Lesser General Public License 00017 * along with this library; if not, write to the Free Software Foundation, 00018 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 00019 * */ 00020 #ifndef __EG_LPNUM_H__ 00021 #define __EG_LPNUM_H__ 00022 /* ========================================================================= */ 00023 /** @defgroup EGlpNum EGlpNum 00024 * 00025 * Here we define a common interface to handle numbers in general, the idea is 00026 * to be able to work with infinite precicion numbers, plain doubles, floats, 00027 * integers, or fixed point numbers, without actually making different codes 00028 * for each of those types, rather we preffer to fix that at compyle time. 00029 * 00030 * @par History: 00031 * Revision 0.0.2 00032 * We start doing the migration to gmp and 'natural' types, this means that we 00033 * drop support for EGrat, this allow us to drop the requirement to use 00034 * pointers, and instead we can just call the functions with the original 00035 * parameters, still we have to be VERY carefull regarding changing 00036 * local/external copies. 00037 * - 2007-10-08 00038 * - Move EGswap, EGabs, EGmin and Egmax to eg_numutil.h 00039 * - 2005-10-31 00040 * - Add EGswap to swap elements of any predefined type. 00041 * - 2005-08-31 00042 * - Add EGmin and EGmax for built in types (i.e. for types where 00043 * the < comparison works as we want). 00044 * - 2005-08-16 00045 * - Streamline mpq_EGlpNumGetStr 00046 * - Minor Fixes for zeroLpNum 00047 * - 2005-07-29 00048 * - Add EGabs definition. 00049 * - 2005-07-24 00050 * - Split eg_lpnum.h into different headers for each type of 00051 * suported numbers. 00052 * - Deprecate EGlpNumCOmpUFrac 00053 * - 2005-05-26 00054 * - Add epsLpNum 00055 * - 2005-05-17 00056 * - Add mpq_EGlpNumReadStrXc(mpq_t,str) 00057 * - 2005-05-16 00058 * - Add mpq_EGlpNumSet_mpf(mpq,mpf) 00059 * - 2005-05-12 00060 * - Add mpf_EGlpNumEpow(num,power) 00061 * - Add function to change precision of the numbers on the fly. 00062 * - Add EGlpNumReadStr to set a number from a given input string. 00063 * - Add EGlpNumGetStr to get (hopefully) the exact representation 00064 * of the given input as string. 00065 * - 2005-05-03 00066 * - Change the structure of the header so it provides an interface 00067 * in which a program can use all types of numbers at the same time, 00068 * this implies that we must define a start-up and clean-up function 00069 * that would initialize all constants for all numbers, and change 00070 * the naming scheme accordingly to support this. 00071 * - Deprecate EGlpNumInitArray 00072 * - Deprecate EGlpNumFreeIntArray 00073 * - Change all static-inline definitions to Statement Exprs style. 00074 * - 2005-04-25 00075 * - Add EGlpNumIsSumLess(a,b,c) 00076 * - Add EGlpNumIsDiffLess(a,b,c) 00077 * - 2005-04-13 00078 * - Add EGlpNumCopyDiffRatio(a,b,c,d) 00079 * - Add EGlpNumIsEqqual(a,b) 00080 * - Add EGlpNumIsEqual(a,b,error) 00081 * - Add EGlpNumCopy(a,b) 00082 * - Add EGlpNumIsLess(a,b) 00083 * - Add EGlpNumToLf(a) 00084 * - Add EGlpNumCopyDiff(a,b,c) 00085 * - Add EGlpNumCopyAbs(a,b) 00086 * - Add EGlpNumSubTo(a,b) 00087 * - Add EGlpNumAddTo(a,b) 00088 * - Add EGlpNumDivTo(a,b) 00089 * - Add EGlpNumMultTo(a,b) 00090 * - Add EGlpNumZero(a) 00091 * - Add EGlpNumOne(a) 00092 * - Add EGlpNumAddInnProdTo(a,b,c) 00093 * - Add EGlpNumSubInnProdTo(a,b,c) 00094 * - Add EGlpNumSign(a) 00095 * - Add EGlpNumCopyNeg(a,b) 00096 * - Add EGlpNumDivUiTo(a,b) 00097 * - Add EGlpNumMultUiTo(a,b) 00098 * - Add EGlpNumIsLeq(a,b) 00099 * - Add EGlpNumCopySqrOver(a,b,c) 00100 * - Add EGlpNumSet(a,b) 00101 * - Add EGlpNumCopyFrac(a,b,c) 00102 * - Add EGlpNumAddUiTo(a,b) 00103 * - Add EGlpNumSubUiTo(a,b) 00104 * - Add EGlpNumCopySum(a,b,c) 00105 * - Add EGlpNumInv(a) 00106 * - Add EGlpNumFloor(a) 00107 * - Add EGlpNumCeil(a) 00108 * - Add EGlpNumIsLessDbl(a,b) 00109 * - Add EGlpNumIsGreaDbl(a,b) 00110 * - Add EGlpNumSetToMaxAbs(a,b) 00111 * - Add EGlpNumAllocArray(size) 00112 * - Add EGlpNumFreeArray(a) 00113 * - Add EGlpNumReallocArray(a,size) 00114 * - Add EGlpNumInitVar(a) 00115 * - Add EGlpNumClearVar(a) 00116 * - Add EGlpNumIsNeqq(a,b) 00117 * - Add EGlpNumIsNeq(a,b,c) 00118 * - Add EGlpNumIsNeqqZero(a) 00119 * - Add EGlpNumIsNeqZero(a,b) 00120 * - Add EGlpNumSetToMinAbs(a,b) 00121 * Revision 0.0.1 00122 * - 2004-07-15 00123 * - Add support for GNU_MP_F types 00124 * - 2004-07-14 00125 * - Add support for GNU_MP_Q types 00126 * - 2004-07-12 00127 * - Add support for EG-rationals 00128 * - 2004-06-21 00129 * - First Implementation/Definition 00130 * */ 00131 00132 /** @file 00133 * @ingroup EGlpNum */ 00134 /** @addtogroup EGlpNum */ 00135 /** @{ */ 00136 /* ========================================================================= */ 00137 #include <stdlib.h> 00138 #include <stdio.h> 00139 #include <string.h> 00140 #include <limits.h> 00141 #include <math.h> 00142 #include <float.h> 00143 /* ========================================================================= */ 00144 /** @name Number Types Definitions: 00145 * Define (as its name suggest) an internal identifier for the given 00146 * type. this definitions are provided to select different types of data at 00147 * compile time, thus allowing us to provide limited template support. */ 00148 /*@{*/ 00149 /** C double type. */ 00150 #define DBL_TYPE 0 00151 /** C float type. */ 00152 #define FLT_TYPE 1 00153 /** C int type. */ 00154 #define INT_TYPE 2 00155 /** EGlib #EGfp10_t type, this is an implementation of fixed precision 00156 * arithmetic with 10 bits for fractional representation. */ 00157 #define FP10_TYPE 3 00158 /** EGlib #EGfp20_t type, this is an implementation of fixed precision 00159 * arithmetic with 20 bits for fractional representation. */ 00160 #define FP20_TYPE 4 00161 /** EGlib #EGfp28_t type, this is an implementation of fixed precision 00162 * arithmetic with 28 bits for fractional representation. */ 00163 #define FP28_TYPE 5 00164 /** EGlib #EGfp25_t type, this is an implementation of fixed precision 00165 * arithmetic with 25 bits for fractional representation. */ 00166 #define FP25_TYPE 6 00167 #if HAVE_GNU_MP 00168 /** GNU_MP library mpz_t type */ 00169 #define GNU_MP_Z 8 00170 /** GNU_MP library mpq_t type */ 00171 #define GNU_MP_Q 9 00172 /** GNU_MP library mpf_t type */ 00173 #define GNU_MP_F 10 00174 #endif 00175 /** C long double type */ 00176 #define LDBL_TYPE 11 00177 /** C long long int type */ 00178 #define LLINT_TYPE 12 00179 /** SoftFloat 128-bit floating point numbner */ 00180 #define FLOAT128_TYPE 13 00181 /*@}*/ 00182 00183 00184 /* ========================================================================= */ 00185 #ifndef EGLPNUM_TYPE 00186 /** 00187 * @brief default type for EGLPNUM_TYPE 00188 * @par Description: 00189 * If the type of number for EGlpNum is not defined beforehand (eg. eg_config.h 00190 * or through make.conf) we default to GNU_MP_Q. 00191 * available options: 00192 * - GNU_MP_Q 00193 * - GNU_MP_F 00194 * - DBL_TYPE 00195 * - LDBL_TYPE 00196 * */ 00197 #define EGLPNUM_TYPE DBL_TYPE 00198 #endif 00199 #include "eg_lpnum.dbl.h" 00200 //#include "eg_lpnum.int.h" 00201 //#include "eg_lpnum.ldbl.h" 00202 //#include "eg_lpnum.llint.h" 00203 //#include "eg_lpnum.fp20.h" 00204 #if HAVE_SOFTFLOAT 00205 //#include "eg_lpnum.float128.h" 00206 //#include "softfloat.h" 00207 #endif 00208 #if HAVE_GNU_MP 00209 //#include "gmp.h" 00210 //#include "eg_lpnum.mpz.h" 00211 //#include "eg_lpnum.mpq.h" 00212 //#include "eg_lpnum.mpf.h" 00213 #endif 00214 #include "eg_macros.h" 00215 #include "eg_mem.h" 00216 00217 /* ========================================================================= */ 00218 /** @brief Debugging verbosity messages deped on the value of DEBUG (defined in 00219 * eg_configure.h) and on the value of EGLPNUM_DEBUGL macro defined here. 00220 * */ 00221 #define EGLPNUM_DEBUGL 100 00222 00223 #ifndef EGLPNUM_MINEPS 00224 /* ========================================================================= */ 00225 /** @brief This constant define the of the acuracy required while 00226 * converting doubles to rationals, a good number is 1e-5. More exactly, we 00227 * stop the continued fraction method whenever the next e_i-[e_i] computed is 00228 * less than EGLPNUM_MINEPS. Note that this value can't be smaller than 00229 * 1/ULONG_MAX, otherwise we will have problems in the confertion step. */ 00230 #define EGLPNUM_MINEPS 0x1ep-20 00231 #else 00232 #if EGLPNUM_MINEPS < 3e-10 00233 #undef EGLPNUM_MINEPS 00234 #define EGLPNUM_MINEPS 3e-10 00235 #endif 00236 #endif 00237 00238 #if HAVE_GNU_MP 00239 /* ========================================================================= */ 00240 /** @brief Set the default number of __BITS__ used in the precision of the 00241 * float point numbers (mpf_t), a normal double use up to 56-64 bits., the 00242 * default precision is set to 128 */ 00243 extern unsigned long int EGLPNUM_PRECISION; 00244 #endif 00245 00246 /* ========================================================================= */ 00247 /** @brief Change the default precision for mpf_t numbers. */ 00248 void EGlpNumSetPrecision (const unsigned prec); 00249 00250 /* ========================================================================= */ 00251 /** @brief Allocate an array of a given type and store (sizeof(size_t) bytes 00252 * before the actual array) the size of the allocated array. 00253 * @param type the type of the array to be returned. 00254 * @param size the length of the array to be returned, note that it can be 00255 * zero, in wich case no memory allocation is made and NULL is returned. */ 00256 #define __EGlpNumAllocArray(type,size) ({\ 00257 size_t __sz = (size);\ 00258 size_t *__utmp = __sz ? (size_t*) EGmalloc (sizeof(type) * __sz + sizeof(size_t)) : 0;\ 00259 if(__sz) __utmp[0] = __sz;\ 00260 (type*)(__sz ? (__utmp+1):0);}) 00261 00262 /* ========================================================================= */ 00263 /** @brief Given an array allocated with __EGlpNumAllocArray, return the size of 00264 * the given array, if the array is null, return zero. 00265 * @param array the array from where we need the size. */ 00266 /* ========================================================================= */ 00267 #define __EGlpNumArraySize(array) ({\ 00268 size_t *__utmp = (size_t*)(array);\ 00269 if(__utmp) __utmp--;\ 00270 __utmp ? __utmp[0]:0;}) 00271 00272 /* ========================================================================= */ 00273 /** @brief, given an array allocated by __EGlpNumAllocArray, free the allocated 00274 * memory. 00275 * @param array the array to be freed, it can be null. The given array will 00276 * always pooint to NULL when this function is done. 00277 * */ 00278 /* ========================================================================= */ 00279 #define __EGlpNumFreeArray(array) ({\ 00280 size_t *__utmp = (size_t*)(array);\ 00281 if(__utmp) free (__utmp-1);\ 00282 (array) = 0;}) 00283 00284 /** @}*/ 00285 /* ========================================================================= */ 00286 #endif
1.4.6