Main Page   Data Structures   File List   Data Fields   Globals  

test_accuracy.c

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

Generated on Tue Jun 17 10:15:51 2003 for SCSLib by doxygen1.2.15