Main Page   Data Structures   File List   Data Fields   Globals  

test_log.c

00001 /*
00002  * Function to compute the logarithm with fully exact rounding
00003  *
00004  * Author : Defour David  (David.Defour@ens-lyon.fr)
00005  *
00006  * Date of creation : 13/05/2002   
00007  * Last Modified    : 17/05/2002
00008  */
00009 #include "scs.h"
00010 #include "scs_private.h"
00011 #include "log.h"
00012 #include <stdio.h>
00013 #ifdef HAVE_GMP_H
00014 #include <gmp.h>
00015 #include "tbx_timing.h"
00016 
00017 
00018 
00019 
00020 #define LOOPS 10000
00021 
00022 mpf_t mpf_sc_ln2_ptr;
00023 mpf_t mpf_constant_poly_ptr[20];
00024 mpf_t mpf_table_ti_ptr[13];
00025 mpf_t mpf_table_inv_wi_ptr[13];
00026 
00027 
00028 /*
00029  * 6) SCS CONSTANTS
00030  */
00031 #ifdef SCS_TYPECPU_SPARC
00032 static const scs
00033 /* 0   */
00034    scs_zer ={{0x00000000, 0x00000000, 0x00000000, 0x00000000},
00035              {{0, 0}},  0,   1 },
00036 /* 1/2 */
00037    scs_half={{0x02000000, 0x00000000, 0x00000000, 0x00000000},
00038              DB_ONE, -1,   1 },
00039 /*  1  */  
00040    scs_one ={{0x00000001, 0x00000000, 0x00000000, 0x00000000},
00041              DB_ONE,  0,   1 },
00042 /*  2  */
00043    scs_two ={{0x00000002, 0x00000000, 0x00000000, 0x00000000},
00044              DB_ONE,  0,   1 };
00045 #else
00046 static const scs
00047 /* 0   */
00048    scs_zer ={{0x00000000, 0x00000000, 0x00000000, 0x00000000,
00049              0x00000000, 0x00000000, 0x00000000, 0x00000000},
00050              {{0, 0}},  0,   1 },
00051 /* 1/2 */
00052    scs_half={{0x20000000, 0x00000000, 0x00000000, 0x00000000,
00053              0x00000000, 0x00000000, 0x00000000, 0x00000000},
00054              DB_ONE, -1,   1 },
00055 /*  1  */  
00056    scs_one ={{0x00000001, 0x00000000, 0x00000000, 0x00000000,
00057              0x00000000, 0x00000000, 0x00000000, 0x00000000},
00058              DB_ONE,  0,   1 },
00059 /*  2  */
00060    scs_two ={{0x00000002, 0x00000000, 0x00000000, 0x00000000,
00061              0x00000000, 0x00000000, 0x00000000, 0x00000000},
00062              DB_ONE,  0,   1 };
00063 #endif
00064 
00065 #define SCS_ZERO (scs_ptr)(&scs_zer)
00066 #define SCS_HALF (scs_ptr)(&scs_half)
00067 #define SCS_ONE  (scs_ptr)(&scs_one)
00068 #define SCS_TWO  (scs_ptr)(&scs_two)
00069 
00070 double sn_log(double); 
00071 double mpf_sn_log(double);
00072 
00073 
00074 /*
00075  *  1) First reduction: exponent extraction      
00076  *         E  
00077  *  x = 2^   .(1+f)    with  0 <= f < 1
00078  *
00079  *  log(x) = E.log(2) + log(1+f) where:
00080  *     - log(2)   is tabulated
00081  *     - log(1+f) need to be evalute 
00082  *  
00083  *
00084  *  2) Avoiding accuracy problem when E=-1 by testing
00085  *   
00086  *    if (1+f >= sqrt(2)) then 
00087  *        1+f = (1+f)/2;  E = E+1; 
00088  *    and,
00089  *        log(x) = (E+1).log(2) + log((1+f)/2)
00090  *
00091  *    so now:      sqrt(2)/2 <= (1+f) < sqrt(2)
00092  *
00093  *
00094  *  3) Second reduction: tabular reduction
00095  *                   -4  
00096  *   wi = 1 + i. 2^
00097  *                                   1  
00098  *   log(1+f) = log(wi) + log ( 1 + --- . (1 + f - wi) ) 
00099  *                                   wi 
00100  *
00101  *   then |(1+f-wi)/wi| <= 2^-5 if we use rounded to nearest.
00102  *
00103  *  4) Computation:
00104  *   a) Table lookup of: 
00105  *        - ti    = log(wi)
00106  *        - inv_wi = 1/(wi)
00107  *   b) Polynomial evaluation of:
00108  *        - P(R) ~ log(1 + R), where R = (1+f-wi) * inv_wi 
00109  *
00110  *                 -5 
00111  *   with  |R| < 2^
00112  *
00113  *
00114  *  5) Reconstruction:
00115  *   log(x) = E.log(2) + t_i + P(R)
00116  *
00117  */
00118 #define SQRT_2 1.4142135623730950489e0 
00119 
00120 
00121 
00122 /*************************************************************
00123  *************************************************************
00124  *               ROUNDED  TO NEAREST
00125  *************************************************************
00126  *************************************************************/
00127 double sn_log(double x){ 
00128   scs_t R, sc_ln2_times_E, res1;
00129   scs_db_number nb, nb2, wi, resd;
00130   int ti, i, E=0;
00131 
00132   nb.d = x;
00133   /* Filter cases */
00134   if (nb.i[HI_ENDIAN] < 0x00100000){        /* x < 2^(-1022)    */
00135     if (((nb.i[HI_ENDIAN] & 0x7fffffff)|nb.i[LO_ENDIAN])==0)
00136       return 1.0/0.0;                       /* log(+/-0) = -Inf */
00137     if (nb.i[HI_ENDIAN] < 0) 
00138       return (x-x)/0;                       /* log(-x) = Nan    */
00139 
00140     /* Subnormal number */
00141     E    -= (SCS_NB_BITS*2); /* keep in mind that x is a subnormal number */ 
00142     nb.d *=SCS_RADIX_TWO_DOUBLE;  /* make x as normal number     */         
00143     /* We may just want add 2 to the scs number.index */
00144     /* may be .... we will see */
00145   }
00146   if (nb.i[HI_ENDIAN] >= 0x7ff00000)
00147     return x+x;                             /* Inf or Nan       */
00148 
00149   /* find n, nb.d such that sqrt(2)/2 < nb.d < sqrt(2) */
00150   E += (nb.i[HI_ENDIAN]>>20)-1023;
00151   nb.i[HI_ENDIAN] =  (nb.i[HI_ENDIAN] & 0x000fffff) | 0x3ff00000;
00152   if (nb.d > SQRT_2){
00153     nb.d *= 0.5;
00154     E++;
00155   }
00156 
00157   /* to normalize nb.d and round to nearest      */
00158   /* + (1-trunc(sqrt(2.)/2 * 2^(4))*2^(-4) )+2.^(-(4+1))*/ 
00159   nb2.d = nb.d + norm_number.d; 
00160   i = (nb2.i[HI_ENDIAN] & 0x000fffff);
00161   i = i >> 16; /* 0<= i <=11 */
00162   
00163 
00164   wi.d = (11+i)*(double)0.6250e-1;
00165 
00166   /* (1+f-w_i) */
00167   nb.d -= wi.d; 
00168 
00169   /* Table reduction */
00170   ti = i; 
00171 
00172   /* R = (1+f-w_i)/w_i */
00173   scs_set_d(R, nb.d);
00174   scs_mul(R, R, table_inv_wi_ptr[i]);
00175 
00176   /* sc_ln2_times_E = E*log(2)  */
00177   scs_set(sc_ln2_times_E, sc_ln2_ptr);
00178 
00179   if (E >= 0){
00180     scs_mul_ui(sc_ln2_times_E, E);
00181   }else{
00182     scs_mul_ui(sc_ln2_times_E, -E);
00183     sc_ln2_times_E->sign = -1;
00184   }
00185 
00186   /*
00187    * Polynomial evaluation of log(1 + R) with an error less than 2^(-130)
00188    */
00189   scs_mul(res1, constant_poly_ptr[0], R);
00190   for(i=1; i<20; i++){
00191     scs_add(res1, constant_poly_ptr[i], res1);
00192     scs_mul(res1, res1, R);
00193   }
00194   scs_add(res1, res1, table_ti_ptr[ti]);
00195   scs_add(res1, res1, sc_ln2_times_E);  
00196   
00197 
00198   scs_get_d(&resd.d, res1);  
00199 
00200   return resd.d;
00201 }
00202 
00203 
00204 
00205 /*************************************************************
00206  *************************************************************
00207  *               ROUNDED  TO NEAREST
00208  *************************************************************
00209  *************************************************************/
00210 double mpf_sn_log(double x){ 
00211     mpf_t R, sc_ln2_times_E, res1;
00212     scs_db_number nb, nb2, wi, resd;
00213     int ti, i, E=0;
00214 
00215     /* Memory allocation */
00216     mpf_init2(R, 210);
00217     mpf_init2(sc_ln2_times_E, 210);
00218     mpf_init2(res1, 210);
00219 
00220   nb.d = x;
00221   /* Filter cases */
00222   if (nb.i[HI_ENDIAN] < 0x00100000){        /* x < 2^(-1022)    */
00223     if (((nb.i[HI_ENDIAN] & 0x7fffffff)|nb.i[LO_ENDIAN])==0)
00224       return 1.0/0.0;                       /* log(+/-0) = -Inf */
00225     if (nb.i[HI_ENDIAN] < 0) 
00226       return (x-x)/0;                       /* log(-x) = Nan    */
00227 
00228     /* Subnormal number */
00229     E    -= (SCS_NB_BITS*2); /* keep in mind that x is a subnormal number */ 
00230     nb.d *=SCS_RADIX_TWO_DOUBLE;  /* make x as normal number     */         
00231     /* We may just want add 2 to the scs number.index */
00232     /* may be .... we will see */
00233   }
00234   if (nb.i[HI_ENDIAN] >= 0x7ff00000)
00235     return x+x;                             /* Inf or Nan       */
00236 
00237 
00238   /* find n, nb.d such that sqrt(2)/2 < nb.d < sqrt(2) */
00239   E += (nb.i[HI_ENDIAN]>>20)-1023;
00240   nb.i[HI_ENDIAN] =  (nb.i[HI_ENDIAN] & 0x000fffff) | 0x3ff00000;
00241   if (nb.d > SQRT_2){
00242     nb.d *= 0.5;
00243     E++;
00244   }
00245 
00246   /* to normalize nb.d and round to nearest      */
00247   /* + (1-trunc(sqrt(2.)/2 * 2^(4))*2^(-4) )+2.^(-(4+1))*/ 
00248   nb2.d = nb.d + norm_number.d; 
00249   i = (nb2.i[HI_ENDIAN] & 0x000fffff);
00250   i = i >> 16; /* 0<= i <=11 */
00251   
00252 
00253   wi.d = (11+i)*(double)0.6250e-1;
00254 
00255   /* (1+f-w_i) */
00256   nb.d -= wi.d; 
00257 
00258 
00259   /* Table reduction */
00260   ti = i; 
00261  
00262   /* R = (1+f-w_i)/w_i */
00263   mpf_set_d(R, nb.d);
00264   mpf_mul(R, R, mpf_table_inv_wi_ptr[i]);
00265  
00266 
00267   /* sc_ln2_times_E = E*log(2)  */
00268   mpf_set(sc_ln2_times_E, mpf_sc_ln2_ptr);
00269 
00270 
00271   if (E >= 0){
00272     mpf_mul_ui(sc_ln2_times_E, sc_ln2_times_E, E);
00273   }else{
00274     mpf_mul_ui(sc_ln2_times_E, sc_ln2_times_E, -E);
00275     mpf_neg(sc_ln2_times_E, sc_ln2_times_E);
00276   }
00277 
00278 
00279   /*
00280    * Polynomial evaluation of log(1 + R) with an error less than 2^(-130)
00281    */
00282   mpf_mul(res1, mpf_constant_poly_ptr[0], R);
00283   for(i=1; i<20; i++){
00284     mpf_add(res1, mpf_constant_poly_ptr[i], res1);
00285     mpf_mul(res1, res1, R);
00286   }
00287   mpf_add(res1, res1, mpf_table_ti_ptr[ti]);
00288   mpf_add(res1, res1, sc_ln2_times_E);  
00289 
00290   resd.d = mpf_get_d(res1);  
00291 
00292   /* Free Memory */
00293   mpf_clear(R);  mpf_clear(sc_ln2_times_E);  mpf_clear(res1);
00294 
00295   return resd.d;
00296 }
00297 
00298 
00299 
00300 void free_mpf_cst(){
00301     int i;
00302 
00303     for (i=0; i<13; i++) mpf_clear(mpf_table_ti_ptr[i]);
00304     for (i=0; i<13; i++) mpf_clear(mpf_table_inv_wi_ptr[i]);
00305     for (i=0; i<20; i++) mpf_clear(mpf_constant_poly_ptr[i]);
00306     mpf_clear(mpf_sc_ln2_ptr);
00307 }
00308 
00309 void init_mpf_cst(){
00310 
00311     /*
00312      * mpf constant initialisation 
00313      */
00314  mpf_init_set_str(mpf_sc_ln2_ptr,".6931471805599453094172321214581765680755001343602552541206800094933936",10);
00315  mpf_init_set_str(mpf_table_ti_ptr[0],"-.3746934494414106936069849078675769724802936835036038412641523288430009",10);
00316  mpf_init_set_str(mpf_table_ti_ptr[1],"-.2876820724517809274392190059938274315035097108977610565066656853492930",10);
00317  mpf_init_set_str(mpf_table_ti_ptr[2],"-.2076393647782445016154410442673876674967325926808139000636745273101098",10);
00318  mpf_init_set_str(mpf_table_ti_ptr[3],"-.1335313926245226231463436209313499745894156734989045739026498785426010",10);
00319  mpf_init_set_str(mpf_table_ti_ptr[4],"-.06453852113757117167292391568399292812890862534975384283537781286190121",10);
00320  mpf_init_set_str(mpf_table_ti_ptr[5],"0.0",10);
00321  mpf_init_set_str(mpf_table_ti_ptr[6],".06062462181643484258060613204042026328620247514472377081451769990871809",10);
00322  mpf_init_set_str(mpf_table_ti_ptr[7],".1177830356563834545387941094705217050684807125647331411073486387948077",10);
00323  mpf_init_set_str(mpf_table_ti_ptr[8],".1718502569266592223400989460551472649353787238581078020552401984357182",10);
00324  mpf_init_set_str(mpf_table_ti_ptr[9],".2231435513142097557662950903098345033746010855480072136712878724873917",10);
00325  mpf_init_set_str(mpf_table_ti_ptr[10],".2719337154836417588316694945329991619825747499635896237113644456014997",10);
00326  mpf_init_set_str(mpf_table_ti_ptr[11],".3184537311185346158102472135905995955952064508566514128565276806503928",10);
00327  mpf_init_set_str(mpf_table_ti_ptr[12],".3629054936893684531378243459774898461403797773994147255159153395094188",10);
00328   
00329  mpf_init_set_str(mpf_table_inv_wi_ptr[0],"1.454545454545454545454545454545454545454545454545454545454545454545455",10);
00330  mpf_init_set_str(mpf_table_inv_wi_ptr[1],"1.333333333333333333333333333333333333333333333333333333333333333333333",10);
00331  mpf_init_set_str(mpf_table_inv_wi_ptr[2],"1.230769230769230769230769230769230769230769230769230769230769230769231",10);
00332  mpf_init_set_str(mpf_table_inv_wi_ptr[3],"1.142857142857142857142857142857142857142857142857142857142857142857143",10);
00333  mpf_init_set_str(mpf_table_inv_wi_ptr[4],"1.066666666666666666666666666666666666666666666666666666666666666666667",10);
00334  mpf_init_set_str(mpf_table_inv_wi_ptr[5],"1.0",10);
00335  mpf_init_set_str(mpf_table_inv_wi_ptr[6],".9411764705882352941176470588235294117647058823529411764705882352941176",10);
00336  mpf_init_set_str(mpf_table_inv_wi_ptr[7],".8888888888888888888888888888888888888888888888888888888888888888888889",10);
00337  mpf_init_set_str(mpf_table_inv_wi_ptr[8],".8421052631578947368421052631578947368421052631578947368421052631578947",10);
00338  mpf_init_set_str(mpf_table_inv_wi_ptr[9],".8000000000000000000000000000000000000000000000000000000000000000000000",10);
00339  mpf_init_set_str(mpf_table_inv_wi_ptr[10],".7619047619047619047619047619047619047619047619047619047619047619047619",10);
00340  mpf_init_set_str(mpf_table_inv_wi_ptr[11],".7272727272727272727272727272727272727272727272727272727272727272727273",10);
00341  mpf_init_set_str(mpf_table_inv_wi_ptr[12],".6956521739130434782608695652173913043478260869565217391304347826086957",10);
00342   
00343 
00344 
00345  mpf_init_set_str(mpf_constant_poly_ptr[0],"-.5023367294567568078075892234129085015737046673117638148835793879900684e-1",10);
00346  mpf_init_set_str(mpf_constant_poly_ptr[1],".5286469364636800162451931056605008699849205395651010015123349866702438e-1",10);
00347  mpf_init_set_str(mpf_constant_poly_ptr[2],"-.5555504165240301671600703643945025506173351812330050928639639013749408e-1",10);
00348  mpf_init_set_str(mpf_constant_poly_ptr[3],".5882304523791230393387514237891031448778320732847321758009922427128675e-1",10);
00349  mpf_init_set_str(mpf_constant_poly_ptr[4],"-.6250000063225403563289268352338121550211573212295444469323271379834995e-1",10);
00350  mpf_init_set_str(mpf_constant_poly_ptr[5],".6666666722317130353982967043683210254854513387493745794307824211890102e-1",10);
00351  mpf_init_set_str(mpf_constant_poly_ptr[6],"-.7142857142809460344869308613352646790467720996336436837300223496769771e-1",10);
00352  mpf_init_set_str(mpf_constant_poly_ptr[7],".7692307692269045639218260467556323317594125412187843530035925183575360e-1",10);
00353  mpf_init_set_str(mpf_constant_poly_ptr[8],"-.8333333333333356037828752693534976595948243524382699349340678490229276e-1",10);
00354  mpf_init_set_str(mpf_constant_poly_ptr[9],".9090909090909107517931824416149710386096279143853141933153563186364239e-1",10);
00355  mpf_init_set_str(mpf_constant_poly_ptr[10],"-.9999999999999999993224242653285578758483071952239027712038970659952908e-1",10);
00356  mpf_init_set_str(mpf_constant_poly_ptr[11],".1111111111111111110676605039270706067112823008194317757408530067415460",10);
00357  mpf_init_set_str(mpf_constant_poly_ptr[12],"-.1250000000000000000000121547531034901577892307419239533855249184555657",10);
00358  mpf_init_set_str(mpf_constant_poly_ptr[13],".1428571428571428571428636714934431388385979630564777290210429555695253",10);
00359  mpf_init_set_str(mpf_constant_poly_ptr[14],"-.1666666666666666666666666654681759449880186388123105519286024077940397",10);
00360  mpf_init_set_str(mpf_constant_poly_ptr[15],".1999999999999999999999999995018689489549856935258632875129946570831422",10);
00361  mpf_init_set_str(mpf_constant_poly_ptr[16],"-.2500000000000000000000000000000541884832055314060603150374550752650562",10);
00362  mpf_init_set_str(mpf_constant_poly_ptr[17],".3333333333333333333333333333333480752710184240674027421784076979756719",10);
00363  mpf_init_set_str(mpf_constant_poly_ptr[18],"-.4999999999999999999999999999999999992783495160517279217122998393096071",10);
00364  mpf_init_set_str(mpf_constant_poly_ptr[19],".9999999999999999999999999999999999999280145118565650165317802247440368",10);
00365 
00366 }
00367 
00368 
00369 /*
00370  * Used to mesure time taken by the instruction "name"
00371  */
00372 #define TST_FCT(char, name)\
00373          TBX_GET_TICK(t1);\
00374          for(i=0; i< LOOPS-1; i++){ \
00375            name;\
00376            } \
00377          TBX_GET_TICK(t2);\
00378          deltat = TBX_TICK_RAW_DIFF(t1, t2); \
00379          printf("%28s :  %lld ticks \n", char, deltat);
00380 
00381 
00382 
00383 int main(){
00384     /* table storing random number */
00385     double  tableau[LOOPS];
00386     double res[LOOPS];
00387     mpf_t  TAB_MPF[LOOPS];
00388     int i;
00389     tbx_tick_t   t1, t2; 
00390     unsigned long long deltat;
00391 
00392     printf("mpf constant initialisation ... \n");
00393     init_mpf_cst();
00394     printf("Finished ... \n");
00395 
00396     srand(42);  
00397     for(i=0; i<LOOPS; i++) tableau[i]=rand();
00398     printf("End of random number generation \n");
00399 
00400     for(i=0; i< LOOPS-1; i++){ 
00401       res[i]=sn_log(tableau[i]);
00402     } 
00403     for(i=0; i< LOOPS-1; i++){ 
00404       res[i]=mpf_sn_log(tableau[i]);
00405     } 
00406 
00407 
00408     TBX_GET_TICK(t1);
00409     for(i=0; i< LOOPS-1; i++){ 
00410       res[i]=sn_log(tableau[i]);
00411     } 
00412     TBX_GET_TICK(t2);
00413     deltat = TBX_TICK_RAW_DIFF(t1, t2); 
00414     printf("%28s :  %lld ticks \n", "scs_log ", deltat);
00415     printf("\n");
00416 
00417 
00418 
00419 
00420     TBX_GET_TICK(t1);
00421     for(i=0; i< LOOPS-1; i++){ 
00422       res[i]=mpf_sn_log(tableau[i]);
00423     } 
00424     TBX_GET_TICK(t2);
00425     deltat = TBX_TICK_RAW_DIFF(t1, t2); 
00426     printf("%28s :  %lld ticks \n", "mpf_log " , deltat);
00427 
00428 
00429     free_mpf_cst();
00430 }
00431 
00432 #else
00433 main(){
00434     fprintf(stderr,"No GMP detected on your system\n");
00435 }
00436 #endif /*HAVE_GMP_H*/

Generated on Tue Jun 17 10:15:51 2003 for SCSLib by doxygen1.2.15