/* Create tuned thresholds for various algorithms.

Copyright 1999-2003, 2005, 2006, 2008-2012 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify

it under the terms of either:

  * the GNU Lesser General Public License as published by the Free

    Software Foundation; either version 3 of the License, or (at your

    option) any later version.

or

  * the GNU General Public License as published by the Free Software

    Foundation; either version 2 of the License, or (at your option) any

    later version.

or both in parallel, as here.

The GNU MP Library is distributed in the hope that it will be useful, but

WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License

for more details.

You should have received copies of the GNU General Public License and the

GNU Lesser General Public License along with the GNU MP Library.  If not,

see https://www.gnu.org/licenses/.  */

/* Usage: tuneup [-t] [-t] [-p precision]

   -t turns on some diagnostic traces, a second -t turns on more traces.

   Notes:

   The code here isn't a vision of loveliness, mainly because it's subject

   to ongoing changes according to new things wanting to be tuned, and

   practical requirements of systems tested.

   Sometimes running the program twice produces slightly different results.

   This is probably because there's so little separating algorithms near

   their crossover, and on that basis it should make little or no difference

   to the final speed of the relevant routines, but nothing has been done to

   check that carefully.

   Algorithm:

   The thresholds are determined as follows.  A crossover may not be a

   single size but rather a range where it oscillates between method A or

   method B faster.  If the threshold is set making B used where A is faster

   (or vice versa) that's bad.  Badness is the percentage time lost and

   total badness is the sum of this over all sizes measured.  The threshold

   is set to minimize total badness.

   Suppose, as sizes increase, method B becomes faster than method A.  The

   effect of the rule is that, as you look at increasing sizes, isolated

   points where B is faster are ignored, but when it's consistently faster,

   or faster on balance, then the threshold is set there.  The same result

   is obtained thinking in the other direction of A becoming faster at

   smaller sizes.

   In practice the thresholds tend to be chosen to bring on the next

   algorithm fairly quickly.

   This rule is attractive because it's got a basis in reason and is fairly

   easy to implement, but no work has been done to actually compare it in

   absolute terms to other possibilities.

   Implementation:

   In a normal library build the thresholds are constants.  To tune them

   selected objects are recompiled with the thresholds as global variables

   instead.  #define TUNE_PROGRAM_BUILD does this, with help from code at

   the end of gmp-impl.h, and rules in tune/Makefile.am.

   MUL_TOOM22_THRESHOLD for example uses a recompiled mpn_mul_n.  The

   threshold is set to "size+1" to avoid karatsuba, or to "size" to use one

   level, but recurse into the basecase.

   MUL_TOOM33_THRESHOLD makes use of the tuned MUL_TOOM22_THRESHOLD value.

   Other routines in turn will make use of both of those.  Naturally the

   dependants must be tuned first.

   In a couple of cases, like DIVEXACT_1_THRESHOLD, there's no recompiling,

   just a threshold based on comparing two routines (mpn_divrem_1 and

   mpn_divexact_1), and no further use of the value determined.

   Flags like USE_PREINV_MOD_1 or JACOBI_BASE_METHOD are even simpler, being

   just comparisons between certain routines on representative data.

   Shortcuts are applied when native (assembler) versions of routines exist.

   For instance a native mpn_sqr_basecase is assumed to be always faster

   than mpn_mul_basecase, with no measuring.

   No attempt is made to tune within assembler routines, for instance

   DIVREM_1_NORM_THRESHOLD.  An assembler mpn_divrem_1 is expected to be

   written and tuned all by hand.  Assembler routines that might have hard

   limits are recompiled though, to make them accept a bigger range of sizes

   than normal, eg. mpn_sqr_basecase to compare against mpn_toom2_sqr.

   Limitations:

   The FFTs aren't subject to the same badness rule as the other thresholds,

   so each k is probably being brought on a touch early.  This isn't likely

   to make a difference, and the simpler probing means fewer tests.

*/

