/* [Description]
Copyright 2010-2017 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of 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.
The GNU MPFR 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 Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
#include <stdlib.h>
#include <time.h>
#define MPFR_NEED_LONGLONG_H
#include "mpfr-impl.h"
#undef _PROTO
#define _PROTO __GMP_PROTO
#include "speed.h"
/* Let f be a function for which we have several implementations f1, f2... */
/* We wish to have a quick overview of which implementation is the best */
/* in function of the point x where f(x) is computed and of the prectision */
/* prec requested by the user. */
/* This is performed by drawing a 2D graphic with color indicating which */
/* method is the best. */
/* For building this graphic, the following structur is used (see the core */
/* of generate_2D_sample for an explanation of each field. */
struct speed_params2D
{
/* x-window: [min_x, max_x] or [2^min_x, 2^max_x] */
/* or [-2^(max_x), -2^(min_x)] U [2^min_x, 2^max_x] */
/* depending on the value of logarithmic_scale_x */
double min_x;
double max_x;
/* prec-window: [min_prec, max_prec] */
mpfr_prec_t min_prec;
mpfr_prec_t max_prec;
/* number of sample points for the x-axis and the prec-axis */
int nb_points_x;
int nb_points_prec;
/* should the sample points be logarithmically scaled or not */
int logarithmic_scale_x;
int logarithmic_scale_prec;
/* list of functions g1, g2... measuring the execution time of f1, f2... */
/* when given a point x and a precision prec stored in s. */
/* We use s->xp to store the significant of x, s->r to store its exponent */
/* s->align_xp to store its sign, and s->size to store prec. */
double (**speed_funcs) (struct speed_params *s);
};
/* Given an array t of nb_functions double indicating the timings of several */
/* implementations, return i, such that t[i] is the best timing. */
int
find_best (double *t, int nb_functions)
{
int i, ibest;
double best;
if (nb_functions<=0)
{
fprintf (stderr, "There is no function\n");
abort ();
}
ibest = 0;
best = t[0];
for (i=1; i<nb_functions; i++)
{
if (t[i]<best)
{
best = t[i];
ibest = i;
}
}
return ibest;
}
void generate_2D_sample (FILE *output, struct speed_params2D param)
{
mpfr_t temp;
double incr_prec;
mpfr_t incr_x;
mpfr_t x, x2;
double prec;
struct speed_params s;
int i;
int test;
int nb_functions;
double *t; /* store the timing of each implementation */
/* We first determine how many implementations we have */
nb_functions = 0;
while (param.speed_funcs[nb_functions] != NULL)
nb_functions++;
t = malloc (nb_functions * sizeof (double));
if (t == NULL)
{
fprintf (stderr, "Can't allocate memory.\n");
abort ();
}
mpfr_init2 (temp, MPFR_SMALL_PRECISION);
/* The precision is sampled from min_prec to max_prec with */
/* approximately nb_points_prec points. If logarithmic_scale_prec */
/* is true, the precision is multiplied by incr_prec at each */
/* step. Otherwise, incr_prec is added at each step. */
if (param.logarithmic_scale_prec)
{
mpfr_set_ui (temp, (unsigned long int)param.max_prec, MPFR_RNDU);
mpfr_div_ui (temp, temp, (unsigned long int)param.min_prec, MPFR_RNDU);
mpfr_root (temp, temp,
(unsigned long int)param.nb_points_prec, MPFR_RNDU);
incr_prec = mpfr_get_d (temp, MPFR_RNDU);
}
else
{
incr_prec = (double)param.max_prec - (double)param.min_prec;
incr_prec = incr_prec/((double)param.nb_points_prec);
}
/* The points x are sampled according to the following rule: */
/* If logarithmic_scale_x = 0: */
/* nb_points_x points are equally distributed between min_x and max_x */
/* If logarithmic_scale_x = 1: */
/* nb_points_x points are sampled from 2^(min_x) to 2^(max_x). At */
/* each step, the current point is multiplied by incr_x. */
/* If logarithmic_scale_x = -1: */
/* nb_points_x/2 points are sampled from -2^(max_x) to -2^(min_x) */
/* (at each step, the current point is divided by incr_x); and */
/* nb_points_x/2 points are sampled from 2^(min_x) to 2^(max_x) */
/* (at each step, the current point is multiplied by incr_x). */
mpfr_init2 (incr_x, param.max_prec);
if (param.logarithmic_scale_x == 0)
{
mpfr_set_d (incr_x,
(param.max_x - param.min_x)/(double)param.nb_points_x,
MPFR_RNDU);
}
else if (param.logarithmic_scale_x == -1)
{
mpfr_set_d (incr_x,
2.*(param.max_x - param.min_x)/(double)param.nb_points_x,
MPFR_RNDU);
mpfr_exp2 (incr_x, incr_x, MPFR_RNDU);
}
else
{ /* other values of param.logarithmic_scale_x are considered as 1 */
mpfr_set_d (incr_x,
(param.max_x - param.min_x)/(double)param.nb_points_x,
MPFR_RNDU);
mpfr_exp2 (incr_x, incr_x, MPFR_RNDU);
}
/* Main loop */
mpfr_init2 (x, param.max_prec);
mpfr_init2 (x2, param.max_prec);
prec = (double)param.min_prec;
while (prec <= param.max_prec)
{
printf ("prec = %d\n", (int)prec);
if (param.logarithmic_scale_x == 0)
mpfr_set_d (temp, param.min_x, MPFR_RNDU);
else if (param.logarithmic_scale_x == -1)
{
mpfr_set_d (temp, param.max_x, MPFR_RNDD);
mpfr_exp2 (temp, temp, MPFR_RNDD);
mpfr_neg (temp, temp, MPFR_RNDU);
}
else
{
mpfr_set_d (temp, param.min_x, MPFR_RNDD);
mpfr_exp2 (temp, temp, MPFR_RNDD);
}
/* We perturb x a little bit, in order to avoid trailing zeros that */
/* might change the behavior of algorithms. */
mpfr_const_pi (x, MPFR_RNDN);
mpfr_div_2ui (x, x, 7, MPFR_RNDN);
mpfr_add_ui (x, x, 1, MPFR_RNDN);
mpfr_mul (x, x, temp, MPFR_RNDN);
test = 1;
while (test)
{
mpfr_fprintf (output, "%e\t", mpfr_get_d (x, MPFR_RNDN));
mpfr_fprintf (output, "%Pu\t", (mpfr_prec_t)prec);
s.r = (mp_limb_t)mpfr_get_exp (x);
s.size = (mpfr_prec_t)prec;
s.align_xp = (mpfr_sgn (x) > 0)?1:2;
mpfr_set_prec (x2, (mpfr_prec_t)prec);
mpfr_set (x2, x, MPFR_RNDU);
s.xp = x2->_mpfr_d;
for (i=0; i<nb_functions; i++)
{
t[i] = speed_measure (param.speed_funcs[i], &s);
mpfr_fprintf (output, "%e\t", t[i]);
}
fprintf (output, "%d\n", 1 + find_best (t, nb_functions));
if (param.logarithmic_scale_x == 0)
{
mpfr_add (x, x, incr_x, MPFR_RNDU);
if (mpfr_cmp_d (x, param.max_x) > 0)
test=0;
}
else
{
if (mpfr_sgn (x) < 0 )
{ /* if x<0, it means that logarithmic_scale_x=-1 */
mpfr_div (x, x, incr_x, MPFR_RNDU);
mpfr_abs (temp, x, MPFR_RNDD);
mpfr_log2 (temp, temp, MPFR_RNDD);
if (mpfr_cmp_d (temp, param.min_x) < 0)
mpfr_neg (x, x, MPFR_RNDN);
}
else
{
mpfr_mul (x, x, incr_x, MPFR_RNDU);
mpfr_set (temp, x, MPFR_RNDD);
mpfr_log2 (temp, temp, MPFR_RNDD);
if (mpfr_cmp_d (temp, param.max_x) > 0)
test=0;
}
}
}
prec = ( (param.logarithmic_scale_prec) ? (prec * incr_prec)
: (prec + incr_prec) );
fprintf (output, "\n");
}
free (t);
mpfr_clear (incr_x);
mpfr_clear (x);
mpfr_clear (x2);
mpfr_clear (temp);
return;
}
#define SPEED_MPFR_FUNC_2D(mean_func) \
do \
{ \
double t; \
unsigned i; \
mpfr_t w, x; \
mp_size_t size; \
\
SPEED_RESTRICT_COND (s->size >= MPFR_PREC_MIN); \
SPEED_RESTRICT_COND (s->size <= MPFR_PREC_MAX); \
\
size = (s->size-1)/GMP_NUMB_BITS+1; \
s->xp[size-1] |= MPFR_LIMB_HIGHBIT; \
MPFR_TMP_INIT1 (s->xp, x, s->size); \
MPFR_SET_EXP (x, (mpfr_exp_t) s->r); \
if (s->align_xp == 2) MPFR_SET_NEG (x); \
\
mpfr_init2 (w, s->size); \
speed_starttime (); \
i = s->reps; \
\
do \
mean_func (w, x, MPFR_RNDN); \
while (--i != 0); \
t = speed_endtime (); \
\
mpfr_clear (w); \
return t; \
} \
while (0)
mpfr_prec_t mpfr_exp_2_threshold;
mpfr_prec_t old_threshold = MPFR_EXP_2_THRESHOLD;
#undef MPFR_EXP_2_THRESHOLD
#define MPFR_EXP_2_THRESHOLD mpfr_exp_2_threshold
#include "exp_2.c"
double
timing_exp1 (struct speed_params *s)
{
mpfr_exp_2_threshold = s->size+1;
SPEED_MPFR_FUNC_2D (mpfr_exp_2);
}
double
timing_exp2 (struct speed_params *s)
{
mpfr_exp_2_threshold = s->size-1;
SPEED_MPFR_FUNC_2D (mpfr_exp_2);
}
double
timing_exp3 (struct speed_params *s)
{
SPEED_MPFR_FUNC_2D (mpfr_exp_3);
}
#include "ai.c"
double
timing_ai1 (struct speed_params *s)
{
SPEED_MPFR_FUNC_2D (mpfr_ai1);
}
double
timing_ai2 (struct speed_params *s)
{
SPEED_MPFR_FUNC_2D (mpfr_ai2);
}
/* These functions are for testing purpose only */
/* They are used to draw which method is actually used */
double
virtual_timing_ai1 (struct speed_params *s)
{
double t;
unsigned i;
mpfr_t w, x;
mp_size_t size;
mpfr_t temp1, temp2;
SPEED_RESTRICT_COND (s->size >= MPFR_PREC_MIN);
SPEED_RESTRICT_COND (s->size <= MPFR_PREC_MAX);
size = (s->size-1)/GMP_NUMB_BITS+1;
s->xp[size-1] |= MPFR_LIMB_HIGHBIT;
MPFR_TMP_INIT1 (s->xp, x, s->size);
MPFR_SET_EXP (x, (mpfr_exp_t) s->r);
if (s->align_xp == 2) MPFR_SET_NEG (x);
mpfr_init2 (w, s->size);
speed_starttime ();
i = s->reps;
mpfr_init2 (temp1, MPFR_SMALL_PRECISION);
mpfr_init2 (temp2, MPFR_SMALL_PRECISION);
mpfr_set (temp1, x, MPFR_SMALL_PRECISION);
mpfr_set_si (temp2, MPFR_AI_THRESHOLD2, MPFR_RNDN);
mpfr_mul_ui (temp2, temp2, (unsigned int)MPFR_PREC (w), MPFR_RNDN);
if (MPFR_IS_NEG (x))
mpfr_mul_si (temp1, temp1, MPFR_AI_THRESHOLD1, MPFR_RNDN);
else
mpfr_mul_si (temp1, temp1, MPFR_AI_THRESHOLD3, MPFR_RNDN);
mpfr_add (temp1, temp1, temp2, MPFR_RNDN);
if (mpfr_cmp_si (temp1, MPFR_AI_SCALE) > 0)
t = 1000.;
else
t = 1.;
mpfr_clear (temp1);
mpfr_clear (temp2);
return t;
}
double
virtual_timing_ai2 (struct speed_params *s)
{
double t;
unsigned i;
mpfr_t w, x;
mp_size_t size;
mpfr_t temp1, temp2;
SPEED_RESTRICT_COND (s->size >= MPFR_PREC_MIN);
SPEED_RESTRICT_COND (s->size <= MPFR_PREC_MAX);
size = (s->size-1)/GMP_NUMB_BITS+1;
s->xp[size-1] |= MPFR_LIMB_HIGHBIT;
MPFR_TMP_INIT1 (s->xp, x, s->size);
MPFR_SET_EXP (x, (mpfr_exp_t) s->r);
if (s->align_xp == 2) MPFR_SET_NEG (x);
mpfr_init2 (w, s->size);
speed_starttime ();
i = s->reps;
mpfr_init2 (temp1, MPFR_SMALL_PRECISION);
mpfr_init2 (temp2, MPFR_SMALL_PRECISION);
mpfr_set (temp1, x, MPFR_SMALL_PRECISION);
mpfr_set_si (temp2, MPFR_AI_THRESHOLD2, MPFR_RNDN);
mpfr_mul_ui (temp2, temp2, (unsigned int)MPFR_PREC (w), MPFR_RNDN);
if (MPFR_IS_NEG (x))
mpfr_mul_si (temp1, temp1, MPFR_AI_THRESHOLD1, MPFR_RNDN);
else
mpfr_mul_si (temp1, temp1, MPFR_AI_THRESHOLD3, MPFR_RNDN);
mpfr_add (temp1, temp1, temp2, MPFR_RNDN);
if (mpfr_cmp_si (temp1, MPFR_AI_SCALE) > 0)
t = 1.;
else
t = 1000.;
mpfr_clear (temp1);
mpfr_clear (temp2);
return t;
}
int
main (void)
{
FILE *output;
struct speed_params2D param;
double (*speed_funcs[3]) (struct speed_params *s);
/* char filename[256] = "virtual_timing_ai.dat"; */
/* speed_funcs[0] = virtual_timing_ai1; */
/* speed_funcs[1] = virtual_timing_ai2; */
char filename[256] = "airy.dat";
speed_funcs[0] = timing_ai1;
speed_funcs[1] = timing_ai2;
speed_funcs[2] = NULL;
output = fopen (filename, "w");
if (output == NULL)
{
fprintf (stderr, "Can't open file '%s' for writing.\n", filename);
abort ();
}
param.min_x = -80;
param.max_x = 60;
param.min_prec = 50;
param.max_prec = 1500;
param.nb_points_x = 200;
param.nb_points_prec = 200;
param.logarithmic_scale_x = 0;
param.logarithmic_scale_prec = 0;
param.speed_funcs = speed_funcs;
generate_2D_sample (output, param);
fclose (output);
mpfr_free_cache ();
return 0;
}