<!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_accuracy.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_accuracy.c</h1><a href="test__accuracy_8c.html">Go to the documentation of this file.</a><div class="fragment"><pre>00001 <font class="comment">/** Test of SCS library against MPFR</font> 00002 <font class="comment">@file tests/test_accuracy.c</font> 00003 <font class="comment"></font> 00004 <font class="comment">@author David Defour David.Defour@ens-lyon.fr</font> 00005 <font class="comment">@author Florent de Dinechin Florent.de.Dinechin@ens-lyon.fr </font> 00006 <font class="comment"></font> 00007 <font class="comment">This file is part of the SCS library.</font> 00008 <font class="comment">*/</font> 00009 00010 <font class="comment">/*</font> 00011 <font class="comment">Copyright (C) 2002 David Defour and Florent de Dinechin</font> 00012 <font class="comment"></font> 00013 <font class="comment"> This library is free software; you can redistribute it and/or</font> 00014 <font class="comment"> modify it under the terms of the GNU Lesser General Public</font> 00015 <font class="comment"> License as published by the Free Software Foundation; either</font> 00016 <font class="comment"> version 2.1 of the License, or (at your option) any later version.</font> 00017 <font class="comment"></font> 00018 <font class="comment"> This library is distributed in the hope that it will be useful,</font> 00019 <font class="comment"> but WITHOUT ANY WARRANTY; without even the implied warranty of</font> 00020 <font class="comment"> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU</font> 00021 <font class="comment"> Lesser General Public License for more details.</font> 00022 <font class="comment"></font> 00023 <font class="comment"> You should have received a copy of the GNU Lesser General Public</font> 00024 <font class="comment"> License along with this library; if not, write to the Free Software</font> 00025 <font class="comment"> Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA</font> 00026 <font class="comment"></font> 00027 <font class="comment"> */</font> 00028 <font class="preprocessor">#include "<a class="code" href="scs_8h.html">scs.h</a>"</font> 00029 <font class="preprocessor">#include "<a class="code" href="scs__private_8h.html">scs_private.h</a>"</font> 00030 <font class="preprocessor">#include <stdio.h></font> 00031 <font class="preprocessor">#include <stdlib.h></font> 00032 <font class="preprocessor">#include <math.h></font> 00033 00034 00035 <font class="comment">/* Compile only if mpfr is present */</font> 00036 <font class="preprocessor">#ifdef HAVE_MPFR_H</font> 00037 <font class="preprocessor"></font> 00038 <font class="preprocessor">#include <mpfr.h></font> 00039 00040 00041 void (* test_scs_fct) () = NULL; 00042 int (* test_mpfr_fct) () = NULL; 00043 mp_rnd_t ROUNDING = GMP_RNDN; 00044 00045 00046 <font class="comment">/* Function defined in scs2MPF.c */</font> 00047 <font class="comment">//void scs2MPF(scs_t, mpf_t);</font> 00048 00049 00050 <font class="keywordtype">void</font> usage(<font class="keywordtype">char</font> *namefct){ 00051 printf(<font class="stringliteral">"Usage :\n \t %s n [k] \n\n"</font>, namefct); 00052 printf(<font class="stringliteral">"Where k is the number of random tests to generate (default 1000) \n"</font>); 00053 printf(<font class="stringliteral">"and n is the function to test \n"</font>); 00054 printf(<font class="stringliteral">" 1 : scs_add \n"</font>); 00055 printf(<font class="stringliteral">" 2 : scs_sub \n"</font>); 00056 printf(<font class="stringliteral">" 3 : scs_mul \n"</font>); 00057 printf(<font class="stringliteral">" 4 : scs_div \n"</font>); 00058 printf(<font class="stringliteral">" 5 : scs_get_d \n"</font>); 00059 printf(<font class="stringliteral">" 6 : scs_get_d_minf \n"</font>); 00060 printf(<font class="stringliteral">" 7 : scs_get_d_pinf \n"</font>); 00061 printf(<font class="stringliteral">" 8 : scs_get_d_zero \n"</font>); 00062 printf(<font class="stringliteral">" 9 : scs_set \n"</font>); 00063 printf(<font class="stringliteral">" 10: scs_square \n"</font>); 00064 00065 printf(<font class="stringliteral">"\n"</font>); 00066 exit (1); 00067 } 00068 00069 00070 <font class="keywordtype">void</font> test_square(<font class="keywordtype">int</font> bcl){ 00071 <a class="code" href="scs_8h.html#a3">scs_t</a> scx, scr, scx_maxerr, scr_maxerr; 00072 mpfr_t mpx, mpr,mpscr,mperr; 00073 <font class="keywordtype">int</font> j; 00074 <font class="keywordtype">double</font> err, maxerr, avgerr; 00075 00076 mpfr_init2(mpx, SCS_NB_WORDS*SCS_NB_BITS*2); 00077 mpfr_init2(mpr, SCS_NB_WORDS*SCS_NB_BITS*2); 00078 mpfr_init2(mpscr, SCS_NB_WORDS*SCS_NB_BITS*2); 00079 mpfr_init2(mperr, SCS_NB_WORDS*SCS_NB_BITS*2); 00080 maxerr = 0; avgerr=0; 00081 00082 <font class="keywordflow">for</font>(j=0; j<bcl; j++){ 00083 scs_rand(scx, 20); 00084 <font class="preprocessor"> #if 1 </font> 00085 <font class="preprocessor"></font> <font class="comment">/* set worst case */</font> 00086 scx -> h_word[0] = 1; 00087 <font class="preprocessor"> #endif</font> 00088 <font class="preprocessor"></font> scs_get_mpfr(scx, mpx); 00089 00090 <a class="code" href="multiplication__scs_8c.html#a3">scs_square</a>(scr, scx); 00091 mpfr_mul(mpr, mpx, mpx, GMP_RNDN); 00092 scs_get_mpfr(scr, mpscr); 00093 00094 mpfr_sub(mperr, mpscr, mpr, GMP_RNDN); 00095 mpfr_div(mperr, mperr, mpr, GMP_RNDN); 00096 mpfr_abs(mperr, mperr, GMP_RNDN); 00097 00098 err = mpfr_get_d(mperr, GMP_RNDN); 00099 <font class="keywordflow">if</font> (err > maxerr){ 00100 maxerr = err; 00101 <a class="code" href="addition__scs_8c.html#a0">scs_set</a>(scx_maxerr, scx); 00102 <a class="code" href="addition__scs_8c.html#a0">scs_set</a>(scr_maxerr, scr); 00103 } 00104 avgerr += err; 00105 } 00106 avgerr = avgerr / bcl; 00107 printf(<font class="stringliteral">"Average error : %e nb bit : %d \n"</font>, 00108 avgerr, -(<font class="keywordtype">int</font>)(log(avgerr)/log(2.))); 00109 printf(<font class="stringliteral">"Max error : %e nb bit : %d \n"</font>, 00110 maxerr, -(<font class="keywordtype">int</font>)(log(maxerr)/log(2.))); 00111 00112 printf(<font class="stringliteral">"Argument giving max error : \n"</font>); scs_get_std(scx_maxerr); 00113 printf(<font class="stringliteral">"Result with max error: \n"</font>); scs_get_std(scr_maxerr); 00114 00115 mpfr_clear(mpx); mpfr_clear(mpr); mpfr_clear(mperr); 00116 mpfr_clear(mpscr); 00117 } 00118 00119 00120 00121 00122 <font class="keywordtype">void</font> test_one_arg(<font class="keywordtype">int</font> bcl){ 00123 <a class="code" href="scs_8h.html#a3">scs_t</a> sc1, scex; 00124 mpfr_t mpf1, mpfex; 00125 <font class="keywordtype">double</font> scsd, mpfrd, scsdmax, mpfrdmax; 00126 <font class="keywordtype">int</font> j, nberr; 00127 00128 00129 mpfr_init2(mpf1, SCS_NB_WORDS*SCS_NB_BITS*2); 00130 mpfr_init2(mpfex, SCS_NB_WORDS*SCS_NB_BITS*2); 00131 scsd = 0; mpfrd = 0; nberr = 0; scsdmax = 0; 00132 00133 <font class="keywordflow">for</font>(j=0; j<bcl; j++){ 00134 scs_rand(sc1, 1500); <font class="comment">/* test some special cases */</font> 00135 test_scs_fct(&scsd, sc1); 00136 00137 scs_get_mpfr(sc1, mpf1); 00138 00139 mpfrd = mpfr_get_d(mpf1, ROUNDING); 00140 00141 <font class="keywordflow">if</font> (mpfrd != scsd){ 00142 <a class="code" href="addition__scs_8c.html#a0">scs_set</a>(scex, sc1); 00143 nberr ++;} 00144 } 00145 printf(<font class="stringliteral">"Number of misrounds: %d\n"</font>, nberr); 00146 <font class="keywordflow">if</font> (nberr){ 00147 printf(<font class="stringliteral">"Original number :\n"</font>); 00148 scs_get_std(scex); 00149 scs_get_mpfr(scex, mpfex); 00150 mpfr_out_str(stdout, 2, 150, mpfex, ROUNDING); printf(<font class="stringliteral">"\n"</font>); 00151 00152 printf(<font class="stringliteral">"SCS rounding : \n"</font>); 00153 test_scs_fct(&scsd, scex); 00154 printf(<font class="stringliteral">"%.40e\n"</font>,scsd); 00155 printf(<font class="stringliteral">"In binary : \n"</font>); 00156 mpfr_set_d(mpf1, scsd, ROUNDING); 00157 mpfr_out_str(stdout, 2, 150, mpf1, ROUNDING); printf(<font class="stringliteral">"\n"</font>); 00158 00159 printf(<font class="stringliteral">"MPFR rounding : \n"</font>); 00160 mpfrd = mpfr_get_d(mpfex, ROUNDING); 00161 printf(<font class="stringliteral">"%.40e\n"</font>,mpfrd); 00162 printf(<font class="stringliteral">"In binary : \n"</font>); 00163 mpfr_set_d(mpf1, mpfrd, ROUNDING); 00164 mpfr_out_str(stdout, 2, 150, mpf1, ROUNDING); printf(<font class="stringliteral">"\n"</font>); 00165 } 00166 00167 mpfr_clear(mpf1); mpfr_clear(mpfex); 00168 } 00169 00170 00171 <font class="keywordtype">void</font> test_two_arg(<font class="keywordtype">int</font> bcl){ 00172 <a class="code" href="scs_8h.html#a3">scs_t</a> sc1, sc2, sc3, scm1, scm2, scm3; 00173 mpfr_t mp1, mp2, mp3, mp4, mp5; 00174 <font class="keywordtype">double</font> d1, d2, max; 00175 <font class="keywordtype">int</font> j; 00176 00177 mpfr_init2(mp1, SCS_NB_WORDS*SCS_NB_BITS*2); 00178 mpfr_init2(mp2, SCS_NB_WORDS*SCS_NB_BITS*2); 00179 mpfr_init2(mp3, SCS_NB_WORDS*SCS_NB_BITS*2); 00180 mpfr_init2(mp4, SCS_NB_WORDS*SCS_NB_BITS*2); 00181 mpfr_init2(mp5, SCS_NB_WORDS*SCS_NB_BITS*2); 00182 d1 = 0; d2 = 0; max = 0; 00183 00184 00185 <font class="keywordflow">for</font>(j=0; j<bcl; j++){ 00186 scs_rand(sc1, 20); scs_rand(sc2, 20); 00187 <font class="comment">/* You get most worst cases by imposing the following: */</font> 00188 <font class="preprocessor"> #if 1</font> 00189 <font class="preprocessor"></font> sc1 -> h_word[0] = 1; 00190 sc2 -> h_word[0] = 1; 00191 sc1 -> h_word[1] = 1; 00192 sc2 -> h_word[1] = 1; 00193 <font class="preprocessor"> #endif</font> 00194 <font class="preprocessor"></font> test_scs_fct(sc3, sc1, sc2); 00195 00196 scs_get_mpfr(sc1, mp1); scs_get_mpfr(sc2, mp2); scs_get_mpfr(sc3, mp3); 00197 00198 test_mpfr_fct(mp4, mp1, mp2, GMP_RNDN); 00199 00200 <font class="comment">/* Error Computation */</font> 00201 mpfr_sub(mp5, mp4, mp3, GMP_RNDN); mpfr_div(mp5, mp5, mp4, GMP_RNDN); mpfr_abs(mp5, mp5, GMP_RNDN); 00202 00203 d2 = mpfr_get_d(mp5, GMP_RNDN); 00204 <font class="keywordflow">if</font> (d2 > max){ 00205 max = d2; 00206 <a class="code" href="addition__scs_8c.html#a0">scs_set</a>(scm1, sc1); 00207 <a class="code" href="addition__scs_8c.html#a0">scs_set</a>(scm2, sc2); 00208 <a class="code" href="addition__scs_8c.html#a0">scs_set</a>(scm3, sc3); } 00209 00210 d1 += d2; 00211 } 00212 00213 printf(<font class="stringliteral">"Average error : %e nb bit : %d \n"</font>, d1/bcl, -(<font class="keywordtype">int</font>)(log(d1/bcl)/log(2.))); 00214 printf(<font class="stringliteral">"Max error : %e nb bit : %d \n"</font>, max, -(<font class="keywordtype">int</font>)(log(max)/log(2.))); 00215 00216 00217 printf(<font class="stringliteral">"First Number : \n"</font>); scs_get_std(scm1); 00218 printf(<font class="stringliteral">"Second Number : \n"</font>); scs_get_std(scm2); 00219 printf(<font class="stringliteral">"Result Number : \n"</font>); scs_get_std(scm3); 00220 00221 mpfr_clear(mp1); mpfr_clear(mp2); mpfr_clear(mp3); 00222 mpfr_clear(mp4); mpfr_clear(mp5); 00223 00224 } 00225 00226 00227 <font class="keywordtype">void</font> test_scs_set_d(<font class="keywordtype">int</font> loops) { 00228 scs_db_number d; 00229 <font class="keywordtype">double</font> d1, d2; 00230 <a class="code" href="scs_8h.html#a3">scs_t</a> s; 00231 <font class="keywordtype">int</font> i,j, errors, ep; 00232 <font class="keywordtype">double</font> errorcases1[20]; 00233 <font class="keywordtype">double</font> errorcases2[20]; 00234 00235 errors=0; ep=0; 00236 00237 <font class="keywordflow">for</font> (i=0; i<loops; i++) { 00238 d.l = (rand() & 0x000000ff); 00239 <font class="keywordflow">for</font>(j=0; j<(<font class="keyword">sizeof</font>(scs_db_number)); j++){ 00240 d.l = d.l << 8; 00241 d.l += (rand() & 0x000000ff ); 00242 } 00243 d1 = d.d; 00244 <a class="code" href="double2scs_8c.html#a0">scs_set_d</a>(s,d1); 00245 <a class="code" href="scs2double_8c.html#a0">scs_get_d</a>(&d2,s); 00246 <font class="keywordflow">if</font>((d1!=d2) && (!(!(d1==d1)))) { 00247 <font class="comment">/* second test to filter out NaNs, for which the first is always</font> 00248 <font class="comment"> true */</font> 00249 errors++; 00250 <font class="keywordflow">if</font>(ep<20) { 00251 errorcases1[ep]=d1; 00252 errorcases2[ep]=d2; 00253 ep++; 00254 } 00255 } 00256 } 00257 printf(<font class="stringliteral">"Number of errors: %d\n"</font>, errors); 00258 <font class="keywordflow">if</font>(ep!=0) { 00259 <font class="keywordflow">for</font>(i=0; i<ep; i++) 00260 printf(<font class="stringliteral">"example:\t%e\t%e\n"</font>, errorcases1[i], errorcases2[i]); 00261 } 00262 } 00263 00264 00265 <font class="comment">/*</font> 00266 <font class="comment"> * Fct de test . </font> 00267 <font class="comment"> */</font> 00268 <font class="keywordtype">int</font> main(<font class="keywordtype">int</font> argc, <font class="keywordtype">char</font> *argv[]){ 00269 <font class="keywordtype">int</font> fct, bcl; 00270 00271 00272 <font class="keywordflow">if</font> (argc < 2) usage(argv[0]); 00273 00274 fct = atoi(argv[1]); 00275 bcl = (argc >= 3)? atoi(argv[2]) : 1000; 00276 00277 00278 <font class="keywordflow">switch</font>(fct){ 00279 <font class="keywordflow">case</font> 1: 00280 test_scs_fct = scs_add; test_mpfr_fct = mpfr_add; 00281 test_two_arg(bcl); 00282 <font class="keywordflow">break</font>; 00283 <font class="keywordflow">case</font> 2: 00284 test_scs_fct = scs_sub; test_mpfr_fct = mpfr_sub; 00285 test_two_arg(bcl); 00286 <font class="keywordflow">break</font>; 00287 <font class="keywordflow">case</font> 3: 00288 test_scs_fct = scs_mul; test_mpfr_fct = mpfr_mul; 00289 test_two_arg(bcl); 00290 <font class="keywordflow">break</font>; 00291 <font class="keywordflow">case</font> 4: 00292 test_scs_fct = scs_div; test_mpfr_fct = mpfr_div; 00293 test_two_arg(bcl); 00294 <font class="keywordflow">break</font>; 00295 <font class="keywordflow">case</font> 5: 00296 test_scs_fct = scs_get_d; ROUNDING = GMP_RNDN; 00297 test_one_arg(bcl); 00298 <font class="keywordflow">break</font>; 00299 <font class="keywordflow">case</font> 6: 00300 test_scs_fct = scs_get_d_minf; ROUNDING = GMP_RNDD; 00301 test_one_arg(bcl); 00302 <font class="keywordflow">break</font>; 00303 <font class="keywordflow">case</font> 7: 00304 test_scs_fct = scs_get_d_pinf; ROUNDING = GMP_RNDU; 00305 test_one_arg(bcl); 00306 <font class="keywordflow">break</font>; 00307 <font class="keywordflow">case</font> 8: 00308 test_scs_fct = scs_get_d_zero; ROUNDING = GMP_RNDZ; 00309 test_one_arg(bcl); 00310 <font class="keywordflow">break</font>; 00311 <font class="keywordflow">case</font> 9: 00312 test_scs_set_d(bcl); 00313 <font class="keywordflow">break</font>; 00314 <font class="keywordflow">case</font> 10: 00315 test_square(bcl); 00316 <font class="keywordflow">break</font>; 00317 <font class="keywordflow">default</font> : 00318 fprintf(stderr,<font class="stringliteral">"Unknown Function \n"</font>); 00319 usage(argv[0]); 00320 } 00321 00322 <font class="keywordflow">return</font> 0; 00323 } 00324 <font class="preprocessor">#else</font> 00325 <font class="preprocessor"></font>main(){ 00326 fprintf(stderr,<font class="stringliteral">"No MPFR detected on your system.\n Please install MPFR, and compile scslib with MPFR support\n (see ./configure --help) \n"</font>); 00327 } 00328 <font class="preprocessor">#endif </font><font class="comment">/*HAVE_LIBMPFR*/</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>