diff options
Diffstat (limited to 'tune')
| -rw-r--r-- | tune/tune-gcd-p.c | 132 |
1 files changed, 93 insertions, 39 deletions
diff --git a/tune/tune-gcd-p.c b/tune/tune-gcd-p.c index da4213531..c71949afb 100644 --- a/tune/tune-gcd-p.c +++ b/tune/tune-gcd-p.c @@ -30,6 +30,58 @@ along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ #include "speed.h" +/* Search for minimum over a range. FIXME: Implement golden-section / + fibonacci search*/ +static int +search (double *minp, double (*f)(void *, int), void *ctx, int start, int end) +{ + int x[4]; + double y[4]; + + int best_i; + + x[0] = start; + x[3] = end; + + y[0] = f(ctx, x[0]); + y[3] = f(ctx, x[3]); + + for (;;) + { + int i; + int length = x[3] - x[0]; + + x[1] = x[0] + length/3; + x[2] = x[0] + 2*length/3; + + y[1] = f(ctx, x[1]); + y[2] = f(ctx, x[2]); + +#if 0 + printf("%d: %f, %d: %f, %d:, %f %d: %f\n", + x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3]); +#endif + for (best_i = 0, i = 1; i < 4; i++) + if (y[i] < y[best_i]) + best_i = i; + + if (length <= 4) + break; + + if (best_i >= 2) + { + x[0] = x[1]; + y[0] = y[1]; + } + else + { + x[3] = x[2]; + y[3] = y[2]; + } + } + *minp = y[best_i]; + return x[best_i]; +} static int compare_double(const void *ap, const void *bp) @@ -66,16 +118,38 @@ median (double *v, size_t n) res = median(time_measurement, 5); \ } while (0) -int -main(int argc, char **argv) +struct bench_data { - gmp_randstate_t rands; mp_size_t n; mp_ptr ap; mp_ptr bp; mp_ptr up; mp_ptr vp; mp_ptr gp; +}; + +static double +bench_gcd (void *ctx, int p) +{ + struct bench_data *data = ctx; + double t; + + p_table[data->n] = p; + TIME(t, { + MPN_COPY (data->up, data->ap, data->n); + MPN_COPY (data->vp, data->bp, data->n); + mpn_gcd (data->gp, data->up, data->n, data->vp, data->n); + }); + + return t; +} + +int +main(int argc, char **argv) +{ + gmp_randstate_t rands; struct bench_data data; + mp_size_t n; + TMP_DECL; /* Unbuffered so if output is redirected to a file it isn't lost if the @@ -87,14 +161,14 @@ main(int argc, char **argv) TMP_MARK; - ap = TMP_ALLOC_LIMBS (P_TABLE_SIZE); - bp = TMP_ALLOC_LIMBS (P_TABLE_SIZE); - up = TMP_ALLOC_LIMBS (P_TABLE_SIZE); - vp = TMP_ALLOC_LIMBS (P_TABLE_SIZE); - gp = TMP_ALLOC_LIMBS (P_TABLE_SIZE); + data.ap = TMP_ALLOC_LIMBS (P_TABLE_SIZE); + data.bp = TMP_ALLOC_LIMBS (P_TABLE_SIZE); + data.up = TMP_ALLOC_LIMBS (P_TABLE_SIZE); + data.vp = TMP_ALLOC_LIMBS (P_TABLE_SIZE); + data.gp = TMP_ALLOC_LIMBS (P_TABLE_SIZE); - mpn_random (ap, P_TABLE_SIZE); - mpn_random (bp, P_TABLE_SIZE); + mpn_random (data.ap, P_TABLE_SIZE); + mpn_random (data.bp, P_TABLE_SIZE); memset (p_table, 0, sizeof(p_table)); @@ -105,40 +179,20 @@ main(int argc, char **argv) double best_time; double lehmer_time; - if (ap[n-1] == 0) - ap[n-1] = 1; + if (data.ap[n-1] == 0) + data.ap[n-1] = 1; - if (bp[n-1] == 0) - bp[n-1] = 1; + if (data.bp[n-1] == 0) + data.bp[n-1] = 1; - p_table[n] = 0; - TIME(lehmer_time, { - MPN_COPY (up, ap, n); - MPN_COPY (vp, bp, n); - mpn_gcd (gp, up, n, vp, n); - }); + data.n = n; - best_time = lehmer_time; - best_p = 0; - - for (p = n * 0.48; p < n * 0.77; p++) - { - double t; + lehmer_time = bench_gcd (&data, 0); - p_table[n] = p; + best_p = search (&best_time, bench_gcd, &data, 10, n-10); + if (best_time > lehmer_time) + best_p = 0; - TIME(t, { - MPN_COPY (up, ap, n); - MPN_COPY (vp, bp, n); - mpn_gcd (gp, up, n, vp, n); - }); - - if (t < best_time) - { - best_time = t; - best_p = p; - } - } printf("%6d %6d %5.3g", n, best_p, (double) best_p / n); if (best_p > 0) { |