00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
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
00047
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
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);
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
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
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
00248
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
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