eg_lpnum.c

Go to the documentation of this file.
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 #include "eg_config.h"
00021 #include "eg_memslab.h"
00022 /** @file
00023  * @ingroup EGlpNum */
00024 /** @addtogroup EGlpNum */
00025 /** @{ */
00026 /* ========================================================================= */
00027 /** @brief This constant define the acuracy required while
00028  * converting doubles to rationals, a good number is 1e-5. More exactly, we 
00029  * stop the continued fraction method whenever the next e_i-[e_i] computed is
00030  * less than EGLPNUM_MINEPS. Note that this value can't be smaller than
00031  * 1/ULONG_MAX, otherwise we will have problems in the confertion step. */
00032 #ifndef EGLPNUM_MINEPS
00033 #define EGLPNUM_MINEPS 0x1ep-20
00034 #else
00035 #if EGLPNUM_MINEPS < 3e-10
00036 #undef EGLPNUM_MINEPS
00037 #define EGLPNUM_MINEPS 3e-10
00038 #endif
00039 #endif
00040 
00041 
00042 /* ========================================================================= */
00043 /** @brief if non-zero, use slab-pool allocator for GMP, otherwise, use malloc/
00044  * realloc / free, */
00045 #ifndef EG_LPNUM_MEMSLAB
00046 #define EG_LPNUM_MEMSLAB 1
00047 #endif
00048 /* ========================================================================= */
00049 /** @brief type-dependant constants and helper numbers @{ */
00050 #ifdef HAVE_LIBGMP
00051 mpz_t __zeroLpNum_mpz__;
00052 mpz_t __oneLpNum_mpz__;
00053 mpz_t __MaxLpNum_mpz__;
00054 mpz_t __MinLpNum_mpz__;
00055 mpq_t __zeroLpNum_mpq__;
00056 mpq_t __oneLpNum_mpq__;
00057 mpq_t __MaxLpNum_mpq__;
00058 mpq_t __MinLpNum_mpq__;
00059 mpf_t __zeroLpNum_mpf__;
00060 mpf_t __MaxLpNum_mpf__;
00061 mpf_t __MinLpNum_mpf__;
00062 mpf_t __oneLpNum_mpf__;
00063 mpf_t mpf_eps;
00064 unsigned long int EGLPNUM_PRECISION = 128;
00065 #endif
00066 #ifdef HAVE_SOFTFLOAT
00067 #if HAVE_SOFTFLOAT
00068 const float128 __zeroLpNum_float128__ = { 0, 0 };
00069 const float128 __oneLpNum_float128__ = {.high = 0x3fff000000000000LL,.low = 0 };
00070 const float128 float128_eps = {.high = 0x3f8f800000000000LL,.low = 0 };
00071 const float128 __MaxLpNum_float128__ = {.high = 0x43feffffffffffffLL,.low=0xf000000000000000LL};
00072 const float128 __MinLpNum_float128__ = {.high = 0xc3feffffffffffffLL,.low=0xf000000000000000LL};
00073 #endif
00074 #endif
00075 /** @} */
00076 
00077 #ifdef HAVE_LIBGMP
00078 /* ========================================================================= */
00079 static int __EGlpNum_setup=0;
00080 /* ========================================================================= */
00081 /** @name data to handle memory allocations within gmp */
00082 /** @{*/
00083 #define __GMP_MEM_VERBOSE 100
00084 #define __GMP_MEM_STATS__ 0
00085 #define __GMP_MEM_MAX__ 256
00086 static const uint8_t _EGgmpPlTable[257] = \
00087   { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,  /* pool for 0-16 bytes */
00088     1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,  /* pool for 17-32 bytes */
00089     2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,  /* pool for 33-64 bytes */
00090     2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
00091     3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,  /* pool for 65-128 bytes */
00092     3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
00093     3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
00094     3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
00095     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,  /* pool for 129-255 bytes */
00096     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
00097     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
00098     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
00099     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
00100     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
00101     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
00102     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
00103 static const size_t _EGgmpPlSz[5] = {16,32,64,128,256};
00104 #define __GMP_MEM_NPOOL__ 5
00105 EGmemSlabPool_t EGgmpPl[__GMP_MEM_NPOOL__];
00106 #if __GMP_MEM_STATS__
00107 static size_t __alloc_sz[__GMP_MEM_NPOOL__] = {0,0,0,0,0};
00108 static size_t __rlloc_sz[__GMP_MEM_NPOOL__] = {0,0,0,0,0};
00109 static size_t __totmem = 0;
00110 static size_t __maxmem = 0;
00111 static size_t __nallocs = 0;
00112 static size_t __nrllocs = 0;
00113 static size_t __nalarge = 0;
00114 static size_t __nrlarge = 0;
00115 static size_t __laaverage=0;
00116 static size_t __lraverage=0;
00117 #endif
00118 /** @}*/
00119 /* ========================================================================= */
00120 #if __GMP_MEM_STATS__
00121 #define _TACC(__m) do{\
00122   __totmem += __m;\
00123   if(__totmem > __maxmem) __maxmem = __totmem;}while(0)
00124 #else
00125 #define _TACC(__m)
00126 #endif
00127 /* ========================================================================= */
00128 #if __GMP_MEM_STATS__
00129 #define _AACC(__a) do{\
00130   __nallocs++;\
00131   if((__a) <= __GMP_MEM_MAX__)\
00132     __alloc_sz[_EGgmpPlTable[__a]]++;\
00133   else{\
00134     __nalarge++;\
00135     __laaverage += (__a);}}while(0)
00136 #else
00137 #define _AACC(__a)
00138 #endif
00139 /* ========================================================================= */
00140 /** @brief dummy malloc function for gmp, at this stage is used for 
00141  * creating an account of the allocation */
00142 static void* __EGgmp_malloc(size_t sz)
00143 {
00144   void*ptr = 0;
00145   _TACC(sz);
00146   _AACC(sz);
00147   if(sz <= __GMP_MEM_MAX__)
00148   {
00149     ptr = EGmemSlabPoolAlloc(EGgmpPl+_EGgmpPlTable[sz]);
00150     MESSAGE(__GMP_MEM_VERBOSE,"alloc %p [%zd]",ptr, sz);
00151     return ptr;
00152     /*memset(ptr,0,sz);*/
00153   }
00154   else
00155   {
00156     ptr = malloc(sz);
00157     if(!ptr) EXIT(1,"No more memory");
00158     MESSAGE(__GMP_MEM_VERBOSE,"alloc %p [%zd]",ptr, sz);
00159     return ptr;
00160   }
00161   /*fprintf(stderr,"allocating %p [%zd]\n",ptr,sz);*/
00162   /*return ptr;*/
00163 }
00164 /* ========================================================================= */
00165 #if __GMP_MEM_STATS__
00166 #define _RACC(__a) do{\
00167   __nrllocs++;\
00168   if(__a<__GMP_MEM_MAX__)\
00169     __rlloc_sz[_EGgmpPlTable[__a]]++;\
00170   else{\
00171     __nrlarge++;\
00172     __lraverage += (__a);}}while(0)
00173 #else
00174 #define _RACC(__a)
00175 #endif
00176 /* ========================================================================= */
00177 /** @brief dummy realloc for gmp */
00178 static void* __EGgmp_realloc(void*ptr,size_t osz,size_t nsz)
00179 {
00180   const size_t a1 = nsz > __GMP_MEM_MAX__ ? __GMP_MEM_NPOOL__ : _EGgmpPlTable[nsz];
00181   const size_t a2 = osz > __GMP_MEM_MAX__ ? __GMP_MEM_NPOOL__ : _EGgmpPlTable[osz];
00182   const size_t msz = nsz > osz ? osz : nsz;
00183   void*rptr=0;
00184   #if __GMP_MEM_STATS__
00185   _RACC(nsz);
00186   __totmem -= osz;
00187   _TACC(nsz);
00188   #endif
00189   if(a1 < __GMP_MEM_NPOOL__)
00190   {
00191     if(a2 == a1)
00192     {
00193       rptr = ptr;
00194       MESSAGE(__GMP_MEM_VERBOSE,"realloc %p [%zd] to %p [%zd]",ptr, osz, rptr, nsz);
00195       return rptr;
00196     }
00197     else
00198     {
00199       rptr = EGmemSlabPoolAlloc(EGgmpPl+a1);
00200       memcpy(rptr,ptr,msz);
00201       if(a2 < __GMP_MEM_NPOOL__)
00202       {
00203         MESSAGE(__GMP_MEM_VERBOSE,"realloc %p [%zd] to %p [%zd]",ptr, osz, rptr, nsz);
00204         EGmemSlabPoolFree(ptr);
00205         return rptr;
00206       }
00207       else
00208       {
00209         MESSAGE(__GMP_MEM_VERBOSE,"realloc %p [%zd] to %p [%zd]",ptr, osz, rptr, nsz);
00210         EGfree(ptr);
00211         return rptr;
00212       }
00213     }
00214   }
00215   else if(a2 < __GMP_MEM_NPOOL__)
00216   {
00217     rptr = EGmalloc(nsz);
00218     memcpy(rptr,ptr,msz);
00219     MESSAGE(__GMP_MEM_VERBOSE,"realloc %p [%zd] to %p [%zd]",ptr, osz, rptr, nsz);
00220     EGmemSlabPoolFree(ptr);
00221     return rptr;
00222   }
00223   else
00224   {
00225     rptr = EGrealloc(ptr,nsz);
00226     MESSAGE(__GMP_MEM_VERBOSE,"realloc %p [%zd] to %p [%zd]",ptr, osz, rptr, nsz);
00227     return rptr;
00228   }
00229   /*fprintf(stderr,"Re-allocating %p [%zd] to %p [%zd]\n",ptr,osz,rptr,nsz);*/
00230   /*return rptr;*/
00231 }
00232 /* ========================================================================= */
00233 /** @brief dummy free for gmp */
00234 static void __EGgmp_free(void*ptr,size_t sz)
00235 {
00236   const size_t a = sz > __GMP_MEM_MAX__ ? __GMP_MEM_NPOOL__ : _EGgmpPlTable[sz];
00237   #if __GMP_MEM_STATS__
00238   __totmem -= sz;
00239   #endif
00240   MESSAGE(__GMP_MEM_VERBOSE,"freeing %p",ptr);
00241   if(a<__GMP_MEM_NPOOL__)
00242   {
00243     EGmemSlabPoolFree(ptr);
00244     return;
00245   }
00246   else
00247   {
00248     EGfree(ptr);
00249     return;
00250   }
00251 }
00252 #endif
00253 
00254 /* ========================================================================= */
00255 /*void EGlpNumStart(void) __attribute__ ((constructor));*/
00256 void EGlpNumStart(void)
00257 {
00258   #ifdef HAVE_LIBGMP
00259   int rval=0;
00260   register int i;
00261   if(__EGlpNum_setup) return;
00262   if(EG_LPNUM_MEMSLAB)
00263   {
00264     fprintf(stderr,"Using EG-GMP mempool\n");
00265     for( i = __GMP_MEM_NPOOL__ ; i-- ; )
00266     {
00267       EGmemSlabPoolInit(EGgmpPl+i,_EGgmpPlSz[i],0,0);
00268       rval = EGmemSlabPoolSetParam(EGgmpPl+i,EG_MSLBP_FREEFREE,0);
00269       EXIT(rval,"Unknown error");
00270     }
00271     mp_set_memory_functions(__EGgmp_malloc, __EGgmp_realloc, __EGgmp_free);
00272   }
00273   else
00274     fprintf(stderr,"No EG-GMP mempool\n");
00275   mpf_set_default_prec (EGLPNUM_PRECISION);
00276   mpz_init (__zeroLpNum_mpz__);
00277   mpz_init (__oneLpNum_mpz__);
00278   mpz_init (__MaxLpNum_mpz__);
00279   mpz_init (__MinLpNum_mpz__);
00280   mpz_set_ui (__zeroLpNum_mpz__, (unsigned long int)0);
00281   mpz_set_ui (__oneLpNum_mpz__, (unsigned long int)1);
00282   mpq_init (__MaxLpNum_mpq__);
00283   mpq_init (__MinLpNum_mpq__);
00284   mpf_init (__MaxLpNum_mpf__);
00285   mpf_init (__MinLpNum_mpf__);
00286   mpf_init (__zeroLpNum_mpf__);
00287   mpf_init (__oneLpNum_mpf__);
00288   mpf_set_ui(__MaxLpNum_mpf__,1UL);
00289   mpf_set_si(__MinLpNum_mpf__,-1L);
00290   mpf_mul_2exp(__MaxLpNum_mpf__,__MaxLpNum_mpf__,4096);
00291   /*mpf_mul_2exp(__MaxLpNum_mpf__,__MaxLpNum_mpf__,ULONG_MAX);*/
00292   mpf_mul_2exp(__MinLpNum_mpf__,__MinLpNum_mpf__,4096);
00293   /*mpf_mul_2exp(__MinLpNum_mpf__,__MinLpNum_mpf__,ULONG_MAX);*/
00294   mpq_set_f(__MaxLpNum_mpq__,__MaxLpNum_mpf__);
00295   mpq_set_f(__MinLpNum_mpq__,__MinLpNum_mpf__);
00296   mpz_set_f(__MaxLpNum_mpz__,__MaxLpNum_mpf__);
00297   mpz_set_f(__MinLpNum_mpz__,__MinLpNum_mpf__);
00298   mpf_set_ui (__oneLpNum_mpf__, (unsigned long int)1);
00299   mpf_set_ui (__zeroLpNum_mpf__, (unsigned long int)0);
00300   mpf_init_set_ui (mpf_eps, (unsigned long int)1);
00301   mpf_div_2exp (mpf_eps, mpf_eps,  (unsigned long int)(EGLPNUM_PRECISION - 1));
00302   mpq_init (__zeroLpNum_mpq__);
00303   mpq_init (__oneLpNum_mpq__);
00304   mpq_set_ui (__oneLpNum_mpq__, (unsigned long int)1, (unsigned long int)1);
00305   mpq_set_ui (__zeroLpNum_mpq__, (unsigned long int)0, (unsigned long int)1);
00306   __EGlpNum_setup=1;
00307   #endif
00308 }
00309 
00310 /* ========================================================================= */
00311 #ifdef HAVE_LIBGMP
00312 void EGlpNumSetPrecision (const unsigned prec)
00313 {
00314   EGLPNUM_PRECISION = prec;
00315   mpf_set_default_prec (EGLPNUM_PRECISION);
00316   mpf_clear (mpf_eps);
00317   mpf_init_set_ui (mpf_eps, (unsigned long int)1);
00318   mpf_div_2exp (mpf_eps, mpf_eps,  (unsigned long int)(EGLPNUM_PRECISION - 1));
00319 }
00320 #endif
00321 
00322 /* ========================================================================= */
00323 /*void EGlpNumExit(void) __attribute__ ((destructor));*/
00324 void EGlpNumClear(void)
00325 {
00326   #ifdef HAVE_LIBGMP
00327   #if __GMP_MEM_STATS__
00328   const char mc[5][3] = {"b ","Kb","Mb","Gb","Tb"};
00329   #endif
00330   int i;
00331   if(!__EGlpNum_setup) return;
00332   mpf_clear (__zeroLpNum_mpf__);
00333   mpf_clear (__oneLpNum_mpf__);
00334   mpf_clear (__MaxLpNum_mpf__);
00335   mpf_clear (__MinLpNum_mpf__);
00336   mpf_clear (mpf_eps);
00337   mpq_clear (__zeroLpNum_mpq__);
00338   mpq_clear (__oneLpNum_mpq__);
00339   mpq_clear (__MinLpNum_mpq__);
00340   mpq_clear (__MaxLpNum_mpq__);
00341   mpz_clear (__zeroLpNum_mpz__);
00342   mpz_clear (__oneLpNum_mpz__);
00343   mpz_clear (__MaxLpNum_mpz__);
00344   mpz_clear (__MinLpNum_mpz__);
00345   if(EG_LPNUM_MEMSLAB)
00346   {
00347     mp_set_memory_functions(0, 0, 0);
00348     for(i = __GMP_MEM_NPOOL__ ; i-- ; )
00349     {
00350       EGmemSlabPoolClear(EGgmpPl+i);
00351     }
00352     #if __GMP_MEM_STATS__
00353     fprintf(stderr,"GMP alloc statistics:\n");
00354     for( i = 0 ; i < __GMP_MEM_NPOOL__ ; i++)
00355     {
00356       if(__maxmem > 1024*1024) __maxmem= (__maxmem+1023)/1024;
00357       else break;
00358     }
00359     fprintf(stderr,"\tmaximum memory allocated      : %8.3lf %s\n",
00360             ((double)__maxmem)/1024, mc[i+1]);
00361     fprintf(stderr,"\tmalloc calls                  : %11zd\n",__nallocs);
00362     fprintf(stderr,"\trealloc calls                 : %11zd\n",__nrllocs);
00363     fprintf(stderr,"\tsmall size allocs-reallocs:\n");
00364     for( i = 0 ; i < __GMP_MEM_NPOOL__ ; i++)
00365     {
00366       if(__alloc_sz[i] || __rlloc_sz[i])
00367         fprintf(stderr, "\t%4d %11zd (%5.2lf%%) %11zd (%5.2lf%%)\n", 
00368                 _EGgmpPlSz[i], __alloc_sz[i], 
00369                 100.0*((double)__alloc_sz[i])/(__nallocs+__nrllocs),
00370                 __rlloc_sz[i], 
00371                 100.0*((double)__rlloc_sz[i])/(__nrllocs+__nallocs));
00372     }
00373     if(__nalarge)
00374       fprintf(stderr, "\tlarge size allocs (cals/avg)  : %11zd %11zd\n",
00375               __nalarge,__laaverage/__nalarge);
00376     if(__nrlarge)
00377       fprintf(stderr, "\tlarge size reallocs (cals/avg): %11zd %11zd\n",
00378               __nrlarge,__lraverage/__nalarge);
00379     #endif
00380     fprintf(stderr,"Disabling EG-GMP mempool\n");
00381   }
00382   __EGlpNum_setup=0;
00383   #endif
00384 }
00385 
00386 #ifdef HAVE_LIBGMP
00387 /* ========================================================================= */
00388 void mpq_EGlpNumSet_mpf (mpq_t var,
00389                          mpf_t flt)
00390 {
00391   /* local variables */
00392   unsigned long int __lsgn = mpf_cmp_ui (flt, (unsigned long int)0) < 0 ? (unsigned long int)1 : (unsigned long int)0;
00393   mpz_t __utmp,
00394     __z[7],
00395     max_den;
00396   long int __lexp = 0;
00397   int i;
00398   unsigned long int uexp;
00399   mpf_t __cvl,__lpnum__;
00400   mpf_init(__lpnum__);
00401   /* check if the given number is zero, if so, set to zero var and return */
00402   if (mpf_cmp_ui (flt, (unsigned long int)0) == 0)
00403   {
00404     mpq_set_ui (var, (unsigned long int)0, (unsigned long int)1);
00405     return;
00406   }
00407   /* if not, then we have some work to do */
00408   /* now we initialize the internal numbers */
00409   mpf_init (__cvl);
00410   mpf_abs (__cvl, flt);
00411   mpz_init_set_ui (__utmp, (unsigned long int)0);
00412   for (i = 7; i--;)
00413     mpz_init_set_ui (__z[i], (unsigned long int)0);
00414   mpz_set_ui (__z[0], (unsigned long int)1);
00415   mpz_set_ui (__z[4], (unsigned long int)1);
00416   /* max_den is the maximum denominator that we want to see, this number should
00417    * be sligtly larger than the square root of 2^EGLPNUM_PRECISION */
00418   mpz_init_set_ui (max_den, (unsigned long int)1);
00419   mpz_mul_2exp (max_den, max_den, EGLPNUM_PRECISION >> 1);
00420   /* first we compute the exponent stored in the limbs */
00421   /* now we compute the 2^n part needed to set this number between 0.5 and 1 */
00422   mpf_get_d_2exp(&__lexp,__cvl);
00423   if (__lexp < 0)
00424   {
00425     uexp = (unsigned long int)(-__lexp);
00426     mpf_mul_2exp (__cvl, __cvl, (unsigned long int) uexp);
00427   }
00428   else
00429   {
00430     uexp = (unsigned long int)(__lexp);
00431     mpf_div_2exp (__cvl, __cvl, (unsigned long int) uexp);
00432   }
00433   /* now we loop until the next t's is more than mpf_eps */
00434   /* the formula is 
00435    * p_i = t_i*p_{i-1} + p_{i-2}, and 
00436    * q_i = t_i*q_{i-1} + q_{i-2} 
00437    * note that |x-p_i/q_i|<1/q_i^2
00438    * for us t_i = __utmp, and the current number is either [0,1,2] in the __z
00439    * array, we use those popsitions ciclicly, and use the four position as a
00440    * temporary number, __z+4 is used to store q's, at the beginning i = 1. */
00441   while (1)
00442   {
00443     if (mpf_cmp (__cvl, mpf_eps) < 0 || (mpz_cmp (__z[4], max_den) > 0))
00444     {
00445       mpz_set (mpq_denref (var), __z[4]);
00446       mpz_set (mpq_numref (var), __z[1]);
00447       break;
00448     }
00449     /* first run */
00450     mpf_ui_div (__cvl, (unsigned long int)1, __cvl);
00451     mpz_set_f (__utmp, __cvl);
00452     mpf_set_z (__lpnum__, __utmp);
00453     mpf_sub (__cvl, __cvl, __lpnum__);
00454     mpz_set (__z[6], __utmp);
00455     mpz_set (__z[2], __z[0]);
00456     mpz_addmul (__z[2], __z[1], __z[6]);
00457     mpz_set (__z[5], __z[3]);
00458     mpz_addmul (__z[5], __z[4], __z[6]);
00459     if (mpf_cmp (__cvl, mpf_eps) < 0 || (mpz_cmp (__z[5], max_den) > 0))
00460     {
00461       mpz_set (mpq_denref (var), __z[5]);
00462       mpz_set (mpq_numref (var), __z[2]);
00463       break;
00464     }
00465     /* second run */
00466     mpf_ui_div (__cvl, (unsigned long int)1, __cvl);
00467     mpz_set_f (__utmp, __cvl);
00468     mpf_set_z (__lpnum__, __utmp);
00469     mpf_sub (__cvl, __cvl, __lpnum__);
00470     mpz_set (__z[6], __utmp);
00471     mpz_set (__z[0], __z[1]);
00472     mpz_addmul (__z[0], __z[2], __z[6]);
00473     mpz_set (__z[3], __z[4]);
00474     mpz_addmul (__z[3], __z[5], __z[6]);
00475     if (mpf_cmp (__cvl, mpf_eps) < 0 || (mpz_cmp (__z[3], max_den) > 0))
00476     {
00477       mpz_set (mpq_denref (var), __z[3]);
00478       mpz_set (mpq_numref (var), __z[0]);
00479       break;
00480     }
00481     /* third run */
00482     mpf_ui_div (__cvl, (unsigned long int)1, __cvl);
00483     mpz_set_f (__utmp, __cvl);
00484     mpf_set_z (__lpnum__, __utmp);
00485     mpf_sub (__cvl, __cvl, __lpnum__);
00486     mpz_set (__z[6], __utmp);
00487     mpz_set (__z[1], __z[2]);
00488     mpz_addmul (__z[1], __z[0], __z[6]);
00489     mpz_set (__z[4], __z[5]);
00490     mpz_addmul (__z[4], __z[3], __z[6]);
00491   }
00492 
00493   /* ending */
00494   mpq_canonicalize (var);
00495   if (__lsgn)
00496     mpq_neg (var, var);
00497   if (__lexp > 0)
00498     mpq_mul_2exp (var, var, (unsigned long int) __lexp);
00499   if (__lexp < 0)
00500     mpq_div_2exp (var, var, (unsigned long int) (-__lexp));
00501   for (i = 7; i--;)
00502     mpz_clear (__z[i]);
00503   mpf_clear (__cvl);
00504   mpz_clear (max_den);
00505   mpz_clear (__utmp);
00506   mpf_clear(__lpnum__);
00507   return;
00508 }
00509 /* ========================================================================= */
00510 /** @brief verbosity of continued fraction conversion method */
00511 #ifndef MPQ_VERBOSE_CNT_FRAC
00512 #define MPQ_VERBOSE_CNT_FRAC 1000
00513 #endif
00514 /* ========================================================================= */
00515 void mpq_EGlpNumSet (mpq_t var,
00516                      double const dbl)
00517 {
00518   /* local variables */
00519   double __dbl = dbl;
00520   unsigned __lsgn = __dbl > 0.0 ? 0U: 1U;
00521   unsigned long __utmp = 0;
00522   int __lexp = 0;
00523   double __cvl = __dbl = fabs (__dbl);
00524   /* we use the first three numbers for p, and the last three numbers for q */
00525   /* first check that the dbl is not zero */
00526   if (__dbl < 1e-151)
00527   {
00528     mpq_set_ui (var, (unsigned long int)0, (unsigned long int)1);
00529     __lsgn = 0;
00530   }
00531   else if (__dbl > 1e151)
00532     mpq_set_d (var, __dbl);
00533   else
00534   {
00535     /* now we initialize the integer numbers */
00536     mpz_t __z[7];
00537     for (__utmp = 7; __utmp--;)
00538       mpz_init (__z[__utmp]);
00539     mpz_set_ui (__z[0], (unsigned long int)1);
00540     mpz_set_ui (__z[4], (unsigned long int)1);
00541     /* now we compute the 2^n part needed to set this number between 0 and 1 */
00542 #define __HI_EXP(x,e,v,lv) {if( x >=v ){ e = e + lv; x /= v;}}
00543     if (__cvl > 1)
00544     {
00545       __HI_EXP (__cvl, __lexp,
00546                 115792089237316195423570985008687907853269984665640564039457584007913129639936.0,
00547                 256);
00548       __HI_EXP (__cvl, __lexp, 340282366920938463463374607431768211456.0, 128);
00549       __HI_EXP (__cvl, __lexp, 18446744073709551616.0, 64);
00550       __HI_EXP (__cvl, __lexp, 4294967296.0, 32);
00551       __HI_EXP (__cvl, __lexp, 65536.0, 16);
00552       __HI_EXP (__cvl, __lexp, 256.0, 8);
00553       __HI_EXP (__cvl, __lexp, 16.0, 4);
00554       __HI_EXP (__cvl, __lexp, 4.0, 2);
00555       __HI_EXP (__cvl, __lexp, 2.0, 1);
00556 #undef __HI_EXP
00557     }
00558     else if (__cvl < 0.5)
00559     {
00560 #define __LO_EXP(x,e,v,lv) {if( x < 1/v ) { e = e - lv; x *= v;}}
00561       __LO_EXP (__cvl, __lexp,
00562                 115792089237316195423570985008687907853269984665640564039457584007913129639936.0,
00563                 256);
00564       __LO_EXP (__cvl, __lexp, 340282366920938463463374607431768211456.0, 128);
00565       __LO_EXP (__cvl, __lexp, 18446744073709551616.0, 64);
00566       __LO_EXP (__cvl, __lexp, 4294967296.0, 32);
00567       __LO_EXP (__cvl, __lexp, 65536.0, 16);
00568       __LO_EXP (__cvl, __lexp, 256.0, 8);
00569       __LO_EXP (__cvl, __lexp, 16.0, 4);
00570       __LO_EXP (__cvl, __lexp, 4.0, 2);
00571       __LO_EXP (__cvl, __lexp, 2.0, 1);
00572 #undef __LO_EXP
00573     }
00574     /* now we loop until the next t's is more than EGLPNUM_MINEPS */
00575     /* the formula is 
00576      * p_i = t_i*p_{i-1} + p_{i-2}, and 
00577      * q_i = t_i*q_{i-1} + q_{i-2} 
00578      * note that |x-p_i/q_i|<1/q_i^2
00579      * for us t_i = __utmp, and the current number is either [0,1,2] in the __z
00580      * array, we use those popsitions ciclicly, and use the four position as a
00581      * temporary number, __z+4 is used to store q's, at the beginning i = 1. */
00582     while (1)
00583     {
00584       if (__cvl < EGLPNUM_MINEPS || (mpz_cmp_ui (__z[4], (unsigned long int)0xfffffff) > 0))
00585       {
00586         mpz_set (mpq_denref (var), __z[4]);
00587         mpz_set (mpq_numref (var), __z[1]);
00588         break;
00589       }
00590       MESSAGE(MPQ_VERBOSE_CNT_FRAC,"cur approximate %ld/%ld, error %10.7lg", 
00591               mpz_get_si(__z[1]), mpz_get_si(__z[4]), __cvl);
00592       /* first run */
00593       __cvl = 1 / __cvl;
00594       __dbl = floor(__cvl);
00595       __utmp = (unsigned long)__dbl;
00596       __cvl -= __dbl;
00597       mpz_set_ui (__z[6], __utmp);
00598       mpz_set (__z[2], __z[0]);
00599       mpz_addmul (__z[2], __z[1], __z[6]);
00600       mpz_set (__z[5], __z[3]);
00601       mpz_addmul (__z[5], __z[4], __z[6]);
00602       if (__cvl < EGLPNUM_MINEPS || (mpz_cmp_ui (__z[5], (unsigned long int)0xfffffff) > 0))
00603       {
00604         mpz_set (mpq_denref (var), __z[5]);
00605         mpz_set (mpq_numref (var), __z[2]);
00606         break;
00607       }
00608       MESSAGE(MPQ_VERBOSE_CNT_FRAC,"cur approximate %ld/%ld, error %10.7lg", 
00609               mpz_get_si(__z[2]), mpz_get_si(__z[5]), __cvl);
00610       /* second run */
00611       __cvl = 1 / __cvl;
00612       __dbl = floor(__cvl);
00613       __utmp = (unsigned long)(__dbl);
00614       __cvl -= __dbl;
00615       mpz_set_ui (__z[6], __utmp);
00616       mpz_set (__z[0], __z[1]);
00617       mpz_addmul (__z[0], __z[2], __z[6]);
00618       mpz_set (__z[3], __z[4]);
00619       mpz_addmul (__z[3], __z[5], __z[6]);
00620       if (__cvl < EGLPNUM_MINEPS || (mpz_cmp_ui (__z[3], (unsigned long int)0xfffffff) > 0))
00621       {
00622         mpz_set (mpq_denref (var), __z[3]);
00623         mpz_set (mpq_numref (var), __z[0]);
00624         break;
00625       }
00626       MESSAGE(MPQ_VERBOSE_CNT_FRAC,"cur approximate %ld/%ld, error %10.7lg", 
00627               mpz_get_si(__z[0]), mpz_get_si(__z[3]), __cvl);
00628       /* third run */
00629       __cvl = 1 / __cvl;
00630       __dbl = floor(__cvl);
00631       __utmp = (unsigned long)(__dbl);
00632       __cvl -= __dbl;
00633       mpz_set_ui (__z[6], __utmp);
00634       mpz_set (__z[1], __z[2]);
00635       mpz_addmul (__z[1], __z[0], __z[6]);
00636       mpz_set (__z[4], __z[5]);
00637       mpz_addmul (__z[4], __z[3], __z[6]);
00638     }
00639     for (__utmp = 7; __utmp--;)
00640       mpz_clear (__z[__utmp]);
00641   }
00642   /* ending */
00643   mpq_canonicalize (var);
00644   if (__lsgn)
00645     mpq_neg (var, var);
00646   if (__lexp > 0)
00647     mpq_mul_2exp (var, var, (unsigned long int) __lexp);
00648   if (__lexp < 0)
00649     mpq_div_2exp (var, var, (unsigned long int) (-__lexp));
00650   return;
00651 }
00652 
00653 /* ========================================================================= */
00654 int mpz_EGlpNumReadStr (mpz_t var,
00655                         char const *str)
00656 {
00657   /* local variables */
00658   char unsigned a_sgn = 1;
00659   char unsigned sgn = 0;
00660   char c = 0;
00661   int n_char = 0;
00662   /* now we read the string */
00663   c = str[n_char];
00664   mpz_set_ui (var, (unsigned long int)0);
00665   while ((('0' <= c) && (c <= '9')) ||  /* allow to read digits */
00666          (a_sgn && (c == '+' || c == '-')) /* allow sign for exponent */ )
00667   {
00668     switch (c)
00669     {
00670     case '0':
00671     case '1':
00672     case '2':
00673     case '3':
00674     case '4':
00675     case '5':
00676     case '6':
00677     case '7':
00678     case '8':
00679     case '9':
00680       mpz_mul_ui (var, var, (unsigned long int)10);
00681       mpz_add_ui (var, var,  (unsigned long int)(c - '0'));
00682       a_sgn = 0;
00683       break;
00684     case '-':
00685       sgn = 1;
00686     case '+':
00687       a_sgn = 0;
00688       break;
00689     }
00690     /* advance the reading character */
00691     c = str[++n_char];
00692   }
00693   if (sgn)
00694     mpz_neg (var, var);
00695   return n_char;
00696 }
00697 
00698 /* ========================================================================= */
00699 int mpq_EGlpNumReadStrXc (mpq_t var,
00700                           char const *str)
00701 {
00702   /* local variables */
00703   char unsigned a_dot = 1,
00704     a_exp = 0,
00705     a_exp_sgn = 0,
00706     a_sgn = 1,
00707     a_div = 1;
00708   char c = 0;
00709   int l_exp = 0,
00710     sgn = 0,
00711     exp_sgn = 0;
00712   int n_char = 0,
00713     n_dig = 0,
00714     cn = 0;
00715   mpq_t den[2];
00716   mpq_init (den[0]);
00717   mpq_init (den[1]);
00718   mpq_set_ui (den[1], (unsigned long int)1, (unsigned long int)1);
00719   mpq_set_ui (den[0], (unsigned long int)0, (unsigned long int)1);
00720 
00721   /* now we read the string */
00722   c = str[n_char];
00723   while ((('0' <= c) && (c <= '9')) ||  /* allow to read digits */
00724          (a_dot && (c == '.')) || /* allow to read a dot point */
00725          (a_exp && (c == 'e' || c == 'E')) || /* allow an exponent marker */
00726          (a_sgn && (c == '+' || c == '-')) || /* allow a number sign */
00727          (a_div && c == '/') || /* allow the division sign */
00728          (a_exp_sgn && (c == '+' || c == '-')) /* allow sign for exponent */ )
00729   {
00730     switch (c)
00731     {
00732     case '0':
00733     case '1':
00734     case '2':
00735     case '3':
00736     case '4':
00737     case '5':
00738     case '6':
00739     case '7':
00740     case '8':
00741     case '9':
00742       /* if we haven't read the exponent then the digits bellongs to the mantisa
00743        * */
00744       if (a_exp || n_dig == 0)
00745       {
00746         if (!a_dot)
00747           mpz_mul_ui (mpq_denref (den[cn]), mpq_denref (den[cn]), (unsigned long int)10);
00748         mpz_mul_ui (mpq_numref (den[cn]), mpq_numref (den[cn]), (unsigned long int)10);
00749         mpz_add_ui (mpq_numref (den[cn]), mpq_numref (den[cn]),
00750                      (unsigned long int)(c - '0'));
00751         n_dig++;
00752         a_exp = 1;
00753       }
00754       /* otherwise, if we have read the exponent, the digits should go to the
00755        * exponent */
00756       else
00757       {
00758         l_exp = 10 * l_exp + c - '0';
00759         a_exp_sgn = 0;
00760       }
00761       a_sgn = 0;
00762       break;
00763     case '.':
00764       a_sgn = 0;
00765       a_dot = 0;
00766       a_sgn = 0;
00767       break;
00768     case '-':
00769       if (a_sgn)
00770         sgn = 1;
00771       else
00772         exp_sgn = 1;
00773     case '+':
00774       if (a_sgn)
00775         a_sgn = 0;
00776       if (a_exp_sgn)
00777         a_exp_sgn = 0;
00778       break;
00779     case 'e':
00780     case 'E':
00781       a_sgn = 0;
00782       a_exp = 0;
00783       a_exp_sgn = 1;
00784       break;
00785     case '/':
00786       if (exp_sgn)
00787         l_exp = -l_exp;
00788       if (l_exp > 0)
00789         while (l_exp--)
00790           mpz_mul_ui (mpq_numref (den[0]), mpq_numref (den[0]), (unsigned long int)10);
00791       else if (l_exp < 0)
00792       {
00793         l_exp = -l_exp;
00794         while (l_exp--)
00795           mpz_mul_ui (mpq_denref (den[0]), mpq_denref (den[0]), (unsigned long int)10);
00796       }
00797       if (sgn)
00798         mpz_neg (mpq_numref (den[0]), mpq_numref (den[0]));
00799       mpq_canonicalize (den[0]);
00800       mpq_set_ui (den[1], (unsigned long int)0, (unsigned long int)1);
00801       sgn = 0;
00802       exp_sgn = 0;
00803       l_exp = 0;
00804       a_div = 0;
00805       n_dig = 0;
00806       a_dot = 1;
00807       a_exp = 0;
00808       a_exp_sgn = 0;
00809       a_sgn = 1;
00810       cn = 1;
00811       break;
00812     }
00813     /* advance the reading character */
00814     c = str[++n_char];
00815   }
00816   if (n_char)
00817   {
00818     /* now expand the exponent of the denominator */
00819     if (exp_sgn)
00820       l_exp = -l_exp;
00821     if (l_exp > 0)
00822       while (l_exp--)
00823         mpz_mul_ui (mpq_numref (den[cn]), mpq_numref (den[cn]), (unsigned long int)10);
00824     else if (l_exp < 0)
00825     {
00826       l_exp = -l_exp;
00827       while (l_exp--)
00828         mpz_mul_ui (mpq_denref (den[cn]), mpq_denref (den[cn]), (unsigned long int)10);
00829     }
00830     /* check the sign of the whole number */
00831     if (sgn)
00832       mpz_neg (mpq_numref (den[cn]), mpq_numref (den[cn]));
00833     /* ending */
00834     mpq_canonicalize (den[0]);
00835     mpq_canonicalize (den[1]);
00836     mpq_div (var, den[0], den[1]);
00837   }
00838   mpq_clear (den[0]);
00839   mpq_clear (den[1]);
00840   return n_char;
00841 }
00842 #endif
00843 
00844 /* ========================================================================= */
00845 /** @} */