tuneup.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:65k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* Create tuned thresholds for various algorithms.
  2. Copyright 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2008, 2009, 2010 Free
  3. Software Foundation, Inc.
  4. This file is part of the GNU MP Library.
  5. The GNU MP Library is free software; you can redistribute it and/or modify
  6. it under the terms of the GNU Lesser General Public License as published by
  7. the Free Software Foundation; either version 3 of the License, or (at your
  8. option) any later version.
  9. The GNU MP Library is distributed in the hope that it will be useful, but
  10. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  11. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  12. License for more details.
  13. You should have received a copy of the GNU Lesser General Public License
  14. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  15. /* Usage: tuneup [-t] [-t] [-p precision]
  16.    -t turns on some diagnostic traces, a second -t turns on more traces.
  17.    Notes:
  18.    The code here isn't a vision of loveliness, mainly because it's subject
  19.    to ongoing changes according to new things wanting to be tuned, and
  20.    practical requirements of systems tested.
  21.    Sometimes running the program twice produces slightly different results.
  22.    This is probably because there's so little separating algorithms near
  23.    their crossover, and on that basis it should make little or no difference
  24.    to the final speed of the relevant routines, but nothing has been done to
  25.    check that carefully.
  26.    Algorithm:
  27.    The thresholds are determined as follows.  A crossover may not be a
  28.    single size but rather a range where it oscillates between method A or
  29.    method B faster.  If the threshold is set making B used where A is faster
  30.    (or vice versa) that's bad.  Badness is the percentage time lost and
  31.    total badness is the sum of this over all sizes measured.  The threshold
  32.    is set to minimize total badness.
  33.    Suppose, as sizes increase, method B becomes faster than method A.  The
  34.    effect of the rule is that, as you look at increasing sizes, isolated
  35.    points where B is faster are ignored, but when it's consistently faster,
  36.    or faster on balance, then the threshold is set there.  The same result
  37.    is obtained thinking in the other direction of A becoming faster at
  38.    smaller sizes.
  39.    In practice the thresholds tend to be chosen to bring on the next
  40.    algorithm fairly quickly.
  41.    This rule is attractive because it's got a basis in reason and is fairly
  42.    easy to implement, but no work has been done to actually compare it in
  43.    absolute terms to other possibilities.
  44.    Implementation:
  45.    In a normal library build the thresholds are constants.  To tune them
  46.    selected objects are recompiled with the thresholds as global variables
  47.    instead.  #define TUNE_PROGRAM_BUILD does this, with help from code at
  48.    the end of gmp-impl.h, and rules in tune/Makefile.am.
  49.    MUL_TOOM22_THRESHOLD for example uses a recompiled mpn_mul_n.  The
  50.    threshold is set to "size+1" to avoid karatsuba, or to "size" to use one
  51.    level, but recurse into the basecase.
  52.    MUL_TOOM33_THRESHOLD makes use of the tuned MUL_TOOM22_THRESHOLD value.
  53.    Other routines in turn will make use of both of those.  Naturally the
  54.    dependants must be tuned first.
  55.    In a couple of cases, like DIVEXACT_1_THRESHOLD, there's no recompiling,
  56.    just a threshold based on comparing two routines (mpn_divrem_1 and
  57.    mpn_divexact_1), and no further use of the value determined.
  58.    Flags like USE_PREINV_MOD_1 or JACOBI_BASE_METHOD are even simpler, being
  59.    just comparisons between certain routines on representative data.
  60.    Shortcuts are applied when native (assembler) versions of routines exist.
  61.    For instance a native mpn_sqr_basecase is assumed to be always faster
  62.    than mpn_mul_basecase, with no measuring.
  63.    No attempt is made to tune within assembler routines, for instance
  64.    DIVREM_1_NORM_THRESHOLD.  An assembler mpn_divrem_1 is expected to be
  65.    written and tuned all by hand.  Assembler routines that might have hard
  66.    limits are recompiled though, to make them accept a bigger range of sizes
  67.    than normal, eg. mpn_sqr_basecase to compare against mpn_toom2_sqr.
  68.    Limitations:
  69.    The FFTs aren't subject to the same badness rule as the other thresholds,
  70.    so each k is probably being brought on a touch early.  This isn't likely
  71.    to make a difference, and the simpler probing means fewer tests.
  72. */
  73. #define TUNE_PROGRAM_BUILD  1   /* for gmp-impl.h */
  74. #include "config.h"
  75. #include <math.h>
  76. #include <stdio.h>
  77. #include <stdlib.h>
  78. #include <time.h>
  79. #if HAVE_UNISTD_H
  80. #include <unistd.h>
  81. #endif
  82. #include "gmp.h"
  83. #include "gmp-impl.h"
  84. #include "longlong.h"
  85. #include "tests.h"
  86. #include "speed.h"
  87. #if !HAVE_DECL_OPTARG
  88. extern char *optarg;
  89. extern int optind, opterr;
  90. #endif
  91. #define DEFAULT_MAX_SIZE   1000  /* limbs */
  92. #if WANT_FFT
  93. mp_size_t  option_fft_max_size = 50000;  /* limbs */
  94. #else
  95. mp_size_t  option_fft_max_size = 0;
  96. #endif
  97. int        option_trace = 0;
  98. int        option_fft_trace = 0;
  99. struct speed_params  s;
  100. struct dat_t {
  101.   mp_size_t  size;
  102.   double     d;
  103. } *dat = NULL;
  104. int  ndat = 0;
  105. int  allocdat = 0;
  106. /* This is not defined if mpn_sqr_basecase doesn't declare a limit.  In that
  107.    case use zero here, which for params.max_size means no limit.  */
  108. #ifndef TUNE_SQR_TOOM2_MAX
  109. #define TUNE_SQR_TOOM2_MAX  0
  110. #endif
  111. mp_size_t  mul_toom22_threshold         = MP_SIZE_T_MAX;
  112. mp_size_t  mul_toom33_threshold         = MUL_TOOM33_THRESHOLD_LIMIT;
  113. mp_size_t  mul_toom44_threshold         = MUL_TOOM44_THRESHOLD_LIMIT;
  114. mp_size_t  mul_toom6h_threshold         = MUL_TOOM6H_THRESHOLD_LIMIT;
  115. mp_size_t  mul_toom8h_threshold         = MUL_TOOM8H_THRESHOLD_LIMIT;
  116. mp_size_t  mul_toom32_to_toom43_threshold = MP_SIZE_T_MAX;
  117. mp_size_t  mul_toom32_to_toom53_threshold = MP_SIZE_T_MAX;
  118. mp_size_t  mul_toom42_to_toom53_threshold = MP_SIZE_T_MAX;
  119. mp_size_t  mul_toom42_to_toom63_threshold = MP_SIZE_T_MAX;
  120. mp_size_t  mul_fft_threshold            = MP_SIZE_T_MAX;
  121. mp_size_t  mul_fft_modf_threshold       = MP_SIZE_T_MAX;
  122. mp_size_t  sqr_basecase_threshold       = MP_SIZE_T_MAX;
  123. mp_size_t  sqr_toom2_threshold
  124.   = (TUNE_SQR_TOOM2_MAX == 0 ? MP_SIZE_T_MAX : TUNE_SQR_TOOM2_MAX);
  125. mp_size_t  sqr_toom3_threshold          = SQR_TOOM3_THRESHOLD_LIMIT;
  126. mp_size_t  sqr_toom4_threshold          = SQR_TOOM4_THRESHOLD_LIMIT;
  127. mp_size_t  sqr_toom6_threshold          = SQR_TOOM6_THRESHOLD_LIMIT;
  128. mp_size_t  sqr_toom8_threshold          = SQR_TOOM8_THRESHOLD_LIMIT;
  129. mp_size_t  sqr_fft_threshold            = MP_SIZE_T_MAX;
  130. mp_size_t  sqr_fft_modf_threshold       = MP_SIZE_T_MAX;
  131. mp_size_t  mullo_basecase_threshold     = MP_SIZE_T_MAX;
  132. mp_size_t  mullo_dc_threshold           = MP_SIZE_T_MAX;
  133. mp_size_t  mullo_mul_n_threshold        = MP_SIZE_T_MAX;
  134. mp_size_t  mulmod_bnm1_threshold        = MP_SIZE_T_MAX;
  135. mp_size_t  sqrmod_bnm1_threshold        = MP_SIZE_T_MAX;
  136. mp_size_t  div_sb_preinv_threshold      = MP_SIZE_T_MAX;
  137. mp_size_t  dc_div_qr_threshold          = MP_SIZE_T_MAX;
  138. mp_size_t  dc_divappr_q_threshold       = MP_SIZE_T_MAX;
  139. mp_size_t  mu_div_qr_threshold          = MP_SIZE_T_MAX;
  140. mp_size_t  mu_divappr_q_threshold       = MP_SIZE_T_MAX;
  141. mp_size_t  mupi_div_qr_threshold        = MP_SIZE_T_MAX;
  142. mp_size_t  mu_div_q_threshold           = MP_SIZE_T_MAX;
  143. mp_size_t  dc_bdiv_qr_threshold         = MP_SIZE_T_MAX;
  144. mp_size_t  dc_bdiv_q_threshold          = MP_SIZE_T_MAX;
  145. mp_size_t  mu_bdiv_qr_threshold         = MP_SIZE_T_MAX;
  146. mp_size_t  mu_bdiv_q_threshold          = MP_SIZE_T_MAX;
  147. mp_size_t  inv_mulmod_bnm1_threshold    = MP_SIZE_T_MAX;
  148. mp_size_t  inv_newton_threshold         = MP_SIZE_T_MAX;
  149. mp_size_t  inv_appr_threshold           = MP_SIZE_T_MAX;
  150. mp_size_t  binv_newton_threshold        = MP_SIZE_T_MAX;
  151. mp_size_t  redc_1_to_redc_2_threshold   = MP_SIZE_T_MAX;
  152. mp_size_t  redc_1_to_redc_n_threshold   = MP_SIZE_T_MAX;
  153. mp_size_t  redc_2_to_redc_n_threshold   = MP_SIZE_T_MAX;
  154. mp_size_t  powm_threshold               = MP_SIZE_T_MAX;
  155. mp_size_t  matrix22_strassen_threshold  = MP_SIZE_T_MAX;
  156. mp_size_t  hgcd_threshold               = MP_SIZE_T_MAX;
  157. mp_size_t  gcd_accel_threshold          = MP_SIZE_T_MAX;
  158. mp_size_t  gcd_dc_threshold             = MP_SIZE_T_MAX;
  159. mp_size_t  gcdext_dc_threshold          = MP_SIZE_T_MAX;
  160. mp_size_t  divrem_1_norm_threshold      = MP_SIZE_T_MAX;
  161. mp_size_t  divrem_1_unnorm_threshold    = MP_SIZE_T_MAX;
  162. mp_size_t  mod_1_norm_threshold         = MP_SIZE_T_MAX;
  163. mp_size_t  mod_1_unnorm_threshold       = MP_SIZE_T_MAX;
  164. mp_size_t  mod_1n_to_mod_1_1_threshold  = MP_SIZE_T_MAX;
  165. mp_size_t  mod_1u_to_mod_1_1_threshold  = MP_SIZE_T_MAX;
  166. mp_size_t  mod_1_1_to_mod_1_2_threshold = MP_SIZE_T_MAX;
  167. mp_size_t  mod_1_2_to_mod_1_4_threshold = MP_SIZE_T_MAX;
  168. mp_size_t  preinv_mod_1_to_mod_1_threshold = MP_SIZE_T_MAX;
  169. mp_size_t  divrem_2_threshold           = MP_SIZE_T_MAX;
  170. mp_size_t  get_str_dc_threshold         = MP_SIZE_T_MAX;
  171. mp_size_t  get_str_precompute_threshold = MP_SIZE_T_MAX;
  172. mp_size_t  set_str_dc_threshold         = MP_SIZE_T_MAX;
  173. mp_size_t  set_str_precompute_threshold = MP_SIZE_T_MAX;
  174. mp_size_t  fft_modf_sqr_threshold = MP_SIZE_T_MAX;
  175. mp_size_t  fft_modf_mul_threshold = MP_SIZE_T_MAX;
  176. struct param_t {
  177.   const char        *name;
  178.   speed_function_t  function;
  179.   speed_function_t  function2;
  180.   double            step_factor;    /* how much to step relatively */
  181.   int               step;           /* how much to step absolutely */
  182.   double            function_fudge; /* multiplier for "function" speeds */
  183.   int               stop_since_change;
  184.   double            stop_factor;
  185.   mp_size_t         min_size;
  186.   int               min_is_always;
  187.   mp_size_t         max_size;
  188.   mp_size_t         check_size;
  189.   mp_size_t         size_extra;
  190. #define DATA_HIGH_LT_R  1
  191. #define DATA_HIGH_GE_R  2
  192.   int               data_high;
  193.   int               noprint;
  194. };
  195. /* These are normally undefined when false, which suits "#if" fine.
  196.    But give them zero values so they can be used in plain C "if"s.  */
  197. #ifndef UDIV_PREINV_ALWAYS
  198. #define UDIV_PREINV_ALWAYS 0
  199. #endif
  200. #ifndef HAVE_NATIVE_mpn_divexact_1
  201. #define HAVE_NATIVE_mpn_divexact_1 0
  202. #endif
  203. #ifndef HAVE_NATIVE_mpn_divrem_1
  204. #define HAVE_NATIVE_mpn_divrem_1 0
  205. #endif
  206. #ifndef HAVE_NATIVE_mpn_divrem_2
  207. #define HAVE_NATIVE_mpn_divrem_2 0
  208. #endif
  209. #ifndef HAVE_NATIVE_mpn_mod_1
  210. #define HAVE_NATIVE_mpn_mod_1 0
  211. #endif
  212. #ifndef HAVE_NATIVE_mpn_modexact_1_odd
  213. #define HAVE_NATIVE_mpn_modexact_1_odd 0
  214. #endif
  215. #ifndef HAVE_NATIVE_mpn_preinv_divrem_1
  216. #define HAVE_NATIVE_mpn_preinv_divrem_1 0
  217. #endif
  218. #ifndef HAVE_NATIVE_mpn_preinv_mod_1
  219. #define HAVE_NATIVE_mpn_preinv_mod_1 0
  220. #endif
  221. #ifndef HAVE_NATIVE_mpn_sqr_basecase
  222. #define HAVE_NATIVE_mpn_sqr_basecase 0
  223. #endif
  224. #define MAX3(a,b,c)  MAX (MAX (a, b), c)
  225. mp_limb_t
  226. randlimb_norm (void)
  227. {
  228.   mp_limb_t  n;
  229.   mpn_random (&n, 1);
  230.   n |= GMP_NUMB_HIGHBIT;
  231.   return n;
  232. }
  233. #define GMP_NUMB_HALFMASK  ((CNST_LIMB(1) << (GMP_NUMB_BITS/2)) - 1)
  234. mp_limb_t
  235. randlimb_half (void)
  236. {
  237.   mp_limb_t  n;
  238.   mpn_random (&n, 1);
  239.   n &= GMP_NUMB_HALFMASK;
  240.   n += (n==0);
  241.   return n;
  242. }
  243. /* Add an entry to the end of the dat[] array, reallocing to make it bigger
  244.    if necessary.  */
  245. void
  246. add_dat (mp_size_t size, double d)
  247. {
  248. #define ALLOCDAT_STEP  500
  249.   ASSERT_ALWAYS (ndat <= allocdat);
  250.   if (ndat == allocdat)
  251.     {
  252.       dat = (struct dat_t *) __gmp_allocate_or_reallocate
  253.         (dat, allocdat * sizeof(dat[0]),
  254.          (allocdat+ALLOCDAT_STEP) * sizeof(dat[0]));
  255.       allocdat += ALLOCDAT_STEP;
  256.     }
  257.   dat[ndat].size = size;
  258.   dat[ndat].d = d;
  259.   ndat++;
  260. }
  261. /* Return the threshold size based on the data accumulated. */
  262. mp_size_t
  263. analyze_dat (int final)
  264. {
  265.   double  x, min_x;
  266.   int     j, min_j;
  267.   /* If the threshold is set at dat[0].size, any positive values are bad. */
  268.   x = 0.0;
  269.   for (j = 0; j < ndat; j++)
  270.     if (dat[j].d > 0.0)
  271.       x += dat[j].d;
  272.   if (option_trace >= 2 && final)
  273.     {
  274.       printf ("n");
  275.       printf ("x is the sum of the badness from setting thresh at given sizen");
  276.       printf ("  (minimum x is sought)n");
  277.       printf ("size=%ld  first x=%.4fn", (long) dat[j].size, x);
  278.     }
  279.   min_x = x;
  280.   min_j = 0;
  281.   /* When stepping to the next dat[j].size, positive values are no longer
  282.      bad (so subtracted), negative values become bad (so add the absolute
  283.      value, meaning subtract). */
  284.   for (j = 0; j < ndat; x -= dat[j].d, j++)
  285.     {
  286.       if (option_trace >= 2 && final)
  287.         printf ("size=%ld  x=%.4fn", (long) dat[j].size, x);
  288.       if (x < min_x)
  289.         {
  290.           min_x = x;
  291.           min_j = j;
  292.         }
  293.     }
  294.   return min_j;
  295. }
  296. /* Measuring for recompiled mpn/generic/divrem_1.c and mpn/generic/mod_1.c */
  297. mp_limb_t mpn_divrem_1_tune
  298.   __GMP_PROTO ((mp_ptr, mp_size_t, mp_srcptr, mp_size_t, mp_limb_t));
  299. mp_limb_t mpn_mod_1_tune
  300.    __GMP_PROTO ((mp_srcptr, mp_size_t, mp_limb_t));
  301. double
  302. speed_mpn_mod_1_tune (struct speed_params *s)
  303. {
  304.   SPEED_ROUTINE_MPN_MOD_1 (mpn_mod_1_tune);
  305. }
  306. double
  307. speed_mpn_divrem_1_tune (struct speed_params *s)
  308. {
  309.   SPEED_ROUTINE_MPN_DIVREM_1 (mpn_divrem_1_tune);
  310. }
  311. double
  312. tuneup_measure (speed_function_t fun,
  313.                 const struct param_t *param,
  314.                 struct speed_params *s)
  315. {
  316.   static struct param_t  dummy;
  317.   double   t;
  318.   TMP_DECL;
  319.   if (! param)
  320.     param = &dummy;
  321.   s->size += param->size_extra;
  322.   TMP_MARK;
  323.   SPEED_TMP_ALLOC_LIMBS (s->xp, s->size, 0);
  324.   SPEED_TMP_ALLOC_LIMBS (s->yp, s->size, 0);
  325.   mpn_random (s->xp, s->size);
  326.   mpn_random (s->yp, s->size);
  327.   switch (param->data_high) {
  328.   case DATA_HIGH_LT_R:
  329.     s->xp[s->size-1] %= s->r;
  330.     s->yp[s->size-1] %= s->r;
  331.     break;
  332.   case DATA_HIGH_GE_R:
  333.     s->xp[s->size-1] |= s->r;
  334.     s->yp[s->size-1] |= s->r;
  335.     break;
  336.   }
  337.   t = speed_measure (fun, s);
  338.   s->size -= param->size_extra;
  339.   TMP_FREE;
  340.   return t;
  341. }
  342. #define PRINT_WIDTH  31
  343. void
  344. print_define_start (const char *name)
  345. {
  346.   printf ("#define %-*s  ", PRINT_WIDTH, name);
  347.   if (option_trace)
  348.     printf ("...n");
  349. }
  350. void
  351. print_define_end_remark (const char *name, mp_size_t value, const char *remark)
  352. {
  353.   if (option_trace)
  354.     printf ("#define %-*s  ", PRINT_WIDTH, name);
  355.   if (value == MP_SIZE_T_MAX)
  356.     printf ("MP_SIZE_T_MAX");
  357.   else
  358.     printf ("%5ld", (long) value);
  359.   if (remark != NULL)
  360.     printf ("  /* %s */", remark);
  361.   printf ("n");
  362.   fflush (stdout);
  363. }
  364. void
  365. print_define_end (const char *name, mp_size_t value)
  366. {
  367.   const char  *remark;
  368.   if (value == MP_SIZE_T_MAX)
  369.     remark = "never";
  370.   else if (value == 0)
  371.     remark = "always";
  372.   else
  373.     remark = NULL;
  374.   print_define_end_remark (name, value, remark);
  375. }
  376. void
  377. print_define (const char *name, mp_size_t value)
  378. {
  379.   print_define_start (name);
  380.   print_define_end (name, value);
  381. }
  382. void
  383. print_define_remark (const char *name, mp_size_t value, const char *remark)
  384. {
  385.   print_define_start (name);
  386.   print_define_end_remark (name, value, remark);
  387. }
  388. void
  389. one (mp_size_t *threshold, struct param_t *param)
  390. {
  391.   int  since_positive, since_thresh_change;
  392.   int  thresh_idx, new_thresh_idx;
  393. #define DEFAULT(x,n)  do { if (! (x))  (x) = (n); } while (0)
  394.   DEFAULT (param->function_fudge, 1.0);
  395.   DEFAULT (param->function2, param->function);
  396.   DEFAULT (param->step_factor, 0.01);  /* small steps by default */
  397.   DEFAULT (param->step, 1);            /* small steps by default */
  398.   DEFAULT (param->stop_since_change, 80);
  399.   DEFAULT (param->stop_factor, 1.2);
  400.   DEFAULT (param->min_size, 10);
  401.   DEFAULT (param->max_size, DEFAULT_MAX_SIZE);
  402.   if (param->check_size != 0)
  403.     {
  404.       double   t1, t2;
  405.       s.size = param->check_size;
  406.       *threshold = s.size+1;
  407.       t1 = tuneup_measure (param->function, param, &s);
  408.       *threshold = s.size;
  409.       t2 = tuneup_measure (param->function2, param, &s);
  410.       if (t1 == -1.0 || t2 == -1.0)
  411.         {
  412.           printf ("Oops, can't run both functions at size %ldn",
  413.                   (long) s.size);
  414.           abort ();
  415.         }
  416.       t1 *= param->function_fudge;
  417.       /* ask that t2 is at least 4% below t1 */
  418.       if (t1 < t2*1.04)
  419.         {
  420.           if (option_trace)
  421.             printf ("function2 never enough faster: t1=%.9f t2=%.9fn", t1, t2);
  422.           *threshold = MP_SIZE_T_MAX;
  423.           if (! param->noprint)
  424.             print_define (param->name, *threshold);
  425.           return;
  426.         }
  427.       if (option_trace >= 2)
  428.         printf ("function2 enough faster at size=%ld: t1=%.9f t2=%.9fn",
  429.                 (long) s.size, t1, t2);
  430.     }
  431.   if (! param->noprint || option_trace)
  432.     print_define_start (param->name);
  433.   ndat = 0;
  434.   since_positive = 0;
  435.   since_thresh_change = 0;
  436.   thresh_idx = 0;
  437.   if (option_trace >= 2)
  438.     {
  439.       printf ("             algorithm-A  algorithm-B   ratio  possiblen");
  440.       printf ("              (seconds)    (seconds)    diff    threshn");
  441.     }
  442.   for (s.size = param->min_size;
  443.        s.size < param->max_size;
  444.        s.size += MAX ((mp_size_t) floor (s.size * param->step_factor), param->step))
  445.     {
  446.       double   ti, tiplus1, d;
  447.       /*
  448.         FIXME: check minimum size requirements are met, possibly by just
  449.         checking for the -1 returns from the speed functions.
  450.       */
  451.       /* using method A at this size */
  452.       *threshold = s.size+1;
  453.       ti = tuneup_measure (param->function, param, &s);
  454.       if (ti == -1.0)
  455.         abort ();
  456.       ti *= param->function_fudge;
  457.       /* using method B at this size */
  458.       *threshold = s.size;
  459.       tiplus1 = tuneup_measure (param->function2, param, &s);
  460.       if (tiplus1 == -1.0)
  461.         abort ();
  462.       /* Calculate the fraction by which the one or the other routine is
  463.          slower.  */
  464.       if (tiplus1 >= ti)
  465.         d = (tiplus1 - ti) / tiplus1;  /* negative */
  466.       else
  467.         d = (tiplus1 - ti) / ti;       /* positive */
  468.       add_dat (s.size, d);
  469.       new_thresh_idx = analyze_dat (0);
  470.       if (option_trace >= 2)
  471.         printf ("size=%ld  %.9f  %.9f  % .4f %c  %ldn",
  472.                 (long) s.size, ti, tiplus1, d,
  473.                 ti > tiplus1 ? '#' : ' ',
  474.                 (long) dat[new_thresh_idx].size);
  475.       /* Stop if the last time method i was faster was more than a
  476.          certain number of measurements ago.  */
  477. #define STOP_SINCE_POSITIVE  200
  478.       if (d >= 0)
  479.         since_positive = 0;
  480.       else
  481.         if (++since_positive > STOP_SINCE_POSITIVE)
  482.           {
  483.             if (option_trace >= 1)
  484.               printf ("stopped due to since_positive (%d)n",
  485.                       STOP_SINCE_POSITIVE);
  486.             break;
  487.           }
  488.       /* Stop if method A has become slower by a certain factor. */
  489.       if (ti >= tiplus1 * param->stop_factor)
  490.         {
  491.           if (option_trace >= 1)
  492.             printf ("stopped due to ti >= tiplus1 * factor (%.1f)n",
  493.                     param->stop_factor);
  494.           break;
  495.         }
  496.       /* Stop if the threshold implied hasn't changed in a certain
  497.          number of measurements.  (It's this condition that usually
  498.          stops the loop.) */
  499.       if (thresh_idx != new_thresh_idx)
  500.         since_thresh_change = 0, thresh_idx = new_thresh_idx;
  501.       else
  502.         if (++since_thresh_change > param->stop_since_change)
  503.           {
  504.             if (option_trace >= 1)
  505.               printf ("stopped due to since_thresh_change (%d)n",
  506.                       param->stop_since_change);
  507.             break;
  508.           }
  509.       /* Stop if the threshold implied is more than a certain number of
  510.          measurements ago.  */
  511. #define STOP_SINCE_AFTER   500
  512.       if (ndat - thresh_idx > STOP_SINCE_AFTER)
  513.         {
  514.           if (option_trace >= 1)
  515.             printf ("stopped due to ndat - thresh_idx > amount (%d)n",
  516.                     STOP_SINCE_AFTER);
  517.           break;
  518.         }
  519.       /* Stop when the size limit is reached before the end of the
  520.          crossover, but only show this as an error for >= the default max
  521.          size.  FIXME: Maybe should make it a param choice whether this is
  522.          an error.  */
  523.       if (s.size >= param->max_size && param->max_size >= DEFAULT_MAX_SIZE)
  524.         {
  525.           fprintf (stderr, "%sn", param->name);
  526.           fprintf (stderr, "sizes %ld to %ld total %d measurementsn",
  527.                    (long) dat[0].size, (long) dat[ndat-1].size, ndat);
  528.           fprintf (stderr, "    max size reached before end of crossovern");
  529.           break;
  530.         }
  531.     }
  532.   if (option_trace >= 1)
  533.     printf ("sizes %ld to %ld total %d measurementsn",
  534.             (long) dat[0].size, (long) dat[ndat-1].size, ndat);
  535.   *threshold = dat[analyze_dat (1)].size;
  536.   if (param->min_is_always)
  537.     {
  538.       if (*threshold == param->min_size)
  539.         *threshold = 0;
  540.     }
  541.   if (! param->noprint || option_trace)
  542.     print_define_end (param->name, *threshold);
  543. }
  544. /* Special probing for the fft thresholds.  The size restrictions on the
  545.    FFTs mean the graph of time vs size has a step effect.  See this for
  546.    example using
  547.        ./speed -s 4096-16384 -t 128 -P foo mpn_mul_fft.8 mpn_mul_fft.9
  548.        gnuplot foo.gnuplot
  549.    The current approach is to compare routines at the midpoint of relevant
  550.    steps.  Arguably a more sophisticated system of threshold data is wanted
  551.    if this step effect remains. */
  552. struct fft_param_t {
  553.   const char        *table_name;
  554.   const char        *threshold_name;
  555.   const char        *modf_threshold_name;
  556.   mp_size_t         *p_threshold;
  557.   mp_size_t         *p_modf_threshold;
  558.   mp_size_t         first_size;
  559.   mp_size_t         max_size;
  560.   speed_function_t  function;
  561.   speed_function_t  mul_modf_function;
  562.   speed_function_t  mul_function;
  563.   mp_size_t         sqr;
  564. };
  565. /* mpn_mul_fft requires pl a multiple of 2^k limbs, but with
  566.    N=pl*BIT_PER_MP_LIMB it internally also pads out so N/2^k is a multiple
  567.    of 2^(k-1) bits. */
  568. mp_size_t
  569. fft_step_size (int k)
  570. {
  571.   mp_size_t  step;
  572.   step = MAX ((mp_size_t) 1 << (k-1), GMP_LIMB_BITS) / GMP_LIMB_BITS;
  573.   step *= (mp_size_t) 1 << k;
  574.   if (step <= 0)
  575.     {
  576.       printf ("Can't handle k=%dn", k);
  577.       abort ();
  578.     }
  579.   return step;
  580. }
  581. mp_size_t
  582. fft_next_size (mp_size_t pl, int k)
  583. {
  584.   mp_size_t  m = fft_step_size (k);
  585. /*    printf ("[k=%d %ld] %ld ->", k, m, pl); */
  586.   if (pl == 0 || (pl & (m-1)) != 0)
  587.     pl = (pl | (m-1)) + 1;
  588. /*    printf (" %ldn", pl); */
  589.   return pl;
  590. }
  591. #define NMAX_DEFAULT 1000000
  592. #define MAX_REPS 25
  593. #define MIN_REPS 5
  594. static inline size_t
  595. mpn_mul_fft_lcm (size_t a, unsigned int k)
  596. {
  597.   unsigned int l = k;
  598.   while (a % 2 == 0 && k > 0)
  599.     {
  600.       a >>= 1;
  601.       k--;
  602.     }
  603.   return a << l;
  604. }
  605. mp_size_t
  606. fftfill (mp_size_t pl, int k, int sqr)
  607. {
  608.   mp_size_t maxLK;
  609.   mp_bitcnt_t N, Nprime, nprime, M;
  610.   N = pl * GMP_NUMB_BITS;
  611.   M = N >> k;
  612.   maxLK = mpn_mul_fft_lcm ((unsigned long) GMP_NUMB_BITS, k);
  613.   Nprime = (1 + (2 * M + k + 2) / maxLK) * maxLK;
  614.   nprime = Nprime / GMP_NUMB_BITS;
  615.   if (nprime >= (sqr ? SQR_FFT_MODF_THRESHOLD : MUL_FFT_MODF_THRESHOLD))
  616.     {
  617.       size_t K2;
  618.       for (;;)
  619. {
  620.   K2 = 1L << mpn_fft_best_k (nprime, sqr);
  621.   if ((nprime & (K2 - 1)) == 0)
  622.     break;
  623.   nprime = (nprime + K2 - 1) & -K2;
  624.   Nprime = nprime * GMP_LIMB_BITS;
  625. }
  626.     }
  627.   ASSERT_ALWAYS (nprime < pl);
  628.   return Nprime;
  629. }
  630. static int
  631. compare_double (const void *ap, const void *bp)
  632. {
  633.   double a = * (const double *) ap;
  634.   double b = * (const double *) bp;
  635.   if (a < b)
  636.     return -1;
  637.   else if (a > b)
  638.     return 1;
  639.   else
  640.     return 0;
  641. }
  642. double
  643. median (double *times, int n)
  644. {
  645.   qsort (times, n, sizeof (double), compare_double);
  646.   return times[n/2];
  647. }
  648. #define FFT_CACHE_SIZE 25
  649. typedef struct fft_cache
  650. {
  651.   mp_size_t n;
  652.   double time;
  653. } fft_cache_t;
  654. fft_cache_t fft_cache[FFT_CACHE_SIZE];
  655. double
  656. cached_measure (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n, int k,
  657. int n_measurements)
  658. {
  659.   int i;
  660.   double t, ttab[MAX_REPS];
  661.   if (fft_cache[k].n == n)
  662.     return fft_cache[k].time;
  663.   for (i = 0; i < n_measurements; i++)
  664.     {
  665.       speed_starttime ();
  666.       mpn_mul_fft (rp, n, ap, n, bp, n, k);
  667.       ttab[i] = speed_endtime ();
  668.     }
  669.   t = median (ttab, n_measurements);
  670.   fft_cache[k].n = n;
  671.   fft_cache[k].time = t;
  672.   return t;
  673. }
  674. #define INSERT_FFTTAB(idx, nval, kval)
  675.   do {
  676.     fft_tab[idx].n = nval;
  677.     fft_tab[idx].k = kval;
  678.     fft_tab[idx+1].n = -1; /* sentinel */
  679.     fft_tab[idx+1].k = -1;
  680.   } while (0)
  681. int
  682. fftmes (mp_size_t nmin, mp_size_t nmax, int initial_k, struct fft_param_t *p, int idx, int print)
  683. {
  684.   mp_size_t n, n1, prev_n1;
  685.   int k, best_k, last_best_k, kmax;
  686.   int eff, prev_eff;
  687.   double t0, t1;
  688.   int n_measurements;
  689.   mp_limb_t *ap, *bp, *rp;
  690.   mp_size_t alloc;
  691.   char *linepref;
  692.   struct fft_table_nk *fft_tab;
  693.   fft_tab = mpn_fft_table3[p->sqr];
  694.   for (k = 0; k < FFT_CACHE_SIZE; k++)
  695.     fft_cache[k].n = 0;
  696.   if (nmin < (p->sqr ? SQR_FFT_MODF_THRESHOLD : MUL_FFT_MODF_THRESHOLD))
  697.     {
  698.       nmin = (p->sqr ? SQR_FFT_MODF_THRESHOLD : MUL_FFT_MODF_THRESHOLD);
  699.     }
  700.   if (print)
  701.     printf ("#define %s%*s", p->table_name, 38, "");
  702.   if (idx == 0)
  703.     {
  704.       INSERT_FFTTAB (0, nmin, initial_k);
  705.       if (print)
  706. {
  707.   printf ("\n  { ");
  708.   printf ("{%7u,%2u}", fft_tab[0].n, fft_tab[0].k);
  709.   linepref = "    ";
  710. }
  711.       idx = 1;
  712.     }
  713.   ap = malloc (sizeof (mp_limb_t));
  714.   if (p->sqr)
  715.     bp = ap;
  716.   else
  717.     bp = malloc (sizeof (mp_limb_t));
  718.   rp = malloc (sizeof (mp_limb_t));
  719.   alloc = 1;
  720.   /* Round n to comply to initial k value */
  721.   n = (nmin + ((1ul << initial_k) - 1)) & (MP_SIZE_T_MAX << initial_k);
  722.   n_measurements = (18 - initial_k) | 1;
  723.   n_measurements = MAX (n_measurements, MIN_REPS);
  724.   n_measurements = MIN (n_measurements, MAX_REPS);
  725.   last_best_k = initial_k;
  726.   best_k = initial_k;
  727.   while (n < nmax)
  728.     {
  729.       int start_k, end_k;
  730.       /* Assume the current best k is best until we hit its next FFT step.  */
  731.       t0 = 99999;
  732.       prev_n1 = n + 1;
  733.       start_k = MAX (4, best_k - 4);
  734.       end_k = MIN (24, best_k + 4);
  735.       for (k = start_k; k <= end_k; k++)
  736. {
  737.           n1 = mpn_fft_next_size (prev_n1, k);
  738.   eff = 200 * (n1 * GMP_NUMB_BITS >> k) / fftfill (n1, k, p->sqr);
  739.   if (eff < 70) /* avoid measuring too slow fft:s */
  740.     continue;
  741.   if (n1 > alloc)
  742.     {
  743.       alloc = n1;
  744.       if (p->sqr)
  745. {
  746.   ap = realloc (ap, sizeof (mp_limb_t));
  747.   rp = realloc (rp, sizeof (mp_limb_t));
  748.   ap = bp = realloc (ap, alloc * sizeof (mp_limb_t));
  749.   mpn_random (ap, alloc);
  750.   rp = realloc (rp, alloc * sizeof (mp_limb_t));
  751. }
  752.       else
  753. {
  754.   ap = realloc (ap, sizeof (mp_limb_t));
  755.   bp = realloc (bp, sizeof (mp_limb_t));
  756.   rp = realloc (rp, sizeof (mp_limb_t));
  757.   ap = realloc (ap, alloc * sizeof (mp_limb_t));
  758.   mpn_random (ap, alloc);
  759.   bp = realloc (bp, alloc * sizeof (mp_limb_t));
  760.   mpn_random (bp, alloc);
  761.   rp = realloc (rp, alloc * sizeof (mp_limb_t));
  762. }
  763.     }
  764.   t1 = cached_measure (rp, ap, bp, n1, k, n_measurements);
  765.   if (t1 * n_measurements > 0.3)
  766.     n_measurements -= 2;
  767.   n_measurements = MAX (n_measurements, MIN_REPS);
  768.   if (t1 < t0)
  769.     {
  770.       best_k = k;
  771.       t0 = t1;
  772.     }
  773. }
  774.       n1 = mpn_fft_next_size (prev_n1, best_k);
  775.       if (last_best_k != best_k)
  776. {
  777.   ASSERT_ALWAYS ((prev_n1 & ((1ul << last_best_k) - 1)) == 1);
  778.   if (idx >= FFT_TABLE3_SIZE)
  779.     {
  780.       printf ("FFT table exhausted, increase FFT_TABLE3_SIZE in gmp-impl.hn");
  781.       abort ();
  782.     }
  783.   INSERT_FFTTAB (idx, prev_n1 >> last_best_k, best_k);
  784.   if (print)
  785.     {
  786.       printf (", ");
  787.       if (idx % 4 == 0)
  788. printf ("\n    ");
  789.       printf ("{%7u,%2u}", fft_tab[idx].n, fft_tab[idx].k);
  790.     }
  791.   if (option_trace >= 2)
  792.     {
  793.       printf ("{%lu,%u}n", prev_n1, best_k);
  794.       fflush (stdout);
  795.     }
  796.   last_best_k = best_k;
  797.   idx++;
  798. }
  799.       for (;;)
  800. {
  801.   prev_n1 = n1;
  802.   prev_eff = fftfill (prev_n1, best_k, p->sqr);
  803.   n1 = mpn_fft_next_size (prev_n1 + 1, best_k);
  804.   eff = fftfill (n1, best_k, p->sqr);
  805.   if (eff != prev_eff)
  806.     break;
  807. }
  808.       n = prev_n1;
  809.     }
  810.   kmax = sizeof (mp_size_t) * 4; /* GMP_MP_SIZE_T_BITS / 2 */
  811.   kmax = MIN (kmax, 25-1);
  812.   for (k = last_best_k + 1; k <= kmax; k++)
  813.     {
  814.       if (idx >= FFT_TABLE3_SIZE)
  815. {
  816.   printf ("FFT table exhausted, increase FFT_TABLE3_SIZE in gmp-impl.hn");
  817.   abort ();
  818. }
  819.       INSERT_FFTTAB (idx, ((1ul << (2*k-2)) + 1) >> (k-1), k);
  820.       if (print)
  821. {
  822.   printf (", ");
  823.   if (idx % 4 == 0)
  824.     printf ("\n    ");
  825.   printf ("{%7u,%2u}", fft_tab[idx].n, fft_tab[idx].k);
  826. }
  827.       idx++;
  828.     }
  829.   if (print)
  830.     printf (" }n");
  831.   free (ap);
  832.   if (! p->sqr)
  833.     free (bp);
  834.   free (rp);
  835.   return idx;
  836. }
  837. void
  838. fft (struct fft_param_t *p)
  839. {
  840.   mp_size_t  size;
  841.   int        k, idx, initial_k;
  842.   /*** Generate MUL_FFT_MODF_THRESHOLD / SQR_FFT_MODF_THRESHOLD ***/
  843. #if 1
  844.   {
  845.     /* Use plain one() mechanism, for some reasonable initial values of k.  The
  846.        advantage is that we don't depend on mpn_fft_table3, which can therefore
  847.        leave it completely uninitialized.  */
  848.     static struct param_t param;
  849.     mp_size_t thres, best_thres;
  850.     int best_k;
  851.     char buf[20];
  852.     best_thres = MP_SIZE_T_MAX;
  853.     best_k = -1;
  854.     for (k = 5; k <= 7; k++)
  855.       {
  856. param.name = p->modf_threshold_name;
  857. param.min_size = 100;
  858. param.max_size = 2000;
  859. param.function  = p->mul_function;
  860. param.step_factor = 0.0;
  861. param.step = 4;
  862. param.function2 = p->mul_modf_function;
  863. param.noprint = 1;
  864. s.r = k;
  865. one (&thres, &param);
  866. if (thres < best_thres)
  867.   {
  868.     best_thres = thres;
  869.     best_k = k;
  870.   }
  871.       }
  872.     *(p->p_modf_threshold) = best_thres;
  873.     sprintf (buf, "k = %d", best_k);
  874.     print_define_remark (p->modf_threshold_name, best_thres, buf);
  875.     initial_k = best_k;
  876.   }
  877. #else
  878.   size = p->first_size;
  879.   for (;;)
  880.     {
  881.       double  tk, tm;
  882.       size = mpn_fft_next_size (size+1, mpn_fft_best_k (size+1, p->sqr));
  883.       k = mpn_fft_best_k (size, p->sqr);
  884.       if (size >= p->max_size)
  885.         break;
  886.       s.size = size + fft_step_size (k) / 2;
  887.       s.r = k;
  888.       tk = tuneup_measure (p->mul_modf_function, NULL, &s);
  889.       if (tk == -1.0)
  890.         abort ();
  891.       tm = tuneup_measure (p->mul_function, NULL, &s);
  892.       if (tm == -1.0)
  893.         abort ();
  894.       if (option_trace >= 2)
  895.         printf ("at %ld   size=%ld  k=%d  %.9f   size=%ld modf %.9fn",
  896.                 (long) size,
  897.                 (long) size + fft_step_size (k) / 2, k, tk,
  898.                 (long) s.size, tm);
  899.       if (tk < tm)
  900.         {
  901.   *p->p_modf_threshold = s.size;
  902.   print_define (p->modf_threshold_name, *p->p_modf_threshold);
  903.   break;
  904.         }
  905.     }
  906.   initial_k = ?;
  907. #endif
  908.   /*** Generate MUL_FFT_TABLE3 / SQR_FFT_TABLE3 ***/
  909.   idx = fftmes (*p->p_modf_threshold, p->max_size, initial_k, p, 0, 1);
  910.   printf ("#define %s_SIZE %dn", p->table_name, idx);
  911.   /*** Generate MUL_FFT_THRESHOLD / SQR_FFT_THRESHOLD ***/
  912.   size = 2 * *p->p_modf_threshold; /* OK? */
  913.   for (;;)
  914.     {
  915.       double  tk, tm;
  916.       mp_size_t mulmod_size, mul_size;;
  917.       if (size >= p->max_size)
  918.         break;
  919.       mulmod_size = mpn_mulmod_bnm1_next_size (2 * (size + 1)) / 2;
  920.       mul_size = (size + mulmod_size) / 2; /* middle of step */
  921.       s.size = mulmod_size;
  922.       tk = tuneup_measure (p->function, NULL, &s);
  923.       if (tk == -1.0)
  924.         abort ();
  925.       s.size = mul_size;
  926.       tm = tuneup_measure (p->mul_function, NULL, &s);
  927.       if (tm == -1.0)
  928.         abort ();
  929.       if (option_trace >= 2)
  930.         printf ("at %ld   size=%ld  %.9f   size=%ld mul %.9fn",
  931.                 (long) size,
  932.                 (long) mulmod_size, tk,
  933.                 (long) mul_size, tm);
  934.       size = mulmod_size;
  935.       if (tk < tm)
  936.         {
  937.   *p->p_threshold = s.size;
  938.   print_define (p->threshold_name, *p->p_threshold);
  939.   break;
  940.         }
  941.     }
  942. }
  943. /* Start karatsuba from 4, since the Cray t90 ieee code is much faster at 2,
  944.    giving wrong results.  */
  945. void
  946. tune_mul_n (void)
  947. {
  948.   static struct param_t  param;
  949.   param.function = speed_mpn_mul_n;
  950.   param.name = "MUL_TOOM22_THRESHOLD";
  951.   param.min_size = MAX (4, MPN_TOOM22_MUL_MINSIZE);
  952.   param.max_size = MUL_TOOM22_THRESHOLD_LIMIT-1;
  953.   one (&mul_toom22_threshold, &param);
  954.   param.name = "MUL_TOOM33_THRESHOLD";
  955.   param.min_size = MAX (mul_toom22_threshold, MPN_TOOM33_MUL_MINSIZE);
  956.   param.max_size = MUL_TOOM33_THRESHOLD_LIMIT-1;
  957.   one (&mul_toom33_threshold, &param);
  958.   param.name = "MUL_TOOM44_THRESHOLD";
  959.   param.min_size = MAX (mul_toom33_threshold, MPN_TOOM44_MUL_MINSIZE);
  960.   param.max_size = MUL_TOOM44_THRESHOLD_LIMIT-1;
  961.   one (&mul_toom44_threshold, &param);
  962.   param.name = "MUL_TOOM6H_THRESHOLD";
  963.   param.min_size = MAX (mul_toom44_threshold, MPN_TOOM6H_MUL_MINSIZE);
  964.   param.max_size = MUL_TOOM6H_THRESHOLD_LIMIT-1;
  965.   one (&mul_toom6h_threshold, &param);
  966.   param.name = "MUL_TOOM8H_THRESHOLD";
  967.   param.min_size = MAX (mul_toom6h_threshold, MPN_TOOM8H_MUL_MINSIZE);
  968.   param.max_size = MUL_TOOM8H_THRESHOLD_LIMIT-1;
  969.   one (&mul_toom8h_threshold, &param);
  970.   /* disabled until tuned */
  971.   MUL_FFT_THRESHOLD = MP_SIZE_T_MAX;
  972. }
  973. void
  974. tune_mul (void)
  975. {
  976.   static struct param_t  param;
  977.   mp_size_t thres;
  978.   param.noprint = 1;
  979.   param.function = speed_mpn_toom32_for_toom43_mul;
  980.   param.function2 = speed_mpn_toom43_for_toom32_mul;
  981.   param.name = "MUL_TOOM32_TO_TOOM43_THRESHOLD";
  982.   param.min_size = MPN_TOOM43_MUL_MINSIZE;
  983.   one (&thres, &param);
  984.   mul_toom32_to_toom43_threshold = 17*thres/24;
  985.   print_define ("MUL_TOOM32_TO_TOOM43_THRESHOLD", mul_toom32_to_toom43_threshold);
  986.   param.function = speed_mpn_toom32_for_toom53_mul;
  987.   param.function2 = speed_mpn_toom53_for_toom32_mul;
  988.   param.name = "MUL_TOOM32_TO_TOOM53_THRESHOLD";
  989.   param.min_size = MPN_TOOM53_MUL_MINSIZE;
  990.   one (&thres, &param);
  991.   mul_toom32_to_toom53_threshold = 19*thres/30;
  992.   print_define ("MUL_TOOM32_TO_TOOM53_THRESHOLD", mul_toom32_to_toom53_threshold);
  993.   param.function = speed_mpn_toom42_for_toom53_mul;
  994.   param.function2 = speed_mpn_toom53_for_toom42_mul;
  995.   param.name = "MUL_TOOM42_TO_TOOM53_THRESHOLD";
  996.   param.min_size = MPN_TOOM53_MUL_MINSIZE;
  997.   one (&thres, &param);
  998.   mul_toom42_to_toom53_threshold = 11*thres/20;
  999.   print_define ("MUL_TOOM42_TO_TOOM53_THRESHOLD", mul_toom42_to_toom53_threshold);
  1000.   param.function = speed_mpn_toom42_mul;
  1001.   param.function2 = speed_mpn_toom63_mul;
  1002.   param.name = "MUL_TOOM42_TO_TOOM63_THRESHOLD";
  1003.   param.min_size = MPN_TOOM63_MUL_MINSIZE;
  1004.   one (&thres, &param);
  1005.   mul_toom42_to_toom63_threshold = thres/2;
  1006.   print_define ("MUL_TOOM42_TO_TOOM63_THRESHOLD", mul_toom42_to_toom63_threshold);
  1007. }
  1008. void
  1009. tune_mullo (void)
  1010. {
  1011.   static struct param_t  param;
  1012.   param.function = speed_mpn_mullo_n;
  1013.   param.name = "MULLO_BASECASE_THRESHOLD";
  1014.   param.min_size = 1;
  1015.   param.min_is_always = 1;
  1016.   param.max_size = MULLO_BASECASE_THRESHOLD_LIMIT-1;
  1017.   param.stop_factor = 1.5;
  1018.   param.noprint = 1;
  1019.   one (&mullo_basecase_threshold, &param);
  1020.   param.name = "MULLO_DC_THRESHOLD";
  1021.   param.min_size = 8;
  1022.   param.min_is_always = 0;
  1023.   param.max_size = 1000;
  1024.   one (&mullo_dc_threshold, &param);
  1025.   if (mullo_basecase_threshold >= mullo_dc_threshold)
  1026.     {
  1027.       print_define ("MULLO_BASECASE_THRESHOLD", mullo_dc_threshold);
  1028.       print_define_remark ("MULLO_DC_THRESHOLD", 0, "never mpn_mullo_basecase");
  1029.     }
  1030.   else
  1031.     {
  1032.       print_define ("MULLO_BASECASE_THRESHOLD", mullo_basecase_threshold);
  1033.       print_define ("MULLO_DC_THRESHOLD", mullo_dc_threshold);
  1034.     }
  1035. #if WANT_FFT
  1036.   param.name = "MULLO_MUL_N_THRESHOLD";
  1037.   param.min_size = mullo_dc_threshold;
  1038.   param.max_size = 2 * mul_fft_threshold;
  1039.   param.noprint = 0;
  1040.   param.step_factor = 0.03;
  1041.   one (&mullo_mul_n_threshold, &param);
  1042. #else
  1043.   print_define_remark ("MULLO_MUL_N_THRESHOLD", MP_SIZE_T_MAX,
  1044.                            "without FFT use mullo forever");
  1045. #endif
  1046. }
  1047. void
  1048. tune_mulmod_bnm1 (void)
  1049. {
  1050.   static struct param_t  param;
  1051.   param.name = "MULMOD_BNM1_THRESHOLD";
  1052.   param.function = speed_mpn_mulmod_bnm1;
  1053.   param.min_size = 4;
  1054.   param.max_size = 100;
  1055.   one (&mulmod_bnm1_threshold, &param);
  1056. }
  1057. void
  1058. tune_sqrmod_bnm1 (void)
  1059. {
  1060.   static struct param_t  param;
  1061.   param.name = "SQRMOD_BNM1_THRESHOLD";
  1062.   param.function = speed_mpn_sqrmod_bnm1;
  1063.   param.min_size = 4;
  1064.   param.max_size = 100;
  1065.   one (&sqrmod_bnm1_threshold, &param);
  1066. }
  1067. /* Start the basecase from 3, since 1 is a special case, and if mul_basecase
  1068.    is faster only at size==2 then we don't want to bother with extra code
  1069.    just for that.  Start karatsuba from 4 same as MUL above.  */
  1070. void
  1071. tune_sqr (void)
  1072. {
  1073.   /* disabled until tuned */
  1074.   SQR_FFT_THRESHOLD = MP_SIZE_T_MAX;
  1075.   if (HAVE_NATIVE_mpn_sqr_basecase)
  1076.     {
  1077.       print_define_remark ("SQR_BASECASE_THRESHOLD", 0, "always (native)");
  1078.       sqr_basecase_threshold = 0;
  1079.     }
  1080.   else
  1081.     {
  1082.       static struct param_t  param;
  1083.       param.name = "SQR_BASECASE_THRESHOLD";
  1084.       param.function = speed_mpn_sqr;
  1085.       param.min_size = 3;
  1086.       param.min_is_always = 1;
  1087.       param.max_size = TUNE_SQR_TOOM2_MAX;
  1088.       param.noprint = 1;
  1089.       one (&sqr_basecase_threshold, &param);
  1090.     }
  1091.   {
  1092.     static struct param_t  param;
  1093.     param.name = "SQR_TOOM2_THRESHOLD";
  1094.     param.function = speed_mpn_sqr;
  1095.     param.min_size = MAX (4, MPN_TOOM2_SQR_MINSIZE);
  1096.     param.max_size = TUNE_SQR_TOOM2_MAX;
  1097.     param.noprint = 1;
  1098.     one (&sqr_toom2_threshold, &param);
  1099.     if (! HAVE_NATIVE_mpn_sqr_basecase
  1100.         && sqr_toom2_threshold < sqr_basecase_threshold)
  1101.       {
  1102.         /* Karatsuba becomes faster than mul_basecase before
  1103.            sqr_basecase does.  Arrange for the expression
  1104.            "BELOW_THRESHOLD (un, SQR_TOOM2_THRESHOLD))" which
  1105.            selects mpn_sqr_basecase in mpn_sqr to be false, by setting
  1106.            SQR_TOOM2_THRESHOLD to zero, making
  1107.            SQR_BASECASE_THRESHOLD the toom2 threshold.  */
  1108.         sqr_basecase_threshold = SQR_TOOM2_THRESHOLD;
  1109.         SQR_TOOM2_THRESHOLD = 0;
  1110.         print_define_remark ("SQR_BASECASE_THRESHOLD", sqr_basecase_threshold,
  1111.                              "toom2");
  1112.         print_define_remark ("SQR_TOOM2_THRESHOLD",SQR_TOOM2_THRESHOLD,
  1113.                              "never sqr_basecase");
  1114.       }
  1115.     else
  1116.       {
  1117.         if (! HAVE_NATIVE_mpn_sqr_basecase)
  1118.           print_define ("SQR_BASECASE_THRESHOLD", sqr_basecase_threshold);
  1119.         print_define ("SQR_TOOM2_THRESHOLD", SQR_TOOM2_THRESHOLD);
  1120.       }
  1121.   }
  1122.   {
  1123.     static struct param_t  param;
  1124.     mp_size_t toom3_start = MAX (sqr_toom2_threshold, sqr_basecase_threshold);
  1125.     param.function = speed_mpn_sqr;
  1126.     param.name = "SQR_TOOM3_THRESHOLD";
  1127.     param.min_size = MAX (toom3_start, MPN_TOOM3_SQR_MINSIZE);
  1128.     param.max_size = SQR_TOOM3_THRESHOLD_LIMIT-1;
  1129.     one (&sqr_toom3_threshold, &param);
  1130.     param.name = "SQR_TOOM4_THRESHOLD";
  1131.     param.min_size = MAX (sqr_toom3_threshold, MPN_TOOM4_SQR_MINSIZE);
  1132.     param.max_size = SQR_TOOM4_THRESHOLD_LIMIT-1;
  1133.     one (&sqr_toom4_threshold, &param);
  1134.     param.name = "SQR_TOOM6_THRESHOLD";
  1135.     param.min_size = MAX (sqr_toom4_threshold, MPN_TOOM6_SQR_MINSIZE);
  1136.     param.max_size = SQR_TOOM6_THRESHOLD_LIMIT-1;
  1137.     one (&sqr_toom6_threshold, &param);
  1138.     param.name = "SQR_TOOM8_THRESHOLD";
  1139.     param.min_size = MAX (sqr_toom6_threshold, MPN_TOOM8_SQR_MINSIZE);
  1140.     param.max_size = SQR_TOOM8_THRESHOLD_LIMIT-1;
  1141.     one (&sqr_toom8_threshold, &param);
  1142.   }
  1143. }
  1144. void
  1145. tune_dc_div (void)
  1146. {
  1147.   s.r = 0; /* clear to make speed function do 2n/n */
  1148.   {
  1149.     static struct param_t  param;
  1150.     param.name = "DC_DIV_QR_THRESHOLD";
  1151.     param.function = speed_mpn_sbpi1_div_qr;
  1152.     param.function2 = speed_mpn_dcpi1_div_qr;
  1153.     param.min_size = 6;
  1154.     one (&dc_div_qr_threshold, &param);
  1155.   }
  1156.   {
  1157.     static struct param_t  param;
  1158.     param.name = "DC_DIVAPPR_Q_THRESHOLD";
  1159.     param.function = speed_mpn_sbpi1_divappr_q;
  1160.     param.function2 = speed_mpn_dcpi1_divappr_q;
  1161.     param.min_size = 6;
  1162.     one (&dc_divappr_q_threshold, &param);
  1163.   }
  1164. }
  1165. static double
  1166. speed_mpn_sbordcpi1_div_qr (struct speed_params *s)
  1167. {
  1168.   if (s->size < DC_DIV_QR_THRESHOLD)
  1169.     return speed_mpn_sbpi1_div_qr (s);
  1170.   else
  1171.     return speed_mpn_dcpi1_div_qr (s);
  1172. }
  1173. void
  1174. tune_mu_div (void)
  1175. {
  1176.   s.r = 0; /* clear to make speed function do 2n/n */
  1177.   {
  1178.     static struct param_t  param;
  1179.     param.name = "MU_DIV_QR_THRESHOLD";
  1180.     param.function = speed_mpn_dcpi1_div_qr;
  1181.     param.function2 = speed_mpn_mu_div_qr;
  1182.     param.min_size = 6;
  1183.     param.max_size = 5000;
  1184.     param.step_factor = 0.02;
  1185.     one (&mu_div_qr_threshold, &param);
  1186.   }
  1187.   {
  1188.     static struct param_t  param;
  1189.     param.name = "MU_DIVAPPR_Q_THRESHOLD";
  1190.     param.function = speed_mpn_dcpi1_divappr_q;
  1191.     param.function2 = speed_mpn_mu_divappr_q;
  1192.     param.min_size = 6;
  1193.     param.max_size = 5000;
  1194.     param.step_factor = 0.02;
  1195.     one (&mu_divappr_q_threshold, &param);
  1196.   }
  1197.   {
  1198.     static struct param_t  param;
  1199.     param.name = "MUPI_DIV_QR_THRESHOLD";
  1200.     param.function = speed_mpn_sbordcpi1_div_qr;
  1201.     param.function2 = speed_mpn_mupi_div_qr;
  1202.     param.min_size = 6;
  1203.     param.min_is_always = 1;
  1204.     param.max_size = 1000;
  1205.     param.step_factor = 0.02;
  1206.     one (&mupi_div_qr_threshold, &param);
  1207.   }
  1208. }
  1209. void
  1210. tune_dc_bdiv (void)
  1211. {
  1212.   s.r = 0; /* clear to make speed function do 2n/n*/
  1213.   {
  1214.     static struct param_t  param;
  1215.     param.name = "DC_BDIV_QR_THRESHOLD";
  1216.     param.function = speed_mpn_sbpi1_bdiv_qr;
  1217.     param.function2 = speed_mpn_dcpi1_bdiv_qr;
  1218.     param.min_size = 4;
  1219.     one (&dc_bdiv_qr_threshold, &param);
  1220.   }
  1221.   {
  1222.     static struct param_t  param;
  1223.     param.name = "DC_BDIV_Q_THRESHOLD";
  1224.     param.function = speed_mpn_sbpi1_bdiv_q;
  1225.     param.function2 = speed_mpn_dcpi1_bdiv_q;
  1226.     param.min_size = 4;
  1227.     one (&dc_bdiv_q_threshold, &param);
  1228.   }
  1229. }
  1230. void
  1231. tune_mu_bdiv (void)
  1232. {
  1233.   s.r = 0; /* clear to make speed function do 2n/n*/
  1234.   {
  1235.     static struct param_t  param;
  1236.     param.name = "MU_BDIV_QR_THRESHOLD";
  1237.     param.function = speed_mpn_dcpi1_bdiv_qr;
  1238.     param.function2 = speed_mpn_mu_bdiv_qr;
  1239.     param.min_size = 4;
  1240.     param.max_size = 5000;
  1241.     param.step_factor = 0.02;
  1242.     one (&mu_bdiv_qr_threshold, &param);
  1243.   }
  1244.   {
  1245.     static struct param_t  param;
  1246.     param.name = "MU_BDIV_Q_THRESHOLD";
  1247.     param.function = speed_mpn_dcpi1_bdiv_q;
  1248.     param.function2 = speed_mpn_mu_bdiv_q;
  1249.     param.min_size = 4;
  1250.     param.max_size = 5000;
  1251.     param.step_factor = 0.02;
  1252.     one (&mu_bdiv_q_threshold, &param);
  1253.   }
  1254. }
  1255. void
  1256. tune_invertappr (void)
  1257. {
  1258.   static struct param_t  param;
  1259.   param.function = speed_mpn_ni_invertappr;
  1260.   param.name = "INV_MULMOD_BNM1_THRESHOLD";
  1261.   param.min_size = 4;
  1262.   one (&inv_mulmod_bnm1_threshold, &param);
  1263.   param.function = speed_mpn_invertappr;
  1264.   param.name = "INV_NEWTON_THRESHOLD";
  1265.   param.min_size = 3;
  1266.   one (&inv_newton_threshold, &param);
  1267. }
  1268. void
  1269. tune_invert (void)
  1270. {
  1271.   static struct param_t  param;
  1272.   param.function = speed_mpn_invert;
  1273.   param.name = "INV_APPR_THRESHOLD";
  1274.   param.min_size = 3;
  1275.   one (&inv_appr_threshold, &param);
  1276. }
  1277. void
  1278. tune_binvert (void)
  1279. {
  1280.   static struct param_t  param;
  1281.   param.function = speed_mpn_binvert;
  1282.   param.name = "BINV_NEWTON_THRESHOLD";
  1283.   param.min_size = 8; /* pointless with smaller operands */
  1284.   one (&binv_newton_threshold, &param);
  1285. }
  1286. void
  1287. tune_redc (void)
  1288. {
  1289. #define TUNE_REDC_2_MAX 100
  1290. #if HAVE_NATIVE_mpn_addmul_2 || HAVE_NATIVE_mpn_redc_2
  1291. #define WANT_REDC_2 1
  1292. #endif
  1293. #if WANT_REDC_2
  1294.   {
  1295.     static struct param_t  param;
  1296.     param.name = "REDC_1_TO_REDC_2_THRESHOLD";
  1297.     param.function = speed_mpn_redc_1;
  1298.     param.function2 = speed_mpn_redc_2;
  1299.     param.max_size = TUNE_REDC_2_MAX;
  1300.     param.noprint = 1;
  1301.     one (&redc_1_to_redc_2_threshold, &param);
  1302.   }
  1303.   {
  1304.     static struct param_t  param;
  1305.     param.name = "REDC_2_TO_REDC_N_THRESHOLD";
  1306.     param.function = speed_mpn_redc_2;
  1307.     param.function2 = speed_mpn_redc_n;
  1308.     param.min_size = 16;
  1309.     param.noprint = 1;
  1310.     one (&redc_2_to_redc_n_threshold, &param);
  1311.   }
  1312.   if (redc_1_to_redc_2_threshold >= TUNE_REDC_2_MAX - 1)
  1313.     {
  1314.       /* Disable REDC_2.  This is not supposed to happen.  */
  1315.       print_define ("REDC_1_TO_REDC_2_THRESHOLD", REDC_2_TO_REDC_N_THRESHOLD);
  1316.       print_define_remark ("REDC_2_TO_REDC_N_THRESHOLD", 0, "anomaly: never REDC_2");
  1317.     }
  1318.   else
  1319.     {
  1320.       print_define ("REDC_1_TO_REDC_2_THRESHOLD", REDC_1_TO_REDC_2_THRESHOLD);
  1321.       print_define ("REDC_2_TO_REDC_N_THRESHOLD", REDC_2_TO_REDC_N_THRESHOLD);
  1322.     }
  1323. #else
  1324.   {
  1325.     static struct param_t  param;
  1326.     param.name = "REDC_1_TO_REDC_N_THRESHOLD";
  1327.     param.function = speed_mpn_redc_1;
  1328.     param.function2 = speed_mpn_redc_n;
  1329.     param.min_size = 16;
  1330.     one (&redc_1_to_redc_n_threshold, &param);
  1331.   }
  1332. #endif
  1333. }
  1334. void
  1335. tune_matrix22_mul (void)
  1336. {
  1337.   static struct param_t  param;
  1338.   param.name = "MATRIX22_STRASSEN_THRESHOLD";
  1339.   param.function = speed_mpn_matrix22_mul;
  1340.   param.min_size = 2;
  1341.   one (&matrix22_strassen_threshold, &param);
  1342. }
  1343. void
  1344. tune_hgcd (void)
  1345. {
  1346.   static struct param_t  param;
  1347.   param.name = "HGCD_THRESHOLD";
  1348.   param.function = speed_mpn_hgcd;
  1349.   /* We seem to get strange results for small sizes */
  1350.   param.min_size = 30;
  1351.   one (&hgcd_threshold, &param);
  1352. }
  1353. void
  1354. tune_gcd_dc (void)
  1355. {
  1356.   static struct param_t  param;
  1357.   param.name = "GCD_DC_THRESHOLD";
  1358.   param.function = speed_mpn_gcd;
  1359.   param.min_size = hgcd_threshold;
  1360.   param.max_size = 3000;
  1361.   param.step_factor = 0.02;
  1362.   one (&gcd_dc_threshold, &param);
  1363. }
  1364. void
  1365. tune_gcdext_dc (void)
  1366. {
  1367.   static struct param_t  param;
  1368.   param.name = "GCDEXT_DC_THRESHOLD";
  1369.   param.function = speed_mpn_gcdext;
  1370.   param.min_size = hgcd_threshold;
  1371.   param.max_size = 3000;
  1372.   param.step_factor = 0.02;
  1373.   one (&gcdext_dc_threshold, &param);
  1374. }
  1375. /* size_extra==1 reflects the fact that with high<divisor one division is
  1376.    always skipped.  Forcing high<divisor while testing ensures consistency
  1377.    while stepping through sizes, ie. that size-1 divides will be done each
  1378.    time.
  1379.    min_size==2 and min_is_always are used so that if plain division is only
  1380.    better at size==1 then don't bother including that code just for that
  1381.    case, instead go with preinv always and get a size saving.  */
  1382. #define DIV_1_PARAMS                    
  1383.   param.check_size = 256;               
  1384.   param.min_size = 2;                   
  1385.   param.min_is_always = 1;              
  1386.   param.data_high = DATA_HIGH_LT_R;     
  1387.   param.size_extra = 1;                 
  1388.   param.stop_factor = 2.0;
  1389. double (*tuned_speed_mpn_divrem_1) __GMP_PROTO ((struct speed_params *));
  1390. void
  1391. tune_divrem_1 (void)
  1392. {
  1393.   /* plain version by default */
  1394.   tuned_speed_mpn_divrem_1 = speed_mpn_divrem_1;
  1395.   /* No support for tuning native assembler code, do that by hand and put
  1396.      the results in the .asm file, there's no need for such thresholds to
  1397.      appear in gmp-mparam.h.  */
  1398.   if (HAVE_NATIVE_mpn_divrem_1)
  1399.     return;
  1400.   if (GMP_NAIL_BITS != 0)
  1401.     {
  1402.       print_define_remark ("DIVREM_1_NORM_THRESHOLD", MP_SIZE_T_MAX,
  1403.                            "no preinv with nails");
  1404.       print_define_remark ("DIVREM_1_UNNORM_THRESHOLD", MP_SIZE_T_MAX,
  1405.                            "no preinv with nails");
  1406.       return;
  1407.     }
  1408.   if (UDIV_PREINV_ALWAYS)
  1409.     {
  1410.       print_define_remark ("DIVREM_1_NORM_THRESHOLD", 0L, "preinv always");
  1411.       print_define ("DIVREM_1_UNNORM_THRESHOLD", 0L);
  1412.       return;
  1413.     }
  1414.   tuned_speed_mpn_divrem_1 = speed_mpn_divrem_1_tune;
  1415.   /* Tune for the integer part of mpn_divrem_1.  This will very possibly be
  1416.      a bit out for the fractional part, but that's too bad, the integer part
  1417.      is more important. */
  1418.   {
  1419.     static struct param_t  param;
  1420.     param.name = "DIVREM_1_NORM_THRESHOLD";
  1421.     DIV_1_PARAMS;
  1422.     s.r = randlimb_norm ();
  1423.     param.function = speed_mpn_divrem_1_tune;
  1424.     one (&divrem_1_norm_threshold, &param);
  1425.   }
  1426.   {
  1427.     static struct param_t  param;
  1428.     param.name = "DIVREM_1_UNNORM_THRESHOLD";
  1429.     DIV_1_PARAMS;
  1430.     s.r = randlimb_half ();
  1431.     param.function = speed_mpn_divrem_1_tune;
  1432.     one (&divrem_1_unnorm_threshold, &param);
  1433.   }
  1434. }
  1435. void
  1436. tune_mod_1 (void)
  1437. {
  1438.   /* No support for tuning native assembler code, do that by hand and put
  1439.      the results in the .asm file, there's no need for such thresholds to
  1440.      appear in gmp-mparam.h.  */
  1441.   if (HAVE_NATIVE_mpn_mod_1)
  1442.     return;
  1443.   if (GMP_NAIL_BITS != 0)
  1444.     {
  1445.       print_define_remark ("MOD_1_NORM_THRESHOLD", MP_SIZE_T_MAX,
  1446.                            "no preinv with nails");
  1447.       print_define_remark ("MOD_1_UNNORM_THRESHOLD", MP_SIZE_T_MAX,
  1448.                            "no preinv with nails");
  1449.       return;
  1450.     }
  1451.   if (UDIV_PREINV_ALWAYS)
  1452.     {
  1453.       print_define ("MOD_1_NORM_THRESHOLD", 0L);
  1454.       print_define ("MOD_1_UNNORM_THRESHOLD", 0L);
  1455.     }
  1456.   else
  1457.     {
  1458.       {
  1459. static struct param_t  param;
  1460. param.name = "MOD_1_NORM_THRESHOLD";
  1461. DIV_1_PARAMS;
  1462. s.r = randlimb_norm ();
  1463. param.function = speed_mpn_mod_1_tune;
  1464. one (&mod_1_norm_threshold, &param);
  1465.       }
  1466.       {
  1467. static struct param_t  param;
  1468. param.name = "MOD_1_UNNORM_THRESHOLD";
  1469. DIV_1_PARAMS;
  1470. s.r = randlimb_half ();
  1471. param.function = speed_mpn_mod_1_tune;
  1472. one (&mod_1_unnorm_threshold, &param);
  1473.       }
  1474.     }
  1475.   {
  1476.     static struct param_t  param;
  1477.     param.check_size = 256;
  1478.     s.r = randlimb_norm ();
  1479.     param.function = speed_mpn_mod_1_tune;
  1480.     param.name = "MOD_1N_TO_MOD_1_1_THRESHOLD";
  1481.     param.min_size = 2;
  1482.     one (&mod_1n_to_mod_1_1_threshold, &param);
  1483.   }
  1484.   {
  1485.     static struct param_t  param;
  1486.     param.check_size = 256;
  1487.     s.r = randlimb_norm () / 5;
  1488.     param.function = speed_mpn_mod_1_tune;
  1489.     param.noprint = 1;
  1490.     param.name = "MOD_1U_TO_MOD_1_1_THRESHOLD";
  1491.     param.min_size = 2;
  1492.     one (&mod_1u_to_mod_1_1_threshold, &param);
  1493.     param.name = "MOD_1_1_TO_MOD_1_2_THRESHOLD";
  1494.     param.min_size = mod_1u_to_mod_1_1_threshold;
  1495.     one (&mod_1_1_to_mod_1_2_threshold, &param);
  1496.     if (mod_1u_to_mod_1_1_threshold + 2 >= mod_1_1_to_mod_1_2_threshold)
  1497.       {
  1498. /* Disable mod_1_1, mod_1_2 is always faster.  Measure when to switch
  1499.    (from mod_1_unnorm) to mod_1_2.  */
  1500. mod_1_1_to_mod_1_2_threshold = 0;
  1501. /* This really measures mod_1u -> mod_1_2 */
  1502. param.min_size = 1;
  1503. one (&mod_1u_to_mod_1_1_threshold, &param);
  1504.       }
  1505.     print_define_remark ("MOD_1U_TO_MOD_1_1_THRESHOLD", mod_1u_to_mod_1_1_threshold, NULL);
  1506.     param.name = "MOD_1_2_TO_MOD_1_4_THRESHOLD";
  1507.     param.min_size = mod_1_1_to_mod_1_2_threshold;
  1508.     one (&mod_1_2_to_mod_1_4_threshold, &param);
  1509.     if (mod_1_1_to_mod_1_2_threshold + 2 >= mod_1_2_to_mod_1_4_threshold)
  1510.       {
  1511. /* Disable mod_1_2, mod_1_4 is always faster.  Measure when to switch
  1512.    (from mod_1_unnorm or mod_1_1) to mod_1_4.  */
  1513. mod_1_2_to_mod_1_4_threshold = 0;
  1514. param.min_size = 1;
  1515. one (&mod_1_1_to_mod_1_2_threshold, &param);
  1516.       }
  1517.     print_define_remark ("MOD_1_1_TO_MOD_1_2_THRESHOLD", mod_1_1_to_mod_1_2_threshold, NULL);
  1518.     print_define_remark ("MOD_1_2_TO_MOD_1_4_THRESHOLD", mod_1_2_to_mod_1_4_threshold, NULL);
  1519.   }
  1520.   {
  1521.     static struct param_t  param;
  1522.     param.check_size = 256;
  1523.     param.name = "PREINV_MOD_1_TO_MOD_1_THRESHOLD";
  1524.     s.r = randlimb_norm ();
  1525.     param.function = speed_mpn_preinv_mod_1;
  1526.     param.function2 = speed_mpn_mod_1_tune;
  1527.     one (&preinv_mod_1_to_mod_1_threshold, &param);
  1528.   }
  1529. }
  1530. /* A non-zero DIVREM_1_UNNORM_THRESHOLD (or DIVREM_1_NORM_THRESHOLD) would
  1531.    imply that udiv_qrnnd_preinv is worth using, but it seems most
  1532.    straightforward to compare mpn_preinv_divrem_1 and mpn_divrem_1_div
  1533.    directly.  */
  1534. void
  1535. tune_preinv_divrem_1 (void)
  1536. {
  1537.   static struct param_t  param;
  1538.   speed_function_t  divrem_1;
  1539.   const char        *divrem_1_name;
  1540.   double            t1, t2;
  1541.   if (GMP_NAIL_BITS != 0)
  1542.     {
  1543.       print_define_remark ("USE_PREINV_DIVREM_1", 0, "no preinv with nails");
  1544.       return;
  1545.     }
  1546.   /* Any native version of mpn_preinv_divrem_1 is assumed to exist because
  1547.      it's faster than mpn_divrem_1.  */
  1548.   if (HAVE_NATIVE_mpn_preinv_divrem_1)
  1549.     {
  1550.       print_define_remark ("USE_PREINV_DIVREM_1", 1, "native");
  1551.       return;
  1552.     }
  1553.   /* If udiv_qrnnd_preinv is the only division method then of course
  1554.      mpn_preinv_divrem_1 should be used.  */
  1555.   if (UDIV_PREINV_ALWAYS)
  1556.     {
  1557.       print_define_remark ("USE_PREINV_DIVREM_1", 1, "preinv always");
  1558.       return;
  1559.     }
  1560.   /* If we've got an assembler version of mpn_divrem_1, then compare against
  1561.      that, not the mpn_divrem_1_div generic C.  */
  1562.   if (HAVE_NATIVE_mpn_divrem_1)
  1563.     {
  1564.       divrem_1 = speed_mpn_divrem_1;
  1565.       divrem_1_name = "mpn_divrem_1";
  1566.     }
  1567.   else
  1568.     {
  1569.       divrem_1 = speed_mpn_divrem_1_div;
  1570.       divrem_1_name = "mpn_divrem_1_div";
  1571.     }
  1572.   param.data_high = DATA_HIGH_LT_R; /* allow skip one division */
  1573.   s.size = 200;                     /* generous but not too big */
  1574.   /* Divisor, nonzero.  Unnormalized so as to exercise the shift!=0 case,
  1575.      since in general that's probably most common, though in fact for a
  1576.      64-bit limb mp_bases[10].big_base is normalized.  */
  1577.   s.r = urandom() & (GMP_NUMB_MASK >> 4);
  1578.   if (s.r == 0) s.r = 123;
  1579.   t1 = tuneup_measure (speed_mpn_preinv_divrem_1, &param, &s);
  1580.   t2 = tuneup_measure (divrem_1, &param, &s);
  1581.   if (t1 == -1.0 || t2 == -1.0)
  1582.     {
  1583.       printf ("Oops, can't measure mpn_preinv_divrem_1 and %s at %ldn",
  1584.               divrem_1_name, (long) s.size);
  1585.       abort ();
  1586.     }
  1587.   if (option_trace >= 1)
  1588.     printf ("size=%ld, mpn_preinv_divrem_1 %.9f, %s %.9fn",
  1589.             (long) s.size, t1, divrem_1_name, t2);
  1590.   print_define_remark ("USE_PREINV_DIVREM_1", (mp_size_t) (t1 < t2), NULL);
  1591. }
  1592. void
  1593. tune_divrem_2 (void)
  1594. {
  1595.   static struct param_t  param;
  1596.   /* No support for tuning native assembler code, do that by hand and put
  1597.      the results in the .asm file, and there's no need for such thresholds
  1598.      to appear in gmp-mparam.h.  */
  1599.   if (HAVE_NATIVE_mpn_divrem_2)
  1600.     return;
  1601.   if (GMP_NAIL_BITS != 0)
  1602.     {
  1603.       print_define_remark ("DIVREM_2_THRESHOLD", MP_SIZE_T_MAX,
  1604.                            "no preinv with nails");
  1605.       return;
  1606.     }
  1607.   if (UDIV_PREINV_ALWAYS)
  1608.     {
  1609.       print_define_remark ("DIVREM_2_THRESHOLD", 0L, "preinv always");
  1610.       return;
  1611.     }
  1612.   /* Tune for the integer part of mpn_divrem_2.  This will very possibly be
  1613.      a bit out for the fractional part, but that's too bad, the integer part
  1614.      is more important.
  1615.      min_size must be >=2 since nsize>=2 is required, but is set to 4 to save
  1616.      code space if plain division is better only at size==2 or size==3. */
  1617.   param.name = "DIVREM_2_THRESHOLD";
  1618.   param.check_size = 256;
  1619.   param.min_size = 4;
  1620.   param.min_is_always = 1;
  1621.   param.size_extra = 2;      /* does qsize==nsize-2 divisions */
  1622.   param.stop_factor = 2.0;
  1623.   s.r = randlimb_norm ();
  1624.   param.function = speed_mpn_divrem_2;
  1625.   one (&divrem_2_threshold, &param);
  1626. }
  1627. /* mpn_divexact_1 is vaguely expected to be used on smallish divisors, so
  1628.    tune for that.  Its speed can differ on odd or even divisor, so take an
  1629.    average threshold for the two.
  1630.    mpn_divrem_1 can vary with high<divisor or not, whereas mpn_divexact_1
  1631.    might not vary that way, but don't test this since high<divisor isn't
  1632.    expected to occur often with small divisors.  */
  1633. void
  1634. tune_divexact_1 (void)
  1635. {
  1636.   static struct param_t  param;
  1637.   mp_size_t  thresh[2], average;
  1638.   int        low, i;
  1639.   /* Any native mpn_divexact_1 is assumed to incorporate all the speed of a
  1640.      full mpn_divrem_1.  */
  1641.   if (HAVE_NATIVE_mpn_divexact_1)
  1642.     {
  1643.       print_define_remark ("DIVEXACT_1_THRESHOLD", 0, "always (native)");
  1644.       return;
  1645.     }
  1646.   ASSERT_ALWAYS (tuned_speed_mpn_divrem_1 != NULL);
  1647.   param.name = "DIVEXACT_1_THRESHOLD";
  1648.   param.data_high = DATA_HIGH_GE_R;
  1649.   param.check_size = 256;
  1650.   param.min_size = 2;
  1651.   param.stop_factor = 1.5;
  1652.   param.function  = tuned_speed_mpn_divrem_1;
  1653.   param.function2 = speed_mpn_divexact_1;
  1654.   param.noprint = 1;
  1655.   print_define_start (param.name);
  1656.   for (low = 0; low <= 1; low++)
  1657.     {
  1658.       s.r = randlimb_half();
  1659.       if (low == 0)
  1660.         s.r |= 1;
  1661.       else
  1662.         s.r &= ~CNST_LIMB(7);
  1663.       one (&thresh[low], &param);
  1664.       if (option_trace)
  1665.         printf ("low=%d thresh %ldn", low, (long) thresh[low]);
  1666.       if (thresh[low] == MP_SIZE_T_MAX)
  1667.         {
  1668.           average = MP_SIZE_T_MAX;
  1669.           goto divexact_1_done;
  1670.         }
  1671.     }
  1672.   if (option_trace)
  1673.     {
  1674.       printf ("average of:");
  1675.       for (i = 0; i < numberof(thresh); i++)
  1676.         printf (" %ld", (long) thresh[i]);
  1677.       printf ("n");
  1678.     }
  1679.   average = 0;
  1680.   for (i = 0; i < numberof(thresh); i++)
  1681.     average += thresh[i];
  1682.   average /= numberof(thresh);
  1683.   /* If divexact turns out to be better as early as 3 limbs, then use it
  1684.      always, so as to reduce code size and conditional jumps.  */
  1685.   if (average <= 3)
  1686.     average = 0;
  1687.  divexact_1_done:
  1688.   print_define_end (param.name, average);
  1689. }
  1690. /* The generic mpn_modexact_1_odd skips a divide step if high<divisor, the
  1691.    same as mpn_mod_1, but this might not be true of an assembler
  1692.    implementation.  The threshold used is an average based on data where a
  1693.    divide can be skipped and where it can't.
  1694.    If modexact turns out to be better as early as 3 limbs, then use it
  1695.    always, so as to reduce code size and conditional jumps.  */
  1696. void
  1697. tune_modexact_1_odd (void)
  1698. {
  1699.   static struct param_t  param;
  1700.   mp_size_t  thresh_lt, thresh_ge, average;
  1701. #if 0
  1702.   /* Any native mpn_modexact_1_odd is assumed to incorporate all the speed
  1703.      of a full mpn_mod_1.  */
  1704.   if (HAVE_NATIVE_mpn_modexact_1_odd)
  1705.     {
  1706.       print_define_remark ("BMOD_1_TO_MOD_1_THRESHOLD", MP_SIZE_T_MAX, "always bmod_1");
  1707.       return;
  1708.     }
  1709. #endif
  1710.   param.name = "BMOD_1_TO_MOD_1_THRESHOLD";
  1711.   param.check_size = 256;
  1712.   param.min_size = 2;
  1713.   param.stop_factor = 1.5;
  1714.   param.function  = speed_mpn_modexact_1c_odd;
  1715.   param.function2 = speed_mpn_mod_1_tune;
  1716.   param.noprint = 1;
  1717.   s.r = randlimb_half () | 1;
  1718.   print_define_start (param.name);
  1719.   param.data_high = DATA_HIGH_LT_R;
  1720.   one (&thresh_lt, &param);
  1721.   if (option_trace)
  1722.     printf ("lt thresh %ldn", (long) thresh_lt);
  1723.   average = thresh_lt;
  1724.   if (thresh_lt != MP_SIZE_T_MAX)
  1725.     {
  1726.       param.data_high = DATA_HIGH_GE_R;
  1727.       one (&thresh_ge, &param);
  1728.       if (option_trace)
  1729.         printf ("ge thresh %ldn", (long) thresh_ge);
  1730.       if (thresh_ge != MP_SIZE_T_MAX)
  1731.         {
  1732.           average = (thresh_ge + thresh_lt) / 2;
  1733.           if (thresh_ge <= 3)
  1734.             average = 0;
  1735.         }
  1736.     }
  1737.   print_define_end (param.name, average);
  1738. }
  1739. void
  1740. tune_jacobi_base (void)
  1741. {
  1742.   static struct param_t  param;
  1743.   double   t1, t2, t3;
  1744.   int      method;
  1745.   s.size = GMP_LIMB_BITS * 3 / 4;
  1746.   t1 = tuneup_measure (speed_mpn_jacobi_base_1, &param, &s);
  1747.   if (option_trace >= 1)
  1748.     printf ("size=%ld, mpn_jacobi_base_1 %.9fn", (long) s.size, t1);
  1749.   t2 = tuneup_measure (speed_mpn_jacobi_base_2, &param, &s);
  1750.   if (option_trace >= 1)
  1751.     printf ("size=%ld, mpn_jacobi_base_2 %.9fn", (long) s.size, t2);
  1752.   t3 = tuneup_measure (speed_mpn_jacobi_base_3, &param, &s);
  1753.   if (option_trace >= 1)
  1754.     printf ("size=%ld, mpn_jacobi_base_3 %.9fn", (long) s.size, t3);
  1755.   if (t1 == -1.0 || t2 == -1.0 || t3 == -1.0)
  1756.     {
  1757.       printf ("Oops, can't measure all mpn_jacobi_base methods at %ldn",
  1758.               (long) s.size);
  1759.       abort ();
  1760.     }
  1761.   if (t1 < t2 && t1 < t3)
  1762.     method = 1;
  1763.   else if (t2 < t3)
  1764.     method = 2;
  1765.   else
  1766.     method = 3;
  1767.   print_define ("JACOBI_BASE_METHOD", method);
  1768. }
  1769. void
  1770. tune_get_str (void)
  1771. {
  1772.   /* Tune for decimal, it being most common.  Some rough testing suggests
  1773.      other bases are different, but not by very much.  */
  1774.   s.r = 10;
  1775.   {
  1776.     static struct param_t  param;
  1777.     GET_STR_PRECOMPUTE_THRESHOLD = 0;
  1778.     param.name = "GET_STR_DC_THRESHOLD";
  1779.     param.function = speed_mpn_get_str;
  1780.     param.min_size = 4;
  1781.     param.max_size = GET_STR_THRESHOLD_LIMIT;
  1782.     one (&get_str_dc_threshold, &param);
  1783.   }
  1784.   {
  1785.     static struct param_t  param;
  1786.     param.name = "GET_STR_PRECOMPUTE_THRESHOLD";
  1787.     param.function = speed_mpn_get_str;
  1788.     param.min_size = GET_STR_DC_THRESHOLD;
  1789.     param.max_size = GET_STR_THRESHOLD_LIMIT;
  1790.     one (&get_str_precompute_threshold, &param);
  1791.   }
  1792. }
  1793. double
  1794. speed_mpn_pre_set_str (struct speed_params *s)
  1795. {
  1796.   unsigned char *str;
  1797.   mp_ptr     wp;
  1798.   mp_size_t  wn;
  1799.   unsigned   i;
  1800.   int        base;
  1801.   double     t;
  1802.   mp_ptr powtab_mem, tp;
  1803.   powers_t powtab[GMP_LIMB_BITS];
  1804.   mp_size_t un;
  1805.   int chars_per_limb;
  1806.   TMP_DECL;
  1807.   SPEED_RESTRICT_COND (s->size >= 1);
  1808.   base = s->r == 0 ? 10 : s->r;
  1809.   SPEED_RESTRICT_COND (base >= 2 && base <= 256);
  1810.   TMP_MARK;
  1811.   str = TMP_ALLOC (s->size);
  1812.   for (i = 0; i < s->size; i++)
  1813.     str[i] = s->xp[i] % base;
  1814.   wn = ((mp_size_t) (s->size / mp_bases[base].chars_per_bit_exactly))
  1815.     / GMP_LIMB_BITS + 2;
  1816.   SPEED_TMP_ALLOC_LIMBS (wp, wn, s->align_wp);
  1817.   /* use this during development to check wn is big enough */
  1818.   /*
  1819.   ASSERT_ALWAYS (mpn_set_str (wp, str, s->size, base) <= wn);
  1820.   */
  1821.   speed_operand_src (s, (mp_ptr) str, s->size/BYTES_PER_MP_LIMB);
  1822.   speed_operand_dst (s, wp, wn);
  1823.   speed_cache_fill (s);
  1824.   chars_per_limb = mp_bases[base].chars_per_limb;
  1825.   un = s->size / chars_per_limb + 1;
  1826.   powtab_mem = TMP_BALLOC_LIMBS (mpn_dc_set_str_powtab_alloc (un));
  1827.   mpn_set_str_compute_powtab (powtab, powtab_mem, un, base);
  1828.   tp = TMP_BALLOC_LIMBS (mpn_dc_set_str_itch (un));
  1829.   speed_starttime ();
  1830.   i = s->reps;
  1831.   do
  1832.     {
  1833.       mpn_pre_set_str (wp, str, s->size, powtab, tp);
  1834.     }
  1835.   while (--i != 0);
  1836.   t = speed_endtime ();
  1837.   TMP_FREE;
  1838.   return t;
  1839. }
  1840. void
  1841. tune_set_str (void)
  1842. {
  1843.   s.r = 10;  /* decimal */
  1844.   {
  1845.     static struct param_t  param;
  1846.     SET_STR_PRECOMPUTE_THRESHOLD = 0;
  1847.     param.step_factor = 0.01;
  1848.     param.name = "SET_STR_DC_THRESHOLD";
  1849.     param.function = speed_mpn_pre_set_str;
  1850.     param.min_size = 100;
  1851.     param.max_size = 50000;
  1852.     one (&set_str_dc_threshold, &param);
  1853.   }
  1854.   {
  1855.     static struct param_t  param;
  1856.     param.step_factor = 0.02;
  1857.     param.name = "SET_STR_PRECOMPUTE_THRESHOLD";
  1858.     param.function = speed_mpn_set_str;
  1859.     param.min_size = SET_STR_DC_THRESHOLD;
  1860.     param.max_size = 100000;
  1861.     one (&set_str_precompute_threshold, &param);
  1862.   }
  1863. }
  1864. void
  1865. tune_fft_mul (void)
  1866. {
  1867.   static struct fft_param_t  param;
  1868.   if (option_fft_max_size == 0)
  1869.     return;
  1870.   param.table_name          = "MUL_FFT_TABLE3";
  1871.   param.threshold_name      = "MUL_FFT_THRESHOLD";
  1872.   param.p_threshold         = &mul_fft_threshold;
  1873.   param.modf_threshold_name = "MUL_FFT_MODF_THRESHOLD";
  1874.   param.p_modf_threshold    = &mul_fft_modf_threshold;
  1875.   param.first_size          = MUL_TOOM33_THRESHOLD / 2;
  1876.   param.max_size            = option_fft_max_size;
  1877.   param.function            = speed_mpn_fft_mul;
  1878.   param.mul_modf_function   = speed_mpn_mul_fft;
  1879.   param.mul_function        = speed_mpn_mul_n;
  1880.   param.sqr = 0;
  1881.   fft (&param);
  1882. }
  1883. void
  1884. tune_fft_sqr (void)
  1885. {
  1886.   static struct fft_param_t  param;
  1887.   if (option_fft_max_size == 0)
  1888.     return;
  1889.   param.table_name          = "SQR_FFT_TABLE3";
  1890.   param.threshold_name      = "SQR_FFT_THRESHOLD";
  1891.   param.p_threshold         = &sqr_fft_threshold;
  1892.   param.modf_threshold_name = "SQR_FFT_MODF_THRESHOLD";
  1893.   param.p_modf_threshold    = &sqr_fft_modf_threshold;
  1894.   param.first_size          = SQR_TOOM3_THRESHOLD / 2;
  1895.   param.max_size            = option_fft_max_size;
  1896.   param.function            = speed_mpn_fft_sqr;
  1897.   param.mul_modf_function   = speed_mpn_mul_fft_sqr;
  1898.   param.mul_function        = speed_mpn_sqr;
  1899.   param.sqr = 1;
  1900.   fft (&param);
  1901. }
  1902. void
  1903. all (void)
  1904. {
  1905.   time_t  start_time, end_time;
  1906.   TMP_DECL;
  1907.   TMP_MARK;
  1908.   SPEED_TMP_ALLOC_LIMBS (s.xp_block, SPEED_BLOCK_SIZE, 0);
  1909.   SPEED_TMP_ALLOC_LIMBS (s.yp_block, SPEED_BLOCK_SIZE, 0);
  1910.   mpn_random (s.xp_block, SPEED_BLOCK_SIZE);
  1911.   mpn_random (s.yp_block, SPEED_BLOCK_SIZE);
  1912.   fprintf (stderr, "Parameters for %sn", GMP_MPARAM_H_SUGGEST);
  1913.   speed_time_init ();
  1914.   fprintf (stderr, "Using: %sn", speed_time_string);
  1915.   fprintf (stderr, "speed_precision %d", speed_precision);
  1916.   if (speed_unittime == 1.0)
  1917.     fprintf (stderr, ", speed_unittime 1 cycle");
  1918.   else
  1919.     fprintf (stderr, ", speed_unittime %.2e secs", speed_unittime);
  1920.   if (speed_cycletime == 1.0 || speed_cycletime == 0.0)
  1921.     fprintf (stderr, ", CPU freq unknownn");
  1922.   else
  1923.     fprintf (stderr, ", CPU freq %.2f MHzn", 1e-6/speed_cycletime);
  1924.   fprintf (stderr, "DEFAULT_MAX_SIZE %d, fft_max_size %ldn",
  1925.            DEFAULT_MAX_SIZE, (long) option_fft_max_size);
  1926.   fprintf (stderr, "n");
  1927.   time (&start_time);
  1928.   {
  1929.     struct tm  *tp;
  1930.     tp = localtime (&start_time);
  1931.     printf ("/* Generated by tuneup.c, %d-%02d-%02d, ",
  1932.             tp->tm_year+1900, tp->tm_mon+1, tp->tm_mday);
  1933. #ifdef __GNUC__
  1934.     /* gcc sub-minor version doesn't seem to come through as a define */
  1935.     printf ("gcc %d.%d */n", __GNUC__, __GNUC_MINOR__);
  1936. #define PRINTED_COMPILER
  1937. #endif
  1938. #if defined (__SUNPRO_C)
  1939.     printf ("Sun C %d.%d */n", __SUNPRO_C / 0x100, __SUNPRO_C % 0x100);
  1940. #define PRINTED_COMPILER
  1941. #endif
  1942. #if ! defined (__GNUC__) && defined (__sgi) && defined (_COMPILER_VERSION)
  1943.     /* gcc defines __sgi and _COMPILER_VERSION on irix 6, avoid that */
  1944.     printf ("MIPSpro C %d.%d.%d */n",
  1945.     _COMPILER_VERSION / 100,
  1946.     _COMPILER_VERSION / 10 % 10,
  1947.     _COMPILER_VERSION % 10);
  1948. #define PRINTED_COMPILER
  1949. #endif
  1950. #if defined (__DECC) && defined (__DECC_VER)
  1951.     printf ("DEC C %d */n", __DECC_VER);
  1952. #define PRINTED_COMPILER
  1953. #endif
  1954. #if ! defined (PRINTED_COMPILER)
  1955.     printf ("system compiler */n");
  1956. #endif
  1957.   }
  1958.   printf ("n");
  1959.   tune_divrem_1 ();
  1960.   tune_mod_1 ();
  1961.   tune_preinv_divrem_1 ();
  1962.   tune_divrem_2 ();
  1963.   tune_divexact_1 ();
  1964.   tune_modexact_1_odd ();
  1965.   printf("n");
  1966.   tune_mul_n ();
  1967.   printf("n");
  1968.   tune_mul ();
  1969.   printf("n");
  1970.   tune_sqr ();
  1971.   printf("n");
  1972.   tune_mulmod_bnm1 ();
  1973.   tune_sqrmod_bnm1 ();
  1974.   printf("n");
  1975.   tune_fft_mul ();
  1976.   printf("n");
  1977.   tune_fft_sqr ();
  1978.   printf ("n");
  1979.   tune_mullo ();
  1980.   printf("n");
  1981.   tune_dc_div ();
  1982.   tune_dc_bdiv ();
  1983.   printf("n");
  1984.   tune_invertappr ();
  1985.   tune_invert ();
  1986.   printf("n");
  1987.   tune_binvert ();
  1988.   tune_redc ();
  1989.   printf("n");
  1990.   tune_mu_div ();
  1991.   tune_mu_bdiv ();
  1992.   printf("n");
  1993.   tune_matrix22_mul ();
  1994.   tune_hgcd ();
  1995.   tune_gcd_dc ();
  1996.   tune_gcdext_dc ();
  1997.   tune_jacobi_base ();
  1998.   printf("n");
  1999.   tune_get_str ();
  2000.   tune_set_str ();
  2001.   printf("n");
  2002.   time (&end_time);
  2003.   printf ("/* Tuneup completed successfully, took %ld seconds */n",
  2004.           (long) (end_time - start_time));
  2005.   TMP_FREE;
  2006. }
  2007. int
  2008. main (int argc, char *argv[])
  2009. {
  2010.   int  opt;
  2011.   /* Unbuffered so if output is redirected to a file it isn't lost if the
  2012.      program is killed part way through.  */
  2013.   setbuf (stdout, NULL);
  2014.   setbuf (stderr, NULL);
  2015.   while ((opt = getopt(argc, argv, "f:o:p:t")) != EOF)
  2016.     {
  2017.       switch (opt) {
  2018.       case 'f':
  2019.         if (optarg[0] == 't')
  2020.           option_fft_trace = 2;
  2021.         else
  2022.           option_fft_max_size = atol (optarg);
  2023.         break;
  2024.       case 'o':
  2025.         speed_option_set (optarg);
  2026.         break;
  2027.       case 'p':
  2028.         speed_precision = atoi (optarg);
  2029.         break;
  2030.       case 't':
  2031.         option_trace++;
  2032.         break;
  2033.       case '?':
  2034.         exit(1);
  2035.       }
  2036.     }
  2037.   all ();
  2038.   exit (0);
  2039. }