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