#define TUNE_PROGRAM_BUILD  1   /* for gmp-impl.h */

#include "config.h"

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

#include <time.h>

#if HAVE_UNISTD_H

#include <unistd.h>

#endif

#include "gmp.h"

#include "gmp-impl.h"

#include "longlong.h"

#include "tests.h"

#include "speed.h"

#if !HAVE_DECL_OPTARG

extern char *optarg;

extern int optind, opterr;

#endif

#define DEFAULT_MAX_SIZE   1000  /* limbs */

#if WANT_FFT

mp_size_t  option_fft_max_size = 50000;  /* limbs */

#else

mp_size_t  option_fft_max_size = 0;

#endif

int        option_trace = 0;

int        option_fft_trace = 0;

struct speed_params  s;

struct dat_t {

  mp_size_t  size;

  double     d;

} *dat = NULL;

int  ndat = 0;

int  allocdat = 0;

/* This is not defined if mpn_sqr_basecase doesn't declare a limit.  In that

   case use zero here, which for params.max_size means no limit.  */

#ifndef TUNE_SQR_TOOM2_MAX

#define TUNE_SQR_TOOM2_MAX  0

#endif

mp_size_t  mul_toom22_threshold         = MP_SIZE_T_MAX;

mp_size_t  mul_toom33_threshold         = MUL_TOOM33_THRESHOLD_LIMIT;

mp_size_t  mul_toom44_threshold         = MUL_TOOM44_THRESHOLD_LIMIT;

mp_size_t  mul_toom6h_threshold         = MUL_TOOM6H_THRESHOLD_LIMIT;

mp_size_t  mul_toom8h_threshold         = MUL_TOOM8H_THRESHOLD_LIMIT;

mp_size_t  mul_toom32_to_toom43_threshold = MP_SIZE_T_MAX;

mp_size_t  mul_toom32_to_toom53_threshold = MP_SIZE_T_MAX;

mp_size_t  mul_toom42_to_toom53_threshold = MP_SIZE_T_MAX;

mp_size_t  mul_toom42_to_toom63_threshold = MP_SIZE_T_MAX;

mp_size_t  mul_toom43_to_toom54_threshold = MP_SIZE_T_MAX;

mp_size_t  mul_fft_threshold            = MP_SIZE_T_MAX;

mp_size_t  mul_fft_modf_threshold       = MP_SIZE_T_MAX;

mp_size_t  sqr_basecase_threshold       = MP_SIZE_T_MAX;

mp_size_t  sqr_toom2_threshold

  = (TUNE_SQR_TOOM2_MAX == 0 ? MP_SIZE_T_MAX : TUNE_SQR_TOOM2_MAX);

mp_size_t  sqr_toom3_threshold          = SQR_TOOM3_THRESHOLD_LIMIT;

mp_size_t  sqr_toom4_threshold          = SQR_TOOM4_THRESHOLD_LIMIT;

mp_size_t  sqr_toom6_threshold          = SQR_TOOM6_THRESHOLD_LIMIT;

mp_size_t  sqr_toom8_threshold          = SQR_TOOM8_THRESHOLD_LIMIT;

mp_size_t  sqr_fft_threshold            = MP_SIZE_T_MAX;

mp_size_t  sqr_fft_modf_threshold       = MP_SIZE_T_MAX;

mp_size_t  mullo_basecase_threshold     = MP_SIZE_T_MAX;

mp_size_t  mullo_dc_threshold           = MP_SIZE_T_MAX;

mp_size_t  mullo_mul_n_threshold        = MP_SIZE_T_MAX;

mp_size_t  sqrlo_basecase_threshold     = MP_SIZE_T_MAX;

mp_size_t  sqrlo_dc_threshold           = MP_SIZE_T_MAX;

mp_size_t  sqrlo_sqr_threshold          = MP_SIZE_T_MAX;

mp_size_t  mulmid_toom42_threshold      = MP_SIZE_T_MAX;

