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