x86_64-gcc.c revision 4969cc9b0ab2905ec478277f50ed3849b37a6c6b
1/* x86_64 BIGNUM accelerator version 0.1, December 2002.
2 *
3 * Implemented by Andy Polyakov <appro@fy.chalmers.se> for the OpenSSL
4 * project.
5 *
6 * Rights for redistribution and usage in source and binary forms are
7 * granted according to the OpenSSL license. Warranty of any kind is
8 * disclaimed.
9 *
10 * Q. Version 0.1? It doesn't sound like Andy, he used to assign real
11 *    versions, like 1.0...
12 * A. Well, that's because this code is basically a quick-n-dirty
13 *    proof-of-concept hack. As you can see it's implemented with
14 *    inline assembler, which means that you're bound to GCC and that
15 *    there might be enough room for further improvement.
16 *
17 * Q. Why inline assembler?
18 * A. x86_64 features own ABI which I'm not familiar with. This is
19 *    why I decided to let the compiler take care of subroutine
20 *    prologue/epilogue as well as register allocation. For reference.
21 *    Win64 implements different ABI for AMD64, different from Linux.
22 *
23 * Q. How much faster does it get?
24 * A. 'apps/openssl speed rsa dsa' output with no-asm:
25 *
26 *	                  sign    verify    sign/s verify/s
27 *	rsa  512 bits   0.0006s   0.0001s   1683.8  18456.2
28 *	rsa 1024 bits   0.0028s   0.0002s    356.0   6407.0
29 *	rsa 2048 bits   0.0172s   0.0005s     58.0   1957.8
30 *	rsa 4096 bits   0.1155s   0.0018s      8.7    555.6
31 *	                  sign    verify    sign/s verify/s
32 *	dsa  512 bits   0.0005s   0.0006s   2100.8   1768.3
33 *	dsa 1024 bits   0.0014s   0.0018s    692.3    559.2
34 *	dsa 2048 bits   0.0049s   0.0061s    204.7    165.0
35 *
36 *    'apps/openssl speed rsa dsa' output with this module:
37 *
38 *	                  sign    verify    sign/s verify/s
39 *	rsa  512 bits   0.0004s   0.0000s   2767.1  33297.9
40 *	rsa 1024 bits   0.0012s   0.0001s    867.4  14674.7
41 *	rsa 2048 bits   0.0061s   0.0002s    164.0   5270.0
42 *	rsa 4096 bits   0.0384s   0.0006s     26.1   1650.8
43 *	                  sign    verify    sign/s verify/s
44 *	dsa  512 bits   0.0002s   0.0003s   4442.2   3786.3
45 *	dsa 1024 bits   0.0005s   0.0007s   1835.1   1497.4
46 *	dsa 2048 bits   0.0016s   0.0020s    620.4    504.6
47 *
48 *    For the reference. IA-32 assembler implementation performs
49 *    very much like 64-bit code compiled with no-asm on the same
50 *    machine.
51 */
52
53#include <openssl/bn.h>
54
55/* TODO(davidben): Get this file working on Windows x64. */
56#if !defined(OPENSSL_NO_ASM) && defined(OPENSSL_X86_64) && defined(__GNUC__)
57
58#include "../internal.h"
59
60
61#undef mul
62#undef mul_add
63
64#define asm __asm__
65
66/*
67 * "m"(a), "+m"(r)	is the way to favor DirectPath µ-code;
68 * "g"(0)		let the compiler to decide where does it
69 *			want to keep the value of zero;
70 */
71#define mul_add(r, a, word, carry)                                     \
72  do {                                                                 \
73    register BN_ULONG high, low;                                       \
74    asm("mulq %3" : "=a"(low), "=d"(high) : "a"(word), "m"(a) : "cc"); \
75    asm("addq %2,%0; adcq %3,%1"                                       \
76        : "+r"(carry), "+d"(high)                                      \
77        : "a"(low), "g"(0)                                             \
78        : "cc");                                                       \
79    asm("addq %2,%0; adcq %3,%1"                                       \
80        : "+m"(r), "+d"(high)                                          \
81        : "r"(carry), "g"(0)                                           \
82        : "cc");                                                       \
83    carry = high;                                                      \
84  } while (0)
85
86#define mul(r, a, word, carry)                                         \
87  do {                                                                 \
88    register BN_ULONG high, low;                                       \
89    asm("mulq %3" : "=a"(low), "=d"(high) : "a"(word), "g"(a) : "cc"); \
90    asm("addq %2,%0; adcq %3,%1"                                       \
91        : "+r"(carry), "+d"(high)                                      \
92        : "a"(low), "g"(0)                                             \
93        : "cc");                                                       \
94    (r) = carry, carry = high;                                         \
95  } while (0)
96#undef sqr
97#define sqr(r0, r1, a) asm("mulq %2" : "=a"(r0), "=d"(r1) : "a"(a) : "cc");
98
99BN_ULONG bn_mul_add_words(BN_ULONG *rp, const BN_ULONG *ap, int num,
100                          BN_ULONG w) {
101  BN_ULONG c1 = 0;
102
103  if (num <= 0) {
104    return (c1);
105  }
106
107  while (num & ~3) {
108    mul_add(rp[0], ap[0], w, c1);
109    mul_add(rp[1], ap[1], w, c1);
110    mul_add(rp[2], ap[2], w, c1);
111    mul_add(rp[3], ap[3], w, c1);
112    ap += 4;
113    rp += 4;
114    num -= 4;
115  }
116  if (num) {
117    mul_add(rp[0], ap[0], w, c1);
118    if (--num == 0) {
119      return c1;
120    }
121    mul_add(rp[1], ap[1], w, c1);
122    if (--num == 0) {
123      return c1;
124    }
125    mul_add(rp[2], ap[2], w, c1);
126    return c1;
127  }
128
129  return c1;
130}
131
132BN_ULONG bn_mul_words(BN_ULONG *rp, const BN_ULONG *ap, int num, BN_ULONG w) {
133  BN_ULONG c1 = 0;
134
135  if (num <= 0) {
136    return c1;
137  }
138
139  while (num & ~3) {
140    mul(rp[0], ap[0], w, c1);
141    mul(rp[1], ap[1], w, c1);
142    mul(rp[2], ap[2], w, c1);
143    mul(rp[3], ap[3], w, c1);
144    ap += 4;
145    rp += 4;
146    num -= 4;
147  }
148  if (num) {
149    mul(rp[0], ap[0], w, c1);
150    if (--num == 0) {
151      return c1;
152    }
153    mul(rp[1], ap[1], w, c1);
154    if (--num == 0) {
155      return c1;
156    }
157    mul(rp[2], ap[2], w, c1);
158  }
159  return c1;
160}
161
162void bn_sqr_words(BN_ULONG *r, const BN_ULONG *a, int n) {
163  if (n <= 0) {
164    return;
165  }
166
167  while (n & ~3) {
168    sqr(r[0], r[1], a[0]);
169    sqr(r[2], r[3], a[1]);
170    sqr(r[4], r[5], a[2]);
171    sqr(r[6], r[7], a[3]);
172    a += 4;
173    r += 8;
174    n -= 4;
175  }
176  if (n) {
177    sqr(r[0], r[1], a[0]);
178    if (--n == 0) {
179      return;
180    }
181    sqr(r[2], r[3], a[1]);
182    if (--n == 0) {
183      return;
184    }
185    sqr(r[4], r[5], a[2]);
186  }
187}
188
189BN_ULONG bn_add_words(BN_ULONG *rp, const BN_ULONG *ap, const BN_ULONG *bp,
190                      int n) {
191  BN_ULONG ret;
192  size_t i = 0;
193
194  if (n <= 0) {
195    return 0;
196  }
197
198  asm volatile (
199      "	subq	%0,%0		\n" /* clear carry */
200      "	jmp	1f		\n"
201      ".p2align 4			\n"
202      "1:	movq	(%4,%2,8),%0	\n"
203      "	adcq	(%5,%2,8),%0	\n"
204      "	movq	%0,(%3,%2,8)	\n"
205      "	lea	1(%2),%2	\n"
206      "	loop	1b		\n"
207      "	sbbq	%0,%0		\n"
208      : "=&r"(ret), "+c"(n), "+r"(i)
209      : "r"(rp), "r"(ap), "r"(bp)
210      : "cc", "memory");
211
212  return ret & 1;
213}
214
215BN_ULONG bn_sub_words(BN_ULONG *rp, const BN_ULONG *ap, const BN_ULONG *bp,
216                      int n) {
217  BN_ULONG ret;
218  size_t i = 0;
219
220  if (n <= 0) {
221    return 0;
222  }
223
224  asm volatile (
225      "	subq	%0,%0		\n" /* clear borrow */
226      "	jmp	1f		\n"
227      ".p2align 4			\n"
228      "1:	movq	(%4,%2,8),%0	\n"
229      "	sbbq	(%5,%2,8),%0	\n"
230      "	movq	%0,(%3,%2,8)	\n"
231      "	lea	1(%2),%2	\n"
232      "	loop	1b		\n"
233      "	sbbq	%0,%0		\n"
234      : "=&r"(ret), "+c"(n), "+r"(i)
235      : "r"(rp), "r"(ap), "r"(bp)
236      : "cc", "memory");
237
238  return ret & 1;
239}
240
241/* mul_add_c(a,b,c0,c1,c2)  -- c+=a*b for three word number c=(c2,c1,c0) */
242/* mul_add_c2(a,b,c0,c1,c2) -- c+=2*a*b for three word number c=(c2,c1,c0) */
243/* sqr_add_c(a,i,c0,c1,c2)  -- c+=a[i]^2 for three word number c=(c2,c1,c0) */
244/* sqr_add_c2(a,i,c0,c1,c2) -- c+=2*a[i]*a[j] for three word number c=(c2,c1,c0)
245 */
246
247/* Keep in mind that carrying into high part of multiplication result can not
248 * overflow, because it cannot be all-ones. */
249#define mul_add_c(a, b, c0, c1, c2)          \
250  do {                                       \
251    BN_ULONG t1, t2;                \
252    asm("mulq %3" : "=a"(t1), "=d"(t2) : "a"(a), "m"(b) : "cc"); \
253    asm("addq %3,%0; adcq %4,%1; adcq %5,%2" \
254        : "+r"(c0), "+r"(c1), "+r"(c2)       \
255        : "r"(t1), "r"(t2), "g"(0)           \
256        : "cc");                             \
257  } while (0)
258
259#define sqr_add_c(a, i, c0, c1, c2)          \
260  do {                                       \
261    BN_ULONG t1, t2;                         \
262    asm("mulq %2" : "=a"(t1), "=d"(t2) : "a"(a[i]) : "cc"); \
263    asm("addq %3,%0; adcq %4,%1; adcq %5,%2" \
264        : "+r"(c0), "+r"(c1), "+r"(c2)       \
265        : "r"(t1), "r"(t2), "g"(0)           \
266        : "cc");                             \
267  } while (0)
268
269#define mul_add_c2(a, b, c0, c1, c2)         \
270  do {                                       \
271    BN_ULONG t1, t2;                                                    \
272    asm("mulq %3" : "=a"(t1), "=d"(t2) : "a"(a), "m"(b) : "cc");        \
273    asm("addq %3,%0; adcq %4,%1; adcq %5,%2" \
274        : "+r"(c0), "+r"(c1), "+r"(c2)       \
275        : "r"(t1), "r"(t2), "g"(0)           \
276        : "cc");                             \
277    asm("addq %3,%0; adcq %4,%1; adcq %5,%2" \
278        : "+r"(c0), "+r"(c1), "+r"(c2)       \
279        : "r"(t1), "r"(t2), "g"(0)           \
280        : "cc");                             \
281  } while (0)
282
283#define sqr_add_c2(a, i, j, c0, c1, c2) mul_add_c2((a)[i], (a)[j], c0, c1, c2)
284
285void bn_mul_comba8(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b) {
286  BN_ULONG c1, c2, c3;
287
288  c1 = 0;
289  c2 = 0;
290  c3 = 0;
291  mul_add_c(a[0], b[0], c1, c2, c3);
292  r[0] = c1;
293  c1 = 0;
294  mul_add_c(a[0], b[1], c2, c3, c1);
295  mul_add_c(a[1], b[0], c2, c3, c1);
296  r[1] = c2;
297  c2 = 0;
298  mul_add_c(a[2], b[0], c3, c1, c2);
299  mul_add_c(a[1], b[1], c3, c1, c2);
300  mul_add_c(a[0], b[2], c3, c1, c2);
301  r[2] = c3;
302  c3 = 0;
303  mul_add_c(a[0], b[3], c1, c2, c3);
304  mul_add_c(a[1], b[2], c1, c2, c3);
305  mul_add_c(a[2], b[1], c1, c2, c3);
306  mul_add_c(a[3], b[0], c1, c2, c3);
307  r[3] = c1;
308  c1 = 0;
309  mul_add_c(a[4], b[0], c2, c3, c1);
310  mul_add_c(a[3], b[1], c2, c3, c1);
311  mul_add_c(a[2], b[2], c2, c3, c1);
312  mul_add_c(a[1], b[3], c2, c3, c1);
313  mul_add_c(a[0], b[4], c2, c3, c1);
314  r[4] = c2;
315  c2 = 0;
316  mul_add_c(a[0], b[5], c3, c1, c2);
317  mul_add_c(a[1], b[4], c3, c1, c2);
318  mul_add_c(a[2], b[3], c3, c1, c2);
319  mul_add_c(a[3], b[2], c3, c1, c2);
320  mul_add_c(a[4], b[1], c3, c1, c2);
321  mul_add_c(a[5], b[0], c3, c1, c2);
322  r[5] = c3;
323  c3 = 0;
324  mul_add_c(a[6], b[0], c1, c2, c3);
325  mul_add_c(a[5], b[1], c1, c2, c3);
326  mul_add_c(a[4], b[2], c1, c2, c3);
327  mul_add_c(a[3], b[3], c1, c2, c3);
328  mul_add_c(a[2], b[4], c1, c2, c3);
329  mul_add_c(a[1], b[5], c1, c2, c3);
330  mul_add_c(a[0], b[6], c1, c2, c3);
331  r[6] = c1;
332  c1 = 0;
333  mul_add_c(a[0], b[7], c2, c3, c1);
334  mul_add_c(a[1], b[6], c2, c3, c1);
335  mul_add_c(a[2], b[5], c2, c3, c1);
336  mul_add_c(a[3], b[4], c2, c3, c1);
337  mul_add_c(a[4], b[3], c2, c3, c1);
338  mul_add_c(a[5], b[2], c2, c3, c1);
339  mul_add_c(a[6], b[1], c2, c3, c1);
340  mul_add_c(a[7], b[0], c2, c3, c1);
341  r[7] = c2;
342  c2 = 0;
343  mul_add_c(a[7], b[1], c3, c1, c2);
344  mul_add_c(a[6], b[2], c3, c1, c2);
345  mul_add_c(a[5], b[3], c3, c1, c2);
346  mul_add_c(a[4], b[4], c3, c1, c2);
347  mul_add_c(a[3], b[5], c3, c1, c2);
348  mul_add_c(a[2], b[6], c3, c1, c2);
349  mul_add_c(a[1], b[7], c3, c1, c2);
350  r[8] = c3;
351  c3 = 0;
352  mul_add_c(a[2], b[7], c1, c2, c3);
353  mul_add_c(a[3], b[6], c1, c2, c3);
354  mul_add_c(a[4], b[5], c1, c2, c3);
355  mul_add_c(a[5], b[4], c1, c2, c3);
356  mul_add_c(a[6], b[3], c1, c2, c3);
357  mul_add_c(a[7], b[2], c1, c2, c3);
358  r[9] = c1;
359  c1 = 0;
360  mul_add_c(a[7], b[3], c2, c3, c1);
361  mul_add_c(a[6], b[4], c2, c3, c1);
362  mul_add_c(a[5], b[5], c2, c3, c1);
363  mul_add_c(a[4], b[6], c2, c3, c1);
364  mul_add_c(a[3], b[7], c2, c3, c1);
365  r[10] = c2;
366  c2 = 0;
367  mul_add_c(a[4], b[7], c3, c1, c2);
368  mul_add_c(a[5], b[6], c3, c1, c2);
369  mul_add_c(a[6], b[5], c3, c1, c2);
370  mul_add_c(a[7], b[4], c3, c1, c2);
371  r[11] = c3;
372  c3 = 0;
373  mul_add_c(a[7], b[5], c1, c2, c3);
374  mul_add_c(a[6], b[6], c1, c2, c3);
375  mul_add_c(a[5], b[7], c1, c2, c3);
376  r[12] = c1;
377  c1 = 0;
378  mul_add_c(a[6], b[7], c2, c3, c1);
379  mul_add_c(a[7], b[6], c2, c3, c1);
380  r[13] = c2;
381  c2 = 0;
382  mul_add_c(a[7], b[7], c3, c1, c2);
383  r[14] = c3;
384  r[15] = c1;
385}
386
387void bn_mul_comba4(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b) {
388  BN_ULONG c1, c2, c3;
389
390  c1 = 0;
391  c2 = 0;
392  c3 = 0;
393  mul_add_c(a[0], b[0], c1, c2, c3);
394  r[0] = c1;
395  c1 = 0;
396  mul_add_c(a[0], b[1], c2, c3, c1);
397  mul_add_c(a[1], b[0], c2, c3, c1);
398  r[1] = c2;
399  c2 = 0;
400  mul_add_c(a[2], b[0], c3, c1, c2);
401  mul_add_c(a[1], b[1], c3, c1, c2);
402  mul_add_c(a[0], b[2], c3, c1, c2);
403  r[2] = c3;
404  c3 = 0;
405  mul_add_c(a[0], b[3], c1, c2, c3);
406  mul_add_c(a[1], b[2], c1, c2, c3);
407  mul_add_c(a[2], b[1], c1, c2, c3);
408  mul_add_c(a[3], b[0], c1, c2, c3);
409  r[3] = c1;
410  c1 = 0;
411  mul_add_c(a[3], b[1], c2, c3, c1);
412  mul_add_c(a[2], b[2], c2, c3, c1);
413  mul_add_c(a[1], b[3], c2, c3, c1);
414  r[4] = c2;
415  c2 = 0;
416  mul_add_c(a[2], b[3], c3, c1, c2);
417  mul_add_c(a[3], b[2], c3, c1, c2);
418  r[5] = c3;
419  c3 = 0;
420  mul_add_c(a[3], b[3], c1, c2, c3);
421  r[6] = c1;
422  r[7] = c2;
423}
424
425void bn_sqr_comba8(BN_ULONG *r, const BN_ULONG *a) {
426  BN_ULONG c1, c2, c3;
427
428  c1 = 0;
429  c2 = 0;
430  c3 = 0;
431  sqr_add_c(a, 0, c1, c2, c3);
432  r[0] = c1;
433  c1 = 0;
434  sqr_add_c2(a, 1, 0, c2, c3, c1);
435  r[1] = c2;
436  c2 = 0;
437  sqr_add_c(a, 1, c3, c1, c2);
438  sqr_add_c2(a, 2, 0, c3, c1, c2);
439  r[2] = c3;
440  c3 = 0;
441  sqr_add_c2(a, 3, 0, c1, c2, c3);
442  sqr_add_c2(a, 2, 1, c1, c2, c3);
443  r[3] = c1;
444  c1 = 0;
445  sqr_add_c(a, 2, c2, c3, c1);
446  sqr_add_c2(a, 3, 1, c2, c3, c1);
447  sqr_add_c2(a, 4, 0, c2, c3, c1);
448  r[4] = c2;
449  c2 = 0;
450  sqr_add_c2(a, 5, 0, c3, c1, c2);
451  sqr_add_c2(a, 4, 1, c3, c1, c2);
452  sqr_add_c2(a, 3, 2, c3, c1, c2);
453  r[5] = c3;
454  c3 = 0;
455  sqr_add_c(a, 3, c1, c2, c3);
456  sqr_add_c2(a, 4, 2, c1, c2, c3);
457  sqr_add_c2(a, 5, 1, c1, c2, c3);
458  sqr_add_c2(a, 6, 0, c1, c2, c3);
459  r[6] = c1;
460  c1 = 0;
461  sqr_add_c2(a, 7, 0, c2, c3, c1);
462  sqr_add_c2(a, 6, 1, c2, c3, c1);
463  sqr_add_c2(a, 5, 2, c2, c3, c1);
464  sqr_add_c2(a, 4, 3, c2, c3, c1);
465  r[7] = c2;
466  c2 = 0;
467  sqr_add_c(a, 4, c3, c1, c2);
468  sqr_add_c2(a, 5, 3, c3, c1, c2);
469  sqr_add_c2(a, 6, 2, c3, c1, c2);
470  sqr_add_c2(a, 7, 1, c3, c1, c2);
471  r[8] = c3;
472  c3 = 0;
473  sqr_add_c2(a, 7, 2, c1, c2, c3);
474  sqr_add_c2(a, 6, 3, c1, c2, c3);
475  sqr_add_c2(a, 5, 4, c1, c2, c3);
476  r[9] = c1;
477  c1 = 0;
478  sqr_add_c(a, 5, c2, c3, c1);
479  sqr_add_c2(a, 6, 4, c2, c3, c1);
480  sqr_add_c2(a, 7, 3, c2, c3, c1);
481  r[10] = c2;
482  c2 = 0;
483  sqr_add_c2(a, 7, 4, c3, c1, c2);
484  sqr_add_c2(a, 6, 5, c3, c1, c2);
485  r[11] = c3;
486  c3 = 0;
487  sqr_add_c(a, 6, c1, c2, c3);
488  sqr_add_c2(a, 7, 5, c1, c2, c3);
489  r[12] = c1;
490  c1 = 0;
491  sqr_add_c2(a, 7, 6, c2, c3, c1);
492  r[13] = c2;
493  c2 = 0;
494  sqr_add_c(a, 7, c3, c1, c2);
495  r[14] = c3;
496  r[15] = c1;
497}
498
499void bn_sqr_comba4(BN_ULONG *r, const BN_ULONG *a) {
500  BN_ULONG c1, c2, c3;
501
502  c1 = 0;
503  c2 = 0;
504  c3 = 0;
505  sqr_add_c(a, 0, c1, c2, c3);
506  r[0] = c1;
507  c1 = 0;
508  sqr_add_c2(a, 1, 0, c2, c3, c1);
509  r[1] = c2;
510  c2 = 0;
511  sqr_add_c(a, 1, c3, c1, c2);
512  sqr_add_c2(a, 2, 0, c3, c1, c2);
513  r[2] = c3;
514  c3 = 0;
515  sqr_add_c2(a, 3, 0, c1, c2, c3);
516  sqr_add_c2(a, 2, 1, c1, c2, c3);
517  r[3] = c1;
518  c1 = 0;
519  sqr_add_c(a, 2, c2, c3, c1);
520  sqr_add_c2(a, 3, 1, c2, c3, c1);
521  r[4] = c2;
522  c2 = 0;
523  sqr_add_c2(a, 3, 2, c3, c1, c2);
524  r[5] = c3;
525  c3 = 0;
526  sqr_add_c(a, 3, c1, c2, c3);
527  r[6] = c1;
528  r[7] = c2;
529}
530
531#endif  /* !NO_ASM && X86_64 && __GNUC__ */
532