mp_size_t  mulmod_bnm1_threshold        = MP_SIZE_T_MAX;

mp_size_t  sqrmod_bnm1_threshold        = MP_SIZE_T_MAX;

mp_size_t  div_qr_2_pi2_threshold       = MP_SIZE_T_MAX;

mp_size_t  dc_div_qr_threshold          = MP_SIZE_T_MAX;

mp_size_t  dc_divappr_q_threshold       = MP_SIZE_T_MAX;

mp_size_t  mu_div_qr_threshold          = MP_SIZE_T_MAX;

mp_size_t  mu_divappr_q_threshold       = MP_SIZE_T_MAX;

mp_size_t  mupi_div_qr_threshold        = MP_SIZE_T_MAX;

mp_size_t  mu_div_q_threshold           = MP_SIZE_T_MAX;

mp_size_t  dc_bdiv_qr_threshold         = MP_SIZE_T_MAX;

mp_size_t  dc_bdiv_q_threshold          = MP_SIZE_T_MAX;

mp_size_t  mu_bdiv_qr_threshold         = MP_SIZE_T_MAX;

mp_size_t  mu_bdiv_q_threshold          = MP_SIZE_T_MAX;

mp_size_t  inv_mulmod_bnm1_threshold    = MP_SIZE_T_MAX;

mp_size_t  inv_newton_threshold         = MP_SIZE_T_MAX;

mp_size_t  inv_appr_threshold           = MP_SIZE_T_MAX;

mp_size_t  binv_newton_threshold        = MP_SIZE_T_MAX;

mp_size_t  redc_1_to_redc_2_threshold   = MP_SIZE_T_MAX;

mp_size_t  redc_1_to_redc_n_threshold   = MP_SIZE_T_MAX;

mp_size_t  redc_2_to_redc_n_threshold   = MP_SIZE_T_MAX;

mp_size_t  matrix22_strassen_threshold  = MP_SIZE_T_MAX;

mp_size_t  hgcd_threshold               = MP_SIZE_T_MAX;

mp_size_t  hgcd_appr_threshold          = MP_SIZE_T_MAX;

mp_size_t  hgcd_reduce_threshold        = MP_SIZE_T_MAX;

mp_size_t  gcd_dc_threshold             = MP_SIZE_T_MAX;

mp_size_t  gcdext_dc_threshold          = MP_SIZE_T_MAX;

int	   div_qr_1n_pi1_method		= 0;

mp_size_t  div_qr_1_norm_threshold      = MP_SIZE_T_MAX;

mp_size_t  div_qr_1_unnorm_threshold    = MP_SIZE_T_MAX;

mp_size_t  divrem_1_norm_threshold      = MP_SIZE_T_MAX;

mp_size_t  divrem_1_unnorm_threshold    = MP_SIZE_T_MAX;

mp_size_t  mod_1_norm_threshold         = MP_SIZE_T_MAX;

mp_size_t  mod_1_unnorm_threshold       = MP_SIZE_T_MAX;

int	   mod_1_1p_method		= 0;

mp_size_t  mod_1n_to_mod_1_1_threshold  = MP_SIZE_T_MAX;

mp_size_t  mod_1u_to_mod_1_1_threshold  = MP_SIZE_T_MAX;

mp_size_t  mod_1_1_to_mod_1_2_threshold = MP_SIZE_T_MAX;

mp_size_t  mod_1_2_to_mod_1_4_threshold = MP_SIZE_T_MAX;

mp_size_t  preinv_mod_1_to_mod_1_threshold = MP_SIZE_T_MAX;

mp_size_t  divrem_2_threshold           = MP_SIZE_T_MAX;

mp_size_t  get_str_dc_threshold         = MP_SIZE_T_MAX;

mp_size_t  get_str_precompute_threshold = MP_SIZE_T_MAX;

mp_size_t  set_str_dc_threshold         = MP_SIZE_T_MAX;

