1 /* tune-gcd-p 2 3 Tune the choice for splitting p in divide-and-conquer gcd. 4 5 Copyright 2008, 2010, 2011 Free Software Foundation, Inc. 6 7 This file is part of the GNU MP Library. 8 9 The GNU MP Library is free software; you can redistribute it and/or modify 10 it under the terms of the GNU Lesser General Public License as published by 11 the Free Software Foundation; either version 3 of the License, or (at your 12 option) any later version. 13 14 The GNU MP Library is distributed in the hope that it will be useful, but 15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17 License for more details. 18 19 You should have received a copy of the GNU Lesser General Public License 20 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 21 22 #define TUNE_GCD_P 1 23 24 #include "../mpn/gcd.c" 25 26 #include <stdio.h> 27 #include <stdlib.h> 28 #include <string.h> 29 #include <time.h> 30 31 #include "speed.h" 32 33 /* Search for minimum over a range. FIXME: Implement golden-section / 34 fibonacci search*/ 35 static int 36 search (double *minp, double (*f)(void *, int), void *ctx, int start, int end) 37 { 38 int x[4]; 39 double y[4]; 40 41 int best_i; 42 43 x[0] = start; 44 x[3] = end; 45 46 y[0] = f(ctx, x[0]); 47 y[3] = f(ctx, x[3]); 48 49 for (;;) 50 { 51 int i; 52 int length = x[3] - x[0]; 53 54 x[1] = x[0] + length/3; 55 x[2] = x[0] + 2*length/3; 56 57 y[1] = f(ctx, x[1]); 58 y[2] = f(ctx, x[2]); 59 60 #if 0 61 printf("%d: %f, %d: %f, %d:, %f %d: %f\n", 62 x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3]); 63 #endif 64 for (best_i = 0, i = 1; i < 4; i++) 65 if (y[i] < y[best_i]) 66 best_i = i; 67 68 if (length <= 4) 69 break; 70 71 if (best_i >= 2) 72 { 73 x[0] = x[1]; 74 y[0] = y[1]; 75 } 76 else 77 { 78 x[3] = x[2]; 79 y[3] = y[2]; 80 } 81 } 82 *minp = y[best_i]; 83 return x[best_i]; 84 } 85 86 static int 87 compare_double(const void *ap, const void *bp) 88 { 89 double a = * (const double *) ap; 90 double b = * (const double *) bp; 91 92 if (a < b) 93 return -1; 94 else if (a > b) 95 return 1; 96 else 97 return 0; 98 } 99 100 static double 101 median (double *v, size_t n) 102 { 103 qsort(v, n, sizeof(*v), compare_double); 104 105 return v[n/2]; 106 } 107 108 #define TIME(res, code) do { \ 109 double time_measurement[5]; \ 110 unsigned time_i; \ 111 \ 112 for (time_i = 0; time_i < 5; time_i++) \ 113 { \ 114 speed_starttime(); \ 115 code; \ 116 time_measurement[time_i] = speed_endtime(); \ 117 } \ 118 res = median(time_measurement, 5); \ 119 } while (0) 120 121 struct bench_data 122 { 123 mp_size_t n; 124 mp_ptr ap; 125 mp_ptr bp; 126 mp_ptr up; 127 mp_ptr vp; 128 mp_ptr gp; 129 }; 130 131 static double 132 bench_gcd (void *ctx, int p) 133 { 134 struct bench_data *data = ctx; 135 double t; 136 137 p_table[data->n] = p; 138 TIME(t, { 139 MPN_COPY (data->up, data->ap, data->n); 140 MPN_COPY (data->vp, data->bp, data->n); 141 mpn_gcd (data->gp, data->up, data->n, data->vp, data->n); 142 }); 143 144 return t; 145 } 146 147 int 148 main(int argc, char **argv) 149 { 150 gmp_randstate_t rands; struct bench_data data; 151 mp_size_t n; 152 153 TMP_DECL; 154 155 /* Unbuffered so if output is redirected to a file it isn't lost if the 156 program is killed part way through. */ 157 setbuf (stdout, NULL); 158 setbuf (stderr, NULL); 159 160 gmp_randinit_default (rands); 161 162 TMP_MARK; 163 164 data.ap = TMP_ALLOC_LIMBS (P_TABLE_SIZE); 165 data.bp = TMP_ALLOC_LIMBS (P_TABLE_SIZE); 166 data.up = TMP_ALLOC_LIMBS (P_TABLE_SIZE); 167 data.vp = TMP_ALLOC_LIMBS (P_TABLE_SIZE); 168 data.gp = TMP_ALLOC_LIMBS (P_TABLE_SIZE); 169 170 mpn_random (data.ap, P_TABLE_SIZE); 171 mpn_random (data.bp, P_TABLE_SIZE); 172 173 memset (p_table, 0, sizeof(p_table)); 174 175 for (n = 100; n < P_TABLE_SIZE; n++) 176 { 177 mp_size_t p; 178 mp_size_t best_p; 179 double best_time; 180 double lehmer_time; 181 182 if (data.ap[n-1] == 0) 183 data.ap[n-1] = 1; 184 185 if (data.bp[n-1] == 0) 186 data.bp[n-1] = 1; 187 188 data.n = n; 189 190 lehmer_time = bench_gcd (&data, 0); 191 192 best_p = search (&best_time, bench_gcd, &data, n/5, 4*n/5); 193 if (best_time > lehmer_time) 194 best_p = 0; 195 196 printf("%6d %6d %5.3g", n, best_p, (double) best_p / n); 197 if (best_p > 0) 198 { 199 double speedup = 100 * (lehmer_time - best_time) / lehmer_time; 200 printf(" %5.3g%%", speedup); 201 if (speedup < 1.0) 202 { 203 printf(" (ignored)"); 204 best_p = 0; 205 } 206 } 207 printf("\n"); 208 209 p_table[n] = best_p; 210 } 211 TMP_FREE; 212 gmp_randclear(rands); 213 return 0; 214 } 215