/*
This file is part of liberferfc.
Copyright 2009-2010 by
Laboratoire de l'Informatique du Parallélisme, UMR CNRS - ENS Lyon -
UCB Lyon 1 - INRIA 5668,
and by LORIA (CNRS, INPL, INRIA, UHP, U-Nancy 2).
It has been written by S. Chevillard.
Liberferfc 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.
Liberferfc 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 liberferfc. If not, see .
*/
#include
#include
#include
#include
#include "erferfc.h"
#include "timer.h"
#define LOG2EINF 1.4426950408889634 /* 0x3ff71547652b82fe */
int choose_method(mpfr_t x, mpfr_prec_t prec);
/* This program compares the three implementations erf1, erf2, and erfc3 */
/* when used to evaluate erf(x) at precision t, for x and t within two */
/* ranges. For each couple (x,t), the best method is determined and a */
/* a line of the form "x t eq1 eq2 eq3 method chosen_method" */
/* is written in file tune.dat, where eqi is the time (in microseconds) */
/* necessary for evaluating erf(x) by method i, and where */
/* method \in {1,2,3} is the number of the best method, and */
/* chosen_method is the method automatically selected by the algorithm. */
/* x will be chosen logarithmically distributed between 2^{MIN_EXPONENT} */
/* and 2^{MAX_EXPONENT+1} with POINTS_X sample points. */
/* t will be chosen logarithmically distributed between MIN_PRECISION */
/* and MAX_PRECISION with POINTS_PREC sample points. */
/* In order to determine the time spent in one execution of the */
/* function f, we successively run f, during at least */
/* TIME_OF_A_SINGLE_SAMPLE microseconds. Then we divide the total time */
/* by the number of times that f has been executed. */
#define MIN_EXPONENT -10
#define MAX_EXPONENT 5
#define MIN_PRECISION 24
#define MAX_PRECISION 100000
#define POINTS_PREC 240
#define POINTS_X 120
#define TIME_OF_A_SINGLE_SAMPLE 1000.
int initialization(gmp_randstate_t state) {
mpz_t zseed;
const char seed[] = "1BLg5Vq84iukX4SPS2MnBo1dhAwFd8FBWD3kztdhpWsHy3zTUluAsG1y7xOgDwOIFKMSIdKg4eM6yxgUKxpPp51iKnZj2bQAB";
gmp_randinit_default(state);
mpz_init(zseed);
mpz_set_str (zseed, seed, 62);
gmp_randseed(state, zseed);
srand(mpz_get_ui(zseed));
mpz_clear(zseed);
return 0;
}
/* Returns the average time needed to evalute erf(x) with a relative error */
/* smaller than 2^{-prec}, by erf1, expressed in microseconds. */
double timing_eq1(mpfr_t x, mp_prec_t prec) {
timer_type timer;
mpfr_t *y;
int count;
timer_init(timer);
count = 0;
while(timer_get(timer) < TIME_OF_A_SINGLE_SAMPLE) {
timer_start(timer);
y = mp_erf1(x, prec);
timer_stop(timer);
mpfr_clear(*y);
free(y);
count++;
}
return timer_get(timer)/(double)count;
}
/* Returns the average time needed to evalute erf(x) with a relative error */
/* smaller than 2^{-prec}, by erf2, expressed in microseconds. */
double timing_eq2(mpfr_t x, mp_prec_t prec) {
timer_type timer;
mpfr_t *y;
int count;
timer_init(timer);
count = 0;
while(timer_get(timer) < TIME_OF_A_SINGLE_SAMPLE) {
timer_start(timer);
y = mp_erf2(x, prec);
timer_stop(timer);
mpfr_clear(*y);
free(y);
count++;
}
return timer_get(timer)/(double)count;
}
/* Returns the average time needed to evalute erf(x) with a relative error */
/* smaller than 2^{-prec}, by erfc3, expressed in microseconds. */
/* If erfc3 does not provide enough accuracy, -1 is returned. */
double timing_eq3(mpfr_t x, mp_prec_t prec) {
timer_type timer;
mpfr_t *y;
mpfr_t tmp, tmp2;
int count;
mp_prec_t s;
mpfr_init2(tmp, 30);
mpfr_init2(tmp2, 53);
mpfr_set_d(tmp2, LOG2EINF, GMP_RNDD);
mpfr_set(tmp, x, GMP_RNDD);
mpfr_sqr(tmp, tmp, GMP_RNDD);
mpfr_mul(tmp, tmp, tmp2, GMP_RNDD);
mpfr_add_si(tmp, tmp, mpfr_get_exp(x)-3, GMP_RNDD);
mpfr_neg(tmp, tmp, GMP_RNDU);
mpfr_add_ui(tmp, tmp, prec, GMP_RNDU);
mpfr_clear(tmp2);
if (mpfr_cmp_ui(tmp, 1)<=0) {
mpfr_clear(tmp); return 0.;
}
if (!mpfr_fits_ulong_p(tmp, GMP_RNDU)) {
printf("Error precision too high in mp_erfc3\n");
exit(1);
}
else s = mpfr_get_ui(tmp, GMP_RNDU);
timer_init(timer);
count = 0;
while(timer_get(timer) < TIME_OF_A_SINGLE_SAMPLE) {
timer_start(timer);
y = mp_erfc3(x, s);
timer_stop(timer);
if (y==NULL) return -1.;
mpfr_clear(*y);
free(y);
count++;
}
return timer_get(timer)/(double)count;
}
/* Returns the average time needed to evalute erf(x) in precision prec */
/* by cr_erf, expressed in microseconds. */
double timing_cr_erf(mpfr_t x, mp_prec_t prec) {
timer_type timer;
mpfr_t y;
int count;
mpfr_init2(y, prec);
timer_init(timer);
count = 0;
while(timer_get(timer) < TIME_OF_A_SINGLE_SAMPLE) {
timer_start(timer);
cr_erf(y, x, GMP_RNDN);
timer_stop(timer);
count++;
}
mpfr_clear(y);
return timer_get(timer)/(double)count;
}
/* Returns the average time needed to evalute erf(x) in precision prec */
/* by mpfr_erf, expressed in microseconds. */
double timing_mpfr_erf(mpfr_t x, mp_prec_t prec) {
timer_type timer;
mpfr_t y;
int count;
mpfr_init2(y, prec);
timer_init(timer);
count = 0;
while(timer_get(timer) < TIME_OF_A_SINGLE_SAMPLE) {
timer_start(timer);
mpfr_erf(y, x, GMP_RNDN);
timer_stop(timer);
count++;
}
mpfr_clear(y);
return timer_get(timer)/(double)count;
}
int main(void) {
double incr_prec;
mpfr_t incr_x;
double eq1, eq2, eq3;
mpfr_t temp, temp2;
mpfr_t x0, x;
mp_prec_t prec;
gmp_randstate_t state;
mpz_t mant;
int i,j, best;
FILE *output;
initialization(state);
output = fopen("tune.dat", "w");
if (output == NULL) {
printf("Impossible to open tune.dat for writing.\n");
exit(1);
}
/* Chooses x0 randomly with exponent MIN_EXPONENT */
mpfr_init2(x0, MAX_PRECISION);
mpz_init(mant);
mpz_urandomb(mant, state, MAX_PRECISION);
mpfr_set_z(x0, mant, GMP_RNDN);
x0->_mpfr_sign = 1;
mpfr_set_exp(x0, MIN_EXPONENT);
mpz_clear(mant);
gmp_randclear(state);
/* Chooses incr_prec such that, starting from MIN_PRECISION */
/* and multiplying POINTS_PREC times by incr_prec, we get */
/* something approximately equal to MAX_PRECISION. */
mpfr_init2(temp, 30);
mpfr_init2(temp2, 30);
mpfr_set_ui(temp, MAX_PRECISION, GMP_RNDU);
mpfr_div_ui(temp, temp, MIN_PRECISION, GMP_RNDU);
mpfr_set_ui(temp2, POINTS_PREC, GMP_RNDD);
mpfr_ui_div(temp2, 1, temp2, GMP_RNDU);
mpfr_pow(temp, temp, temp2, GMP_RNDU);
incr_prec = mpfr_get_d(temp, GMP_RNDU);
/* Chooses incr_x such that, starting from 2^MIN_EXPONENT */
/* and multiplying POINTS_X times by incr_x, we get */
/* something approximately equal to 2^MAX_PRECISION */
mpfr_init2(incr_x, 53);
mpfr_set_ui(incr_x, MAX_EXPONENT + 1 - MIN_EXPONENT, GMP_RNDU);
mpfr_div_ui(incr_x, incr_x, POINTS_X, GMP_RNDU);
mpfr_exp2(incr_x, incr_x, GMP_RNDU);
mpfr_clear(temp);
mpfr_clear(temp2);
/* Main loop */
mpfr_init(x);
prec = MIN_PRECISION;
for(i=0; i