mp_size_t  set_str_precompute_threshold = MP_SIZE_T_MAX;

mp_size_t  fac_odd_threshold            = 0;

mp_size_t  fac_dsc_threshold            = FAC_DSC_THRESHOLD_LIMIT;

mp_size_t  fft_modf_sqr_threshold = MP_SIZE_T_MAX;

mp_size_t  fft_modf_mul_threshold = MP_SIZE_T_MAX;

struct param_t {

  const char        *name;

  speed_function_t  function;

  speed_function_t  function2;

  double            step_factor;    /* how much to step relatively */

  int               step;           /* how much to step absolutely */

  double            function_fudge; /* multiplier for "function" speeds */

  int               stop_since_change;

  double            stop_factor;

  mp_size_t         min_size;

  int               min_is_always;

  mp_size_t         max_size;

  mp_size_t         check_size;

  mp_size_t         size_extra;

#define DATA_HIGH_LT_R  1

#define DATA_HIGH_GE_R  2

  int               data_high;

  int               noprint;

};

/* These are normally undefined when false, which suits "#if" fine.

   But give them zero values so they can be used in plain C "if"s.  */

#ifndef UDIV_PREINV_ALWAYS

#define UDIV_PREINV_ALWAYS 0

#endif

#ifndef HAVE_NATIVE_mpn_divexact_1

#define HAVE_NATIVE_mpn_divexact_1 0

#endif

#ifndef HAVE_NATIVE_mpn_div_qr_1n_pi1

#define HAVE_NATIVE_mpn_div_qr_1n_pi1 0

#endif

#ifndef HAVE_NATIVE_mpn_divrem_1

#define HAVE_NATIVE_mpn_divrem_1 0

#endif

#ifndef HAVE_NATIVE_mpn_divrem_2

#define HAVE_NATIVE_mpn_divrem_2 0

#endif

#ifndef HAVE_NATIVE_mpn_mod_1

#define HAVE_NATIVE_mpn_mod_1 0

#endif

#ifndef HAVE_NATIVE_mpn_mod_1_1p

#define HAVE_NATIVE_mpn_mod_1_1p 0

#endif

#ifndef HAVE_NATIVE_mpn_modexact_1_odd

#define HAVE_NATIVE_mpn_modexact_1_odd 0

#endif

#ifndef HAVE_NATIVE_mpn_preinv_divrem_1

#define HAVE_NATIVE_mpn_preinv_divrem_1 0

#endif

#ifndef HAVE_NATIVE_mpn_preinv_mod_1

#define HAVE_NATIVE_mpn_preinv_mod_1 0

#endif

#ifndef HAVE_NATIVE_mpn_sqr_basecase

#define HAVE_NATIVE_mpn_sqr_basecase 0

#endif

#define MAX3(a,b,c)  MAX (MAX (a, b), c)

mp_limb_t

randlimb_norm (void)

{

  mp_limb_t  n;

  mpn_random (&n, 1);

  n |= GMP_NUMB_HIGHBIT;

  return n;

}

#define GMP_NUMB_HALFMASK  ((CNST_LIMB(1) << (GMP_NUMB_BITS/2)) - 1)

mp_limb_t

randlimb_half (void)

{

  mp_limb_t  n;

  mpn_random (&n, 1);

  n &= GMP_NUMB_HALFMASK;

  n += (n==0);

  return n;

}

/* Add an entry to the end of the dat[] array, reallocing to make it bigger

   if necessary.  */

void

add_dat (mp_size_t size, double d)

{

#define ALLOCDAT_STEP  500

  ASSERT_ALWAYS (ndat <= allocdat);

  if (ndat == allocdat)

    {

      dat = (struct dat_t *) __gmp_allocate_or_reallocate

        (dat, allocdat * sizeof(dat[0]),

         (allocdat+ALLOCDAT_STEP) * sizeof(dat[0]));

      allocdat += ALLOCDAT_STEP;

    }

  dat[ndat].size = size;

  dat[ndat].d = d;

  ndat++;

}

