<!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>