contrib/mpfr/tuneup.c

   1 /* Tune various threshold of MPFR
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 2.1 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #include <stdlib.h>
  24 #include <time.h>
  25
  26 #define MPFR_NEED_LONGLONG_H
  27 #include "mpfr-impl.h"
  28
  29 #undef _PROTO
  30 #define _PROTO __GMP_PROTO
  31 #include "speed.h"
  32
  33 int verbose;
  34
  35 /* s->size: precision of both input and output
  36    s->xp  : Mantissa of first input
  37    s->yp  : mantissa of second input                    */
  38
  39 #define SPEED_MPFR_FUNC(mean_fun) do {               \
  40   unsigned  i;                                       \
  41   mp_ptr    wp;                                      \
  42   double    t;                                       \
  43   mpfr_t    w, x;                                    \
  44   mp_size_t size;                                    \
  45   MPFR_TMP_DECL (marker);                            \
  46                                                      \
  47   SPEED_RESTRICT_COND (s->size >= MPFR_PREC_MIN);    \
  48   SPEED_RESTRICT_COND (s->size <= MPFR_PREC_MAX);    \
  49   MPFR_TMP_MARK (marker);                            \
  50                                                      \
  51   size = (s->size-1)/BITS_PER_MP_LIMB+1;             \
  52   s->xp[size-1] |= MPFR_LIMB_HIGHBIT;                \
  53   MPFR_TMP_INIT1 (s->xp, x, s->size);                \
  54   MPFR_SET_EXP (x, 0);                               \
  55                                                      \
  56   MPFR_TMP_INIT (wp, w, s->size, size);              \
  57                                                      \
  58   speed_operand_src (s, s->xp, size);                \
  59   speed_operand_dst (s, wp, size);                   \
  60   speed_cache_fill (s);                              \
  61                                                      \
  62   speed_starttime ();                                \
  63   i = s->reps;                                       \
  64   do                                                 \
  65     mean_fun (w, x, GMP_RNDN);                       \
  66   while (--i != 0);                                  \
  67   t = speed_endtime ();                              \
  68                                                      \
  69   MPFR_TMP_FREE (marker);                            \
  70   return t;                                          \
  71 } while (0)
  72
  73 #define SPEED_MPFR_OP(mean_fun) do {                 \
  74   unsigned  i;                                       \
  75   mp_ptr    wp;                                      \
  76   double    t;                                       \
  77   mpfr_t    w, x, y;                                 \
  78   mp_size_t size;                                    \
  79   MPFR_TMP_DECL (marker);                            \
  80                                                      \
  81   SPEED_RESTRICT_COND (s->size >= MPFR_PREC_MIN);    \
  82   SPEED_RESTRICT_COND (s->size <= MPFR_PREC_MAX);    \
  83   MPFR_TMP_MARK (marker);                            \
  84                                                      \
  85   size = (s->size-1)/BITS_PER_MP_LIMB+1;             \
  86   s->xp[size-1] |= MPFR_LIMB_HIGHBIT;                \
  87   MPFR_TMP_INIT1 (s->xp, x, s->size);                \
  88   MPFR_SET_EXP (x, 0);                               \
  89   s->yp[size-1] |= MPFR_LIMB_HIGHBIT;                \
  90   MPFR_TMP_INIT1 (s->yp, y, s->size);                \
  91   MPFR_SET_EXP (y, 0);                               \
  92                                                      \
  93   MPFR_TMP_INIT (wp, w, s->size, size);              \
  94                                                      \
  95   speed_operand_src (s, s->xp, size);                \
  96   speed_operand_src (s, s->yp, size);                \
  97   speed_operand_dst (s, wp, size);                   \
  98   speed_cache_fill (s);                              \
  99                                                      \
 100   speed_starttime ();                                \
 101   i = s->reps;                                       \
 102   do                                                 \
 103     mean_fun (w, x, y, GMP_RNDN);                    \
 104   while (--i != 0);                                  \
 105   t = speed_endtime ();                              \
 106                                                      \
 107   MPFR_TMP_FREE (marker);                            \
 108   return t;                                          \
 109 } while (0)
 110
 111
 112 /* First we include all the functions we want to tune inside this program.
 113    We can't use GNU MPFR library since the THRESHOLD can't vary */
 114
 115 /* Setup mpfr_exp_2 */
 116 mp_prec_t mpfr_exp_2_threshold;
 117 #undef  MPFR_EXP_2_THRESHOLD
 118 #define MPFR_EXP_2_THRESHOLD mpfr_exp_2_threshold
 119 #include "exp_2.c"
 120 static double speed_mpfr_exp_2 (struct speed_params *s) {
 121   SPEED_MPFR_FUNC (mpfr_exp_2);
 122 }
 123
 124 /* Setup mpfr_exp */
 125 mp_prec_t mpfr_exp_threshold;
 126 #undef  MPFR_EXP_THRESHOLD
 127 #define MPFR_EXP_THRESHOLD mpfr_exp_threshold
 128 #include "exp.c"
 129 static double speed_mpfr_exp (struct speed_params *s) {
 130   SPEED_MPFR_FUNC (mpfr_exp);
 131 }
 132
 133 /* Setup mpfr_mul */
 134 mp_prec_t mpfr_mul_threshold;
 135 #undef  MPFR_MUL_THRESHOLD
 136 #define MPFR_MUL_THRESHOLD mpfr_mul_threshold
 137 #include "mul.c"
 138 static double speed_mpfr_mul (struct speed_params *s) {
 139   SPEED_MPFR_OP (mpfr_mul);
 140 }
 141
 142
 143
 144 /************************************************
 145  * Common functions (inspired by GMP function)  *
 146  ************************************************/
 147 static int
 148 analyze_data (double *dat, int ndat) {
 149   double  x, min_x;
 150   int     j, min_j;
 151
 152   x = 0.0;
 153   for (j = 0; j < ndat; j++)
 154     if (dat[j] > 0.0)
 155       x += dat[j];
 156
 157   min_x = x;
 158   min_j = 0;
 159
 160   for (j = 0; j < ndat; x -= dat[j], j++)
 161     {
 162       if (x < min_x)
 163         {
 164           min_x = x;
 165           min_j = j;
 166         }
 167     }
 168   return min_j;
 169 }
 170
 171 #define THRESHOLD_WINDOW 16
 172 #define THRESHOLD_FINAL_WINDOW 128
 173 static double domeasure (mp_prec_t *threshold,
 174                          double (*func) (struct speed_params *),
 175                          mp_prec_t p)
 176 {
 177   struct speed_params s;
 178   mp_size_t size;
 179   double t1, t2, d;
 180
 181   s.align_xp = s.align_yp = s.align_wp = 64;
 182   s.size = p;
 183   size = (p - 1)/BITS_PER_MP_LIMB+1;
 184   s.xp = malloc (2*size*sizeof (mp_limb_t));
 185   if (s.xp == NULL)
 186     {
 187       fprintf (stderr, "Can't allocate memory.\n");
 188       abort ();
 189     }
 190   mpn_random (s.xp, size);
 191   s.yp = s.xp + size;
 192   mpn_random (s.yp, size);
 193   *threshold = MPFR_PREC_MAX;
 194   t1 = speed_measure (func, &s);
 195   if (t1 == -1.0)
 196     {
 197       fprintf (stderr, "Failed to measure function 1!\n");
 198       abort ();
 199     }
 200   *threshold = 1;
 201   t2 = speed_measure (func, &s);
 202   if (t2 == -1.0)
 203     {
 204       fprintf (stderr, "Failed to measure function 2!\n");
 205       abort ();
 206     }
 207   free (s.xp);
 208   /* t1 is the time of the first algo (used for low prec) */
 209   if (t2 >= t1)
 210     d = (t2 - t1) / t2;
 211   else
 212     d = (t2 - t1) / t1;
 213   /* d > 0 if we have to use algo 1.
 214      d < 0 if we have to use algo 2 */
 215   return d;
 216 }
 217
 218 /* Tune a function with a simple THRESHOLD
 219    The function doesn't depend on another threshold.
 220    It assumes that it uses algo1 if p < THRESHOLD
 221    and algo2 otherwise.
 222    if algo2 is better for low prec, and algo1 better for high prec,
 223    the behaviour of this function is undefined. */
 224 static void
 225 tune_simple_func (mp_prec_t *threshold,
 226                   double (*func) (struct speed_params *),
 227                   mp_prec_t pstart)
 228 {
 229   double measure[THRESHOLD_FINAL_WINDOW+1];
 230   double d;
 231   mp_prec_t pstep;
 232   int i, numpos, numneg, try;
 233   mp_prec_t pmin, pmax, p;
 234
 235   /* first look for a lower bound within 10% */
 236   pmin = p = pstart;
 237   d = domeasure (threshold, func, pmin);
 238   if (d < 0.0) {
 239     if (verbose)
 240       printf ("Oops: even for %lu, algo 2 seems to be faster!\n",
 241               (unsigned long) pmin);
 242     *threshold = MPFR_PREC_MIN;
 243     return;
 244   }
 245   if (d >= 1.00)
 246     for (;;) {
 247       d = domeasure (threshold, func, pmin);
 248       if (d < 1.00)
 249         break;
 250       p = pmin;
 251       pmin += pmin/2;
 252     }
 253   pmin = p;
 254   for (;;) {
 255     d = domeasure (threshold, func, pmin);
 256     if (d < 0.10)
 257       break;
 258     pmin += BITS_PER_MP_LIMB;
 259   }
 260
 261   /* then look for a upper bound within 20% */
 262   pmax = pmin * 2;
 263   for (;;) {
 264     d = domeasure (threshold, func, pmax);
 265     if (d < -0.20)
 266       break;
 267     pmax += pmin / 2; /* don't increase too rapidly */
 268   }
 269
 270   /* The threshold is between pmin and pmax. Affine them */
 271   try = 0;
 272   while ((pmax-pmin) >= THRESHOLD_FINAL_WINDOW)
 273     {
 274       pstep = MAX(MIN(BITS_PER_MP_LIMB/2,(pmax-pmin)/(2*THRESHOLD_WINDOW)),1);
 275       if (verbose)
 276         printf ("Pmin = %8lu Pmax = %8lu Pstep=%lu\n", pmin, pmax, pstep);
 277       p = (pmin + pmax) / 2;
 278       for (i = numpos = numneg = 0 ; i < THRESHOLD_WINDOW + 1 ; i++) {
 279         measure[i] = domeasure (threshold, func,
 280                                 p+(i-THRESHOLD_WINDOW/2)*pstep);
 281         if (measure[i] > 0)
 282           numpos ++;
 283         else if (measure[i] < 0)
 284           numneg ++;
 285       }
 286       if (numpos > numneg)
 287         /* We use more often algo 1 than algo 2 */
 288         pmin = p - THRESHOLD_WINDOW/2*pstep;
 289       else if (numpos < numneg)
 290         pmax = p + THRESHOLD_WINDOW/2*pstep;
 291       else
 292         /* numpos == numneg ... */
 293         if (++ try > 2) {
 294           *threshold = p;
 295           if (verbose)
 296             printf ("Quick find: %lu\n", *threshold);
 297           return ;
 298         }
 299     }
 300
 301   /* Final tune... */
 302   if (verbose)
 303     printf ("Finalizing in [%lu, %lu]... ", pmin, pmax);
 304   for (i = 0 ; i < THRESHOLD_FINAL_WINDOW+1 ; i++)
 305     measure[i] = domeasure (threshold, func, pmin+i);
 306   i = analyze_data (measure, THRESHOLD_FINAL_WINDOW+1);
 307   *threshold = pmin + i;
 308   if (verbose)
 309     printf ("%lu\n", *threshold);
 310   return;
 311 }
 312
 313
 314
 315 /************************************
 316  * Tune Mulders' mulhigh function   *
 317  ************************************/
 318 #define TOLERANCE 1.02
 319 #ifndef MPFR_MULHIGH_SIZE
 320 # define MPFR_MULHIGH_SIZE 1024
 321 #endif
 322 #ifndef MPFR_SQRHIGH_SIZE
 323 # define MPFR_SQRHIGH_SIZE 1024
 324 #endif
 325 #define MPFR_MULHIGH_TAB_SIZE MPFR_MULHIGH_SIZE
 326 #define MPFR_SQRHIGH_TAB_SIZE MPFR_SQRHIGH_SIZE
 327 #include "mulders.c"
 328
 329 static double speed_mpfr_mulhigh (struct speed_params *s) {
 330   SPEED_ROUTINE_MPN_MUL_N (mpfr_mulhigh_n);
 331 }
 332 static double speed_mpfr_sqrhigh (struct speed_params *s) {
 333   SPEED_ROUTINE_MPN_SQR (mpfr_sqrhigh_n);
 334 }
 335
 336 #define MAX_STEPS 32 /* maximum number of values of k tried for a given n */
 337
 338 /* Tune size N */
 339 static mp_size_t
 340 tune_mulders_upto (mp_size_t n)
 341 {
 342   struct speed_params s;
 343   mp_size_t k, kbest, step;
 344   double t, tbest;
 345   MPFR_TMP_DECL (marker);
 346
 347   if (n == 0)
 348     return -1;
 349
 350   MPFR_TMP_MARK (marker);
 351   s.align_xp = s.align_yp = s.align_wp = 64;
 352   s.size = n;
 353   s.xp   = MPFR_TMP_ALLOC (n*sizeof (mp_limb_t));
 354   s.yp   = MPFR_TMP_ALLOC (n*sizeof (mp_limb_t));
 355   mpn_random (s.xp, n);
 356   mpn_random (s.yp, n);
 357
 358   /* Check k == -1, mpn_mul_basecase */
 359   mulhigh_ktab[n] = -1;
 360   kbest = -1;
 361   tbest =  speed_measure (speed_mpfr_mulhigh, &s);
 362
 363   /* Check k == 0, mpn_mulhigh_n_basecase */
 364   mulhigh_ktab[n] = 0;
 365   t = speed_measure (speed_mpfr_mulhigh, &s);
 366   if (t * TOLERANCE < tbest)
 367     kbest = 0, tbest = t;
 368
 369   /* Check Mulders with cutoff point k */
 370   step = 1 + n / (2 * MAX_STEPS);
 371   for (k = n / 2 + 1 ; k < n ; k += step)
 372     {
 373       mulhigh_ktab[n] = k;
 374       t =  speed_measure (speed_mpfr_mulhigh, &s);
 375       if (t * TOLERANCE < tbest)
 376         kbest = k, tbest = t;
 377     }
 378
 379   mulhigh_ktab[n] = kbest;
 380
 381   MPFR_TMP_FREE (marker);
 382   return kbest;
 383 }
 384
 385 /* Tune size N */
 386 static mp_size_t
 387 tune_sqr_mulders_upto (mp_size_t n)
 388 {
 389   struct speed_params s;
 390   mp_size_t k, kbest, step;
 391   double t, tbest;
 392   MPFR_TMP_DECL (marker);
 393
 394   if (n == 0)
 395     return -1;
 396
 397   MPFR_TMP_MARK (marker);
 398   s.align_xp = s.align_wp = 64;
 399   s.size = n;
 400   s.xp   = MPFR_TMP_ALLOC (n*sizeof (mp_limb_t));
 401   mpn_random (s.xp, n);
 402
 403   /* Check k == -1, mpn_sqr_basecase */
 404   sqrhigh_ktab[n] = -1;
 405   kbest = -1;
 406   tbest =  speed_measure (speed_mpfr_sqrhigh, &s);
 407
 408   /* Check k == 0, mpfr_mulhigh_n_basecase */
 409   sqrhigh_ktab[n] = 0;
 410   t = speed_measure (speed_mpfr_sqrhigh, &s);
 411   if (t * TOLERANCE < tbest)
 412     kbest = 0, tbest = t;
 413
 414   /* Check Mulders */
 415   step = 1 + n / (2 * MAX_STEPS);
 416   for (k = n / 2 + 1 ; k < n ; k += step)
 417     {
 418       sqrhigh_ktab[n] = k;
 419       t =  speed_measure (speed_mpfr_sqrhigh, &s);
 420       if (t * TOLERANCE < tbest)
 421         kbest = k, tbest = t;
 422     }
 423
 424   sqrhigh_ktab[n] = kbest;
 425
 426   MPFR_TMP_FREE (marker);
 427   return kbest;
 428 }
 429
 430 static void
 431 tune_mulders (FILE *f)
 432 {
 433   mp_size_t k;
 434
 435   if (verbose)
 436     printf ("Tuning mpfr_mulhigh_n[%d]", (int) MPFR_MULHIGH_TAB_SIZE);
 437   fprintf (f, "#define MPFR_MULHIGH_TAB  \\\n ");
 438   for (k = 0 ; k < MPFR_MULHIGH_TAB_SIZE ; k++)
 439     {
 440       fprintf (f, "%d", (int) tune_mulders_upto (k));
 441       if (k != MPFR_MULHIGH_TAB_SIZE-1)
 442         fputc (',', f);
 443       if ((k+1) % 16 == 0)
 444         fprintf (f, " \\\n ");
 445       if (verbose)
 446         putchar ('.');
 447     }
 448   fprintf (f, " \n");
 449   if (verbose)
 450     putchar ('\n');
 451 }
 452
 453 static void
 454 tune_sqr_mulders (FILE *f)
 455 {
 456   mp_size_t k;
 457
 458   if (verbose)
 459     printf ("Tuning mpfr_sqrhigh_n[%d]", (int) MPFR_SQRHIGH_TAB_SIZE);
 460   fprintf (f, "#define MPFR_SQRHIGH_TAB  \\\n ");
 461   for (k = 0 ; k < MPFR_SQRHIGH_TAB_SIZE ; k++)
 462     {
 463       fprintf (f, "%d", (int) tune_sqr_mulders_upto (k));
 464       if (k != MPFR_SQRHIGH_TAB_SIZE-1)
 465         fputc (',', f);
 466       if ((k+1) % 16 == 0)
 467         fprintf (f, " \\\n ");
 468       if (verbose)
 469         putchar ('.');
 470     }
 471   fprintf (f, " \n");
 472   if (verbose)
 473     putchar ('\n');
 474 }
 475
 476 /*******************************************************
 477  *            Tune all the threshold of MPFR           *
 478  * Warning: tune the function in their dependent order!*
 479  *******************************************************/
 480 static void
 481 all (const char *filename)
 482 {
 483   FILE *f;
 484   time_t  start_time, end_time;
 485   struct tm  *tp;
 486
 487   f = fopen (filename, "w");
 488   if (f == NULL) {
 489     fprintf (stderr, "Can't open file '%s' for writing.\n", filename);
 490     abort ();
 491   }
 492
 493   speed_time_init ();
 494   if (verbose) {
 495     printf ("Using: %s\n", speed_time_string);
 496     printf ("speed_precision %d", speed_precision);
 497     if (speed_unittime == 1.0)
 498       printf (", speed_unittime 1 cycle");
 499     else
 500       printf (", speed_unittime %.2e secs", speed_unittime);
 501     if (speed_cycletime == 1.0 || speed_cycletime == 0.0)
 502       printf (", CPU freq unknown\n");
 503     else
 504       printf (", CPU freq %.2f MHz\n\n", 1e-6/speed_cycletime);
 505   }
 506
 507   time (&start_time);
 508   tp = localtime (&start_time);
 509   fprintf (f, "/* Generated by MPFR's tuneup.c, %d-%02d-%02d, ",
 510           tp->tm_year+1900, tp->tm_mon+1, tp->tm_mday);
 511
 512 #ifdef __ICC
 513   fprintf (f, "icc %d.%d.%d */\n", __ICC / 100, __ICC / 10 % 10, __ICC % 10);
 514 #elif defined(__GNUC__)
 515 #ifdef __GNUC_PATCHLEVEL__
 516   fprintf (f, "gcc %d.%d.%d */\n", __GNUC__, __GNUC_MINOR__,
 517            __GNUC_PATCHLEVEL__);
 518 #else
 519   fprintf (f, "gcc %d.%d */\n", __GNUC__, __GNUC_MINOR__);
 520 #endif
 521 #elif defined (__SUNPRO_C)
 522   fprintf (f, "Sun C %d.%d */\n", __SUNPRO_C / 0x100, __SUNPRO_C % 0x100);
 523 #elif defined (__sgi) && defined (_COMPILER_VERSION)
 524   fprintf (f, "MIPSpro C %d.%d.%d */\n",
 525            _COMPILER_VERSION / 100,
 526            _COMPILER_VERSION / 10 % 10,
 527            _COMPILER_VERSION % 10);
 528 #elif defined (__DECC) && defined (__DECC_VER)
 529   fprintf (f, "DEC C %d */\n", __DECC_VER);
 530 #else
 531   fprintf (f, "system compiler */\n");
 532 #endif
 533   fprintf (f, "\n");
 534
 535   /* Tune mulhigh */
 536   tune_mulders (f);
 537
 538   /* Tune sqrhigh */
 539   tune_sqr_mulders (f);
 540
 541   /* Tune mpfr_mul (threshold is in limbs, but it doesn't matter too much) */
 542   if (verbose)
 543     printf ("Tuning mpfr_mul...\n");
 544   tune_simple_func (&mpfr_mul_threshold, speed_mpfr_mul,
 545                     2*BITS_PER_MP_LIMB+1);
 546   fprintf (f, "#define MPFR_MUL_THRESHOLD %lu /* limbs */\n",
 547            (unsigned long) (mpfr_mul_threshold - 1) / BITS_PER_MP_LIMB + 1);
 548
 549   /* Tune mpfr_exp_2 */
 550   if (verbose)
 551     printf ("Tuning mpfr_exp_2...\n");
 552   tune_simple_func (&mpfr_exp_2_threshold, speed_mpfr_exp_2,
 553                     MPFR_PREC_MIN);
 554   mpfr_exp_2_threshold = MAX (BITS_PER_MP_LIMB, mpfr_exp_2_threshold);
 555   fprintf (f, "#define MPFR_EXP_2_THRESHOLD %lu /* bits */\n",
 556            (unsigned long) mpfr_exp_2_threshold);
 557
 558   /* Tune mpfr_exp */
 559   if (verbose)
 560     printf ("Tuning mpfr_exp...\n");
 561   tune_simple_func (&mpfr_exp_threshold, speed_mpfr_exp,
 562                     MPFR_PREC_MIN+3*BITS_PER_MP_LIMB);
 563   fprintf (f, "#define MPFR_EXP_THRESHOLD %lu /* bits */\n",
 564            (unsigned long) mpfr_exp_threshold);
 565
 566   /* End of tuning */
 567   time (&end_time);
 568   fprintf (f, "/* Tuneup completed successfully, took %ld seconds */\n",
 569            end_time - start_time);
 570   if (verbose)
 571     printf ("Complete (took %ld seconds).\n", end_time - start_time);
 572
 573   fclose (f);
 574 }
 575
 576
 577 /* Main function */
 578 int main (int argc, char *argv[])
 579 {
 580   /* Unbuffered so if output is redirected to a file it isn't lost if the
 581      program is killed part way through.  */
 582   setbuf (stdout, NULL);
 583   setbuf (stderr, NULL);
 584
 585   verbose = argc > 1;
 586
 587   if (verbose)
 588     printf ("Tuning MPFR (Coffee time?)...\n");
 589
 590   all ("mparam.h");
 591
 592   return 0;
 593 }