/* Return the threshold size based on the data accumulated. */

mp_size_t

analyze_dat (int final)

{

  double  x, min_x;

  int     j, min_j;

  /* If the threshold is set at dat[0].size, any positive values are bad. */

  x = 0.0;

  for (j = 0; j < ndat; j++)

    if (dat[j].d > 0.0)

      x += dat[j].d;

  if (option_trace >= 2 && final)

    {

      printf ("\n");

      printf ("x is the sum of the badness from setting thresh at given size\n");

      printf ("  (minimum x is sought)\n");

      printf ("size=%ld  first x=%.4f\n", (long) dat[j].size, x);

    }

  min_x = x;

  min_j = 0;

  /* When stepping to the next dat[j].size, positive values are no longer

     bad (so subtracted), negative values become bad (so add the absolute

     value, meaning subtract). */

  for (j = 0; j < ndat; x -= dat[j].d, j++)

    {

      if (option_trace >= 2 && final)

        printf ("size=%ld  x=%.4f\n", (long) dat[j].size, x);

      if (x < min_x)

        {

          min_x = x;

          min_j = j;

        }

    }

  return min_j;

}

/* Measuring for recompiled mpn/generic/div_qr_1.c,

 * mpn/generic/divrem_1.c, mpn/generic/mod_1.c and mpz/fac_ui.c */

mp_limb_t mpn_div_qr_1_tune (mp_ptr, mp_limb_t *, mp_srcptr, mp_size_t, mp_limb_t);

#if defined (__cplusplus)

extern "C" {

#endif

mp_limb_t mpn_divrem_1_tune (mp_ptr, mp_size_t, mp_srcptr, mp_size_t, mp_limb_t);

mp_limb_t mpn_mod_1_tune (mp_srcptr, mp_size_t, mp_limb_t);

void mpz_fac_ui_tune (mpz_ptr, unsigned long);

#if defined (__cplusplus)

}

#endif

double

speed_mpn_mod_1_tune (struct speed_params *s)

{

  SPEED_ROUTINE_MPN_MOD_1 (mpn_mod_1_tune);

}

double

speed_mpn_divrem_1_tune (struct speed_params *s)

{

  SPEED_ROUTINE_MPN_DIVREM_1 (mpn_divrem_1_tune);

}

double

speed_mpz_fac_ui_tune (struct speed_params *s)

{

  SPEED_ROUTINE_MPZ_FAC_UI (mpz_fac_ui_tune);

}

double

speed_mpn_div_qr_1_tune (struct speed_params *s)

{

  SPEED_ROUTINE_MPN_DIV_QR_1 (mpn_div_qr_1_tune);

}

double

tuneup_measure (speed_function_t fun,

                const struct param_t *param,

                struct speed_params *s)

{

  static struct param_t  dummy;

  double   t;

  TMP_DECL;

  if (! param)

    param = &dummy;

  s->size += param->size_extra;

  TMP_MARK;

  SPEED_TMP_ALLOC_LIMBS (s->xp, s->size, 0);

  SPEED_TMP_ALLOC_LIMBS (s->yp, s->size, 0);

  mpn_random (s->xp, s->size);

  mpn_random (s->yp, s->size);

  switch (param->data_high) {

  case DATA_HIGH_LT_R:

    s->xp[s->size-1] %= s->r;

    s->yp[s->size-1] %= s->r;

    break;

  case DATA_HIGH_GE_R:

    s->xp[s->size-1] |= s->r;

    s->yp[s->size-1] |= s->r;

    break;

  }

  t = speed_measure (fun, s);

  s->size -= param->size_extra;

  TMP_FREE;

  return t;

}

#define PRINT_WIDTH  31

void

print_define_start (const char *name)

{

  printf ("#define %-*s  ", PRINT_WIDTH, name);

  if (option_trace)

    printf ("...\n");

}

void

print_define_end_remark (const char *name, mp_size_t value, const char *remark)

