00001
00002
00003
00004
00005
00006
00007
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
00030
00031 #ifdef SCS_TYPECPU_SPARC
00032 static const scs
00033
00034 scs_zer ={{0x00000000, 0x00000000, 0x00000000, 0x00000000},
00035 {{0, 0}}, 0, 1 },
00036
00037 scs_half={{0x02000000, 0x00000000, 0x00000000, 0x00000000},
00038 DB_ONE, -1, 1 },
00039
00040 scs_one ={{0x00000001, 0x00000000, 0x00000000, 0x00000000},
00041 DB_ONE, 0, 1 },
00042
00043 scs_two ={{0x00000002, 0x00000000, 0x00000000, 0x00000000},
00044 DB_ONE, 0, 1 };
00045 #else
00046 static const scs
00047
00048 scs_zer ={{0x00000000, 0x00000000, 0x00000000, 0x00000000,
00049 0x00000000, 0x00000000, 0x00000000, 0x00000000},
00050 {{0, 0}}, 0, 1 },
00051
00052 scs_half={{0x20000000, 0x00000000, 0x00000000, 0x00000000,
00053 0x00000000, 0x00000000, 0x00000000, 0x00000000},
00054 DB_ONE, -1, 1 },
00055
00056 scs_one ={{0x00000001, 0x00000000, 0x00000000, 0x00000000,
00057 0x00000000, 0x00000000, 0x00000000, 0x00000000},
00058 DB_ONE, 0, 1 },
00059
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
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118 #define SQRT_2 1.4142135623730950489e0
00119
00120
00121
00122
00123
00124
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
00134 if (nb.i[HI_ENDIAN] < 0x00100000){
00135 if (((nb.i[HI_ENDIAN] & 0x7fffffff)|nb.i[LO_ENDIAN])==0)
00136 return 1.0/0.0;
00137 if (nb.i[HI_ENDIAN] < 0)
00138 return (x-x)/0;
00139
00140
00141 E -= (SCS_NB_BITS*2);
00142 nb.d *=SCS_RADIX_TWO_DOUBLE;
00143
00144
00145 }
00146 if (nb.i[HI_ENDIAN] >= 0x7ff00000)
00147 return x+x;
00148
00149
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
00158
00159 nb2.d = nb.d + norm_number.d;
00160 i = (nb2.i[HI_ENDIAN] & 0x000fffff);
00161 i = i >> 16;
00162
00163
00164 wi.d = (11+i)*(double)0.6250e-1;
00165
00166
00167 nb.d -= wi.d;
00168
00169
00170 ti = i;
00171
00172
00173 scs_set_d(R, nb.d);
00174 scs_mul(R, R, table_inv_wi_ptr[i]);
00175
00176
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
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
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
00216 mpf_init2(R, 210);
00217 mpf_init2(sc_ln2_times_E, 210);
00218 mpf_init2(res1, 210);
00219
00220 nb.d = x;
00221
00222 if (nb.i[HI_ENDIAN] < 0x00100000){
00223 if (((nb.i[HI_ENDIAN] & 0x7fffffff)|nb.i[LO_ENDIAN])==0)
00224 return 1.0/0.0;
00225 if (nb.i[HI_ENDIAN] < 0)
00226 return (x-x)/0;
00227
00228
00229 E -= (SCS_NB_BITS*2);
00230 nb.d *=SCS_RADIX_TWO_DOUBLE;
00231
00232
00233 }
00234 if (nb.i[HI_ENDIAN] >= 0x7ff00000)
00235 return x+x;
00236
00237
00238
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
00247
00248 nb2.d = nb.d + norm_number.d;
00249 i = (nb2.i[HI_ENDIAN] & 0x000fffff);
00250 i = i >> 16;
00251
00252
00253 wi.d = (11+i)*(double)0.6250e-1;
00254
00255
00256 nb.d -= wi.d;
00257
00258
00259
00260 ti = i;
00261
00262
00263 mpf_set_d(R, nb.d);
00264 mpf_mul(R, R, mpf_table_inv_wi_ptr[i]);
00265
00266
00267
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
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
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
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
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
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