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