{

  if (option_trace)

    printf ("#define %-*s  ", PRINT_WIDTH, name);

  if (value == MP_SIZE_T_MAX)

    printf ("MP_SIZE_T_MAX");

  else

    printf ("%5ld", (long) value);

  if (remark != NULL)

    printf ("  /* %s */", remark);

  printf ("\n");

  fflush (stdout);

}

void

print_define_end (const char *name, mp_size_t value)

{

  const char  *remark;

  if (value == MP_SIZE_T_MAX)

    remark = "never";

  else if (value == 0)

    remark = "always";

  else

    remark = NULL;

  print_define_end_remark (name, value, remark);

}

void

print_define (const char *name, mp_size_t value)

{

  print_define_start (name);

  print_define_end (name, value);

}

void

print_define_remark (const char *name, mp_size_t value, const char *remark)

{

  print_define_start (name);

  print_define_end_remark (name, value, remark);

}

void

one (mp_size_t *threshold, struct param_t *param)

{

  int  since_positive, since_thresh_change;

  int  thresh_idx, new_thresh_idx;

#define DEFAULT(x,n)  do { if (! (x))  (x) = (n); } while (0)

  DEFAULT (param->function_fudge, 1.0);

  DEFAULT (param->function2, param->function);

  DEFAULT (param->step_factor, 0.01);  /* small steps by default */

  DEFAULT (param->step, 1);            /* small steps by default */

  DEFAULT (param->stop_since_change, 80);

  DEFAULT (param->stop_factor, 1.2);

  DEFAULT (param->min_size, 10);

  DEFAULT (param->max_size, DEFAULT_MAX_SIZE);

  if (param->check_size != 0)

    {

      double   t1, t2;

      s.size = param->check_size;

      *threshold = s.size+1;

      t1 = tuneup_measure (param->function, param, &s);

      *threshold = s.size;

      t2 = tuneup_measure (param->function2, param, &s);

      if (t1 == -1.0 || t2 == -1.0)

        {

          printf ("Oops, can't run both functions at size %ld\n",

                  (long) s.size);

          abort ();

        }

      t1 *= param->function_fudge;

      /* ask that t2 is at least 4% below t1 */

      if (t1 < t2*1.04)

        {

          if (option_trace)

            printf ("function2 never enough faster: t1=%.9f t2=%.9f\n", t1, t2);

          *threshold = MP_SIZE_T_MAX;

          if (! param->noprint)

            print_define (param->name, *threshold);

          return;

        }

      if (option_trace >= 2)

        printf ("function2 enough faster at size=%ld: t1=%.9f t2=%.9f\n",

                (long) s.size, t1, t2);

    }

  if (! param->noprint || option_trace)

    print_define_start (param->name);

  ndat = 0;

  since_positive = 0;

  since_thresh_change = 0;

  thresh_idx = 0;

  if (option_trace >= 2)

    {

      printf ("             algorithm-A  algorithm-B   ratio  possible\n");

      printf ("              (seconds)    (seconds)    diff    thresh\n");

    }

  for (s.size = param->min_size;

       s.size < param->max_size;

       s.size += MAX ((mp_size_t) floor (s.size * param->step_factor), param->step))

    {

      double   ti, tiplus1, d;

      /*

        FIXME: check minimum size requirements are met, possibly by just

        checking for the -1 returns from the speed functions.

      */

      /* using method A at this size */

      *threshold = s.size+1;

      ti = tuneup_measure (param->function, param, &s);

      if (ti == -1.0)

        abort ();

      ti *= param->function_fudge;

      /* using method B at this size */

      *threshold = s.size;

      tiplus1 = tuneup_measure (param->function2, param, &s);

      if (tiplus1 == -1.0)

        abort ();

      /* Calculate the fraction by which the one or the other routine is

         slower.  */

      if (tiplus1 >= ti)

        d = (tiplus1 - ti) / tiplus1;  /* negative */

      else

        d = (tiplus1 - ti) / ti;       /* positive */

      add_dat (s.size, d);

      new_thresh_idx = analyze_dat (0);

      if (option_trace >= 2)

        printf ("size=%ld  %.9f  %.9f  % .4f %c  %ld\n",

                (long) s.size, ti, tiplus1, d,

                ti > tiplus1 ? '#' : ' ',

                (long) dat[new_thresh_idx].size);

      /* Stop if the last time method i was faster was more than a

         certain number of measurements ago.  */

#define STOP_SINCE_POSITIVE  200

      if (d >= 0)

        since_positive = 0;

      else

        if (++since_positive > STOP_SINCE_POSITIVE)

          {

            if (option_trace >= 1)

              printf ("stopped due to since_positive (%d)\n",

                      STOP_SINCE_POSITIVE);

            break;

          }

      /* Stop if method A has become slower by a certain factor. */

      if (ti >= tiplus1 * param->stop_factor)

        {

          if (option_trace >= 1)

            printf ("stopped due to ti >= tiplus1 * factor (%.1f)\n",

                    param->stop_factor);

          break;

        }

      /* Stop if the threshold implied hasn't changed in a certain

         number of measurements.  (It's this condition that usually

         stops the loop.) */

      if (thresh_idx != new_thresh_idx)

        since_thresh_change = 0, thresh_idx = new_thresh_idx;

      else

        if (++since_thresh_change > param->stop_since_change)

          {

            if (option_trace >= 1)

              printf ("stopped due to since_thresh_change (%d)\n",

                      param->stop_since_change);

            break;

          }

      /* Stop if the threshold implied is more than a certain number of

         measurements ago.  */

#define STOP_SINCE_AFTER   500

      if (ndat - thresh_idx > STOP_SINCE_AFTER)

        {

          if (option_trace >= 1)

            printf ("stopped due to ndat - thresh_idx > amount (%d)\n",

                    STOP_SINCE_AFTER);

          break;

        }

      /* Stop when the size limit is reached before the end of the

         crossover, but only show this as an error for >= the default max

         size.  FIXME: Maybe should make it a param choice whether this is

         an error.  */

      if (s.size >= param->max_size && param->max_size >= DEFAULT_MAX_SIZE)

        {

          fprintf (stderr, "%s\n", param->name);

          fprintf (stderr, "sizes %ld to %ld total %d measurements\n",

                   (long) dat[0].size, (long) dat[ndat-1].size, ndat);

          fprintf (stderr, "    max size reached before end of crossover\n");

          break;

        }

    }

  if (option_trace >= 1)

    printf ("sizes %ld to %ld total %d measurements\n",

            (long) dat[0].size, (long) dat[ndat-1].size, ndat);

  *threshold = dat[analyze_dat (1)].size;

  if (param->min_is_always)

    {

      if (*threshold == param->min_size)

        *threshold = 0;

    }

  if (! param->noprint || option_trace)

    print_define_end (param->name, *threshold);

}

/* Special probing for the fft thresholds.  The size restrictions on the

   FFTs mean the graph of time vs size has a step effect.  See this for

   example using

       ./speed -s 4096-16384 -t 128 -P foo mpn_mul_fft.8 mpn_mul_fft.9

       gnuplot foo.gnuplot

   The current approach is to compare routines at the midpoint of relevant

   steps.  Arguably a more sophisticated system of threshold data is wanted

   if this step effect remains. */

struct fft_param_t {

  const char        *table_name;

  const char        *threshold_name;

  const char        *modf_threshold_name;

  mp_size_t         *p_threshold;

  mp_size_t         *p_modf_threshold;

  mp_size_t         first_size;

  mp_size_t         max_size;

  speed_function_t  function;

  speed_function_t  mul_modf_function;

  speed_function_t  mul_function;

  mp_size_t         sqr;

};

/* mpn_mul_fft requires pl a multiple of 2^k limbs, but with

   N=pl*BIT_PER_MP_LIMB it internally also pads out so N/2^k is a multiple

   of 2^(k-1) bits. */

mp_size_t

fft_step_size (int k)

{

  mp_size_t  step;

  step = MAX ((mp_size_t) 1 << (k-1), GMP_LIMB_BITS) / GMP_LIMB_BITS;

  step *= (mp_size_t) 1 << k;

  if (step <= 0)

    {

      printf ("Can't handle k=%d\n", k);

      abort ();

    }

  return step;

}

mp_size_t

fft_next_size (mp_size_t pl, int k)

{

  mp_size_t  m = fft_step_size (k);

/*    printf ("[k=%d %ld] %ld ->", k, m, pl); */

  if (pl == 0 || (pl & (m-1)) != 0)

    pl = (pl | (m-1)) + 1;

/*    printf (" %ld\n", pl); */

  return pl;

}

#define NMAX_DEFAULT 1000000

#define MAX_REPS 25

#define MIN_REPS 5

static inline size_t

mpn_mul_fft_lcm (size_t a, unsigned int k)

{

  unsigned int l = k;

  while (a % 2 == 0 && k > 0)

    {

      a >>= 1;

      k--;

    }

  return a << l;

}

mp_size_t

fftfill (mp_size_t pl, int k, int sqr)

{

  mp_size_t maxLK;

  mp_bitcnt_t N, Nprime, nprime, M;

  N = pl * GMP_NUMB_BITS;

  M = N >> k;

  maxLK = mpn_mul_fft_lcm ((unsigned long) GMP_NUMB_BITS, k);

  Nprime = (1 + (2 * M + k + 2) / maxLK) * maxLK;

  nprime = Nprime / GMP_NUMB_BITS;

  if (nprime >= (sqr ? SQR_FFT_MODF_THRESHOLD : MUL_FFT_MODF_THRESHOLD))

    {

      size_t K2;

      for (;;)

	{

	  K2 = 1L << mpn_fft_best_k (nprime, sqr);

	  if ((nprime & (K2 - 1)) == 0)

	    break;

	  nprime = (nprime + K2 - 1) & -K2;

	  Nprime = nprime * GMP_LIMB_BITS;

	}

    }

  ASSERT_ALWAYS (nprime < pl);

  return Nprime;

}

static int

compare_double (const void *ap, const void *bp)

{

  double a = * (const double *) ap;

  double b = * (const double *) bp;

  if (a < b)

    return -1;

  else if (a > b)

    return 1;

  else

    return 0;

}

double

median (double *times, int n)

{

  qsort (times, n, sizeof (double), compare_double);

  return times[n/2];

}

#define FFT_CACHE_SIZE 25

typedef struct fft_cache

{

  mp_size_t n;

  double time;

} fft_cache_t;

fft_cache_t fft_cache[FFT_CACHE_SIZE];

double

cached_measure (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n, int k,

		int n_measurements)

{

  int i;

  double t, ttab[MAX_REPS];

  if (fft_cache[k].n == n)

    return fft_cache[k].time;

  for (i = 0; i < n_measurements; i++)

    {

      speed_starttime ();

      mpn_mul_fft (rp, n, ap, n, bp, n, k);

      ttab[i] = speed_endtime ();

    }

  t = median (ttab, n_measurements);

  fft_cache[k].n = n;

  fft_cache[k].time = t;

  return t;

}

#define INSERT_FFTTAB(idx, nval, kval)					\

  do {									\

    fft_tab[idx].n = nval;						\

    fft_tab[idx].k = kval;						\

    fft_tab[idx+1].n = (1 << 27) - 1;	/* sentinel, 27b wide field */	\

    fft_tab[idx+1].k = (1 <<  5) - 1;					\

  } while (0)

int

fftmes (mp_size_t nmin, mp_size_t nmax, int initial_k, struct fft_param_t *p, int idx, int print)

{

  mp_size_t n, n1, prev_n1;