1 /*
2 * Copyright 2016-2024 The OpenSSL Project Authors. All Rights Reserved.
3 *
4 * Licensed under the Apache License 2.0 (the "License"). You may not use
5 * this file except in compliance with the License. You can obtain a copy
6 * in the file LICENSE in the source distribution or at
7 * https://www.openssl.org/source/license.html
8 */
9
10 /*
11 * This module is meant to be used as template for non-x87 floating-
12 * point assembly modules. The template itself is x86_64-specific
13 * though, as it was debugged on x86_64. So that implementer would
14 * have to recognize platform-specific parts, UxTOy and inline asm,
15 * and act accordingly.
16 *
17 * Huh? x86_64-specific code as template for non-x87? Note seven, which
18 * is not a typo, but reference to 80-bit precision. This module on the
19 * other hand relies on 64-bit precision operations, which are default
20 * for x86_64 code. And since we are at it, just for sense of it,
21 * large-block performance in cycles per processed byte for *this* code
22 * is:
23 * gcc-4.8 icc-15.0 clang-3.4(*)
24 *
25 * Westmere 4.96 5.09 4.37
26 * Sandy Bridge 4.95 4.90 4.17
27 * Haswell 4.92 4.87 3.78
28 * Bulldozer 4.67 4.49 4.68
29 * VIA Nano 7.07 7.05 5.98
30 * Silvermont 10.6 9.61 12.6
31 *
32 * (*) clang managed to discover parallelism and deployed SIMD;
33 *
34 * And for range of other platforms with unspecified gcc versions:
35 *
36 * Freescale e300 12.5
37 * PPC74x0 10.8
38 * POWER6 4.92
39 * POWER7 4.50
40 * POWER8 4.10
41 *
42 * z10 11.2
43 * z196+ 7.30
44 *
45 * UltraSPARC III 16.0
46 * SPARC T4 16.1
47 */
48
49 #if !(defined(__GNUC__) && __GNUC__ >= 2)
50 #error "this is gcc-specific template"
51 #endif
52
53 #include <stdlib.h>
54
55 typedef unsigned char u8;
56 typedef unsigned int u32;
57 typedef unsigned long long u64;
58 typedef union {
59 double d;
60 u64 u;
61 } elem64;
62
63 #define TWO(p) ((double)(1ULL << (p)))
64 #define TWO0 TWO(0)
65 #define TWO32 TWO(32)
66 #define TWO64 (TWO32 * TWO(32))
67 #define TWO96 (TWO64 * TWO(32))
68 #define TWO130 (TWO96 * TWO(34))
69
70 #define EXP(p) ((1023ULL + (p)) << 52)
71
72 #if defined(__x86_64__) || (defined(__PPC__) && defined(__LITTLE_ENDIAN__))
73 #define U8TOU32(p) (*(const u32 *)(p))
74 #define U32TO8(p, v) (*(u32 *)(p) = (v))
75 #elif defined(__PPC__) || defined(__POWERPC__)
76 #define U8TOU32(p) ({u32 ret; asm ("lwbrx %0,0,%1":"=r"(ret):"b"(p)); ret; })
77 #define U32TO8(p, v) asm("stwbrx %0,0,%1" ::"r"(v), "b"(p) : "memory")
78 #elif defined(__s390x__)
79 #define U8TOU32(p) ({u32 ret; asm ("lrv %0,%1":"=d"(ret):"m"(*(u32 *)(p))); ret; })
80 #define U32TO8(p, v) asm("strv %1,%0" : "=m"(*(u32 *)(p)) : "d"(v))
81 #endif
82
83 #ifndef U8TOU32
84 #define U8TOU32(p) ((u32)(p)[0] | (u32)(p)[1] << 8 | (u32)(p)[2] << 16 | (u32)(p)[3] << 24)
85 #endif
86 #ifndef U32TO8
87 #define U32TO8(p, v) ((p)[0] = (u8)(v), (p)[1] = (u8)((v) >> 8), \
88 (p)[2] = (u8)((v) >> 16), (p)[3] = (u8)((v) >> 24))
89 #endif
90
91 typedef struct {
92 elem64 h[4];
93 double r[8];
94 double s[6];
95 } poly1305_internal;
96
97 /* "round toward zero (truncate), mask all exceptions" */
98 #if defined(__x86_64__)
99 static const u32 mxcsr = 0x7f80;
100 #elif defined(__PPC__) || defined(__POWERPC__)
101 static const u64 one = 1;
102 #elif defined(__s390x__)
103 static const u32 fpc = 1;
104 #elif defined(__sparc__)
105 static const u64 fsr = 1ULL << 30;
106 #elif defined(__mips__)
107 static const u32 fcsr = 1;
108 #else
109 #error "unrecognized platform"
110 #endif
111
poly1305_init(void * ctx,const unsigned char key[16])112 int poly1305_init(void *ctx, const unsigned char key[16])
113 {
114 poly1305_internal *st = (poly1305_internal *)ctx;
115 elem64 r0, r1, r2, r3;
116
117 /* h = 0, biased */
118 #if 0
119 st->h[0].d = TWO(52)*TWO0;
120 st->h[1].d = TWO(52)*TWO32;
121 st->h[2].d = TWO(52)*TWO64;
122 st->h[3].d = TWO(52)*TWO96;
123 #else
124 st->h[0].u = EXP(52 + 0);
125 st->h[1].u = EXP(52 + 32);
126 st->h[2].u = EXP(52 + 64);
127 st->h[3].u = EXP(52 + 96);
128 #endif
129
130 if (key) {
131 /*
132 * set "truncate" rounding mode
133 */
134 #if defined(__x86_64__)
135 u32 mxcsr_orig;
136
137 asm volatile("stmxcsr %0" : "=m"(mxcsr_orig));
138 asm volatile("ldmxcsr %0" ::"m"(mxcsr));
139 #elif defined(__PPC__) || defined(__POWERPC__)
140 double fpscr_orig, fpscr = *(double *)&one;
141
142 asm volatile("mffs %0" : "=f"(fpscr_orig));
143 asm volatile("mtfsf 255,%0" ::"f"(fpscr));
144 #elif defined(__s390x__)
145 u32 fpc_orig;
146
147 asm volatile("stfpc %0" : "=m"(fpc_orig));
148 asm volatile("lfpc %0" ::"m"(fpc));
149 #elif defined(__sparc__)
150 u64 fsr_orig;
151
152 asm volatile("stx %%fsr,%0" : "=m"(fsr_orig));
153 asm volatile("ldx %0,%%fsr" ::"m"(fsr));
154 #elif defined(__mips__)
155 u32 fcsr_orig;
156
157 asm volatile("cfc1 %0,$31" : "=r"(fcsr_orig));
158 asm volatile("ctc1 %0,$31" ::"r"(fcsr));
159 #endif
160
161 /* r &= 0xffffffc0ffffffc0ffffffc0fffffff */
162 r0.u = EXP(52 + 0) | (U8TOU32(&key[0]) & 0x0fffffff);
163 r1.u = EXP(52 + 32) | (U8TOU32(&key[4]) & 0x0ffffffc);
164 r2.u = EXP(52 + 64) | (U8TOU32(&key[8]) & 0x0ffffffc);
165 r3.u = EXP(52 + 96) | (U8TOU32(&key[12]) & 0x0ffffffc);
166
167 st->r[0] = r0.d - TWO(52) * TWO0;
168 st->r[2] = r1.d - TWO(52) * TWO32;
169 st->r[4] = r2.d - TWO(52) * TWO64;
170 st->r[6] = r3.d - TWO(52) * TWO96;
171
172 st->s[0] = st->r[2] * (5.0 / TWO130);
173 st->s[2] = st->r[4] * (5.0 / TWO130);
174 st->s[4] = st->r[6] * (5.0 / TWO130);
175
176 /*
177 * base 2^32 -> base 2^16
178 */
179 st->r[1] = (st->r[0] + TWO(52) * TWO(16) * TWO0) - TWO(52) * TWO(16) * TWO0;
180 st->r[0] -= st->r[1];
181
182 st->r[3] = (st->r[2] + TWO(52) * TWO(16) * TWO32) - TWO(52) * TWO(16) * TWO32;
183 st->r[2] -= st->r[3];
184
185 st->r[5] = (st->r[4] + TWO(52) * TWO(16) * TWO64) - TWO(52) * TWO(16) * TWO64;
186 st->r[4] -= st->r[5];
187
188 st->r[7] = (st->r[6] + TWO(52) * TWO(16) * TWO96) - TWO(52) * TWO(16) * TWO96;
189 st->r[6] -= st->r[7];
190
191 st->s[1] = (st->s[0] + TWO(52) * TWO(16) * TWO0 / TWO96) - TWO(52) * TWO(16) * TWO0 / TWO96;
192 st->s[0] -= st->s[1];
193
194 st->s[3] = (st->s[2] + TWO(52) * TWO(16) * TWO32 / TWO96) - TWO(52) * TWO(16) * TWO32 / TWO96;
195 st->s[2] -= st->s[3];
196
197 st->s[5] = (st->s[4] + TWO(52) * TWO(16) * TWO64 / TWO96) - TWO(52) * TWO(16) * TWO64 / TWO96;
198 st->s[4] -= st->s[5];
199
200 /*
201 * restore original FPU control register
202 */
203 #if defined(__x86_64__)
204 asm volatile("ldmxcsr %0" ::"m"(mxcsr_orig));
205 #elif defined(__PPC__) || defined(__POWERPC__)
206 asm volatile("mtfsf 255,%0" ::"f"(fpscr_orig));
207 #elif defined(__s390x__)
208 asm volatile("lfpc %0" ::"m"(fpc_orig));
209 #elif defined(__sparc__)
210 asm volatile("ldx %0,%%fsr" ::"m"(fsr_orig));
211 #elif defined(__mips__)
212 asm volatile("ctc1 %0,$31" ::"r"(fcsr_orig));
213 #endif
214 }
215
216 return 0;
217 }
218
poly1305_blocks(void * ctx,const unsigned char * inp,size_t len,int padbit)219 void poly1305_blocks(void *ctx, const unsigned char *inp, size_t len,
220 int padbit)
221 {
222 poly1305_internal *st = (poly1305_internal *)ctx;
223 elem64 in0, in1, in2, in3;
224 u64 pad = (u64)padbit << 32;
225
226 double x0, x1, x2, x3;
227 double h0lo, h0hi, h1lo, h1hi, h2lo, h2hi, h3lo, h3hi;
228 double c0lo, c0hi, c1lo, c1hi, c2lo, c2hi, c3lo, c3hi;
229
230 const double r0lo = st->r[0];
231 const double r0hi = st->r[1];
232 const double r1lo = st->r[2];
233 const double r1hi = st->r[3];
234 const double r2lo = st->r[4];
235 const double r2hi = st->r[5];
236 const double r3lo = st->r[6];
237 const double r3hi = st->r[7];
238
239 const double s1lo = st->s[0];
240 const double s1hi = st->s[1];
241 const double s2lo = st->s[2];
242 const double s2hi = st->s[3];
243 const double s3lo = st->s[4];
244 const double s3hi = st->s[5];
245
246 /*
247 * set "truncate" rounding mode
248 */
249 #if defined(__x86_64__)
250 u32 mxcsr_orig;
251
252 asm volatile("stmxcsr %0" : "=m"(mxcsr_orig));
253 asm volatile("ldmxcsr %0" ::"m"(mxcsr));
254 #elif defined(__PPC__) || defined(__POWERPC__)
255 double fpscr_orig, fpscr = *(double *)&one;
256
257 asm volatile("mffs %0" : "=f"(fpscr_orig));
258 asm volatile("mtfsf 255,%0" ::"f"(fpscr));
259 #elif defined(__s390x__)
260 u32 fpc_orig;
261
262 asm volatile("stfpc %0" : "=m"(fpc_orig));
263 asm volatile("lfpc %0" ::"m"(fpc));
264 #elif defined(__sparc__)
265 u64 fsr_orig;
266
267 asm volatile("stx %%fsr,%0" : "=m"(fsr_orig));
268 asm volatile("ldx %0,%%fsr" ::"m"(fsr));
269 #elif defined(__mips__)
270 u32 fcsr_orig;
271
272 asm volatile("cfc1 %0,$31" : "=r"(fcsr_orig));
273 asm volatile("ctc1 %0,$31" ::"r"(fcsr));
274 #endif
275
276 /*
277 * load base 2^32 and de-bias
278 */
279 h0lo = st->h[0].d - TWO(52) * TWO0;
280 h1lo = st->h[1].d - TWO(52) * TWO32;
281 h2lo = st->h[2].d - TWO(52) * TWO64;
282 h3lo = st->h[3].d - TWO(52) * TWO96;
283
284 #ifdef __clang__
285 h0hi = 0;
286 h1hi = 0;
287 h2hi = 0;
288 h3hi = 0;
289 #else
290 in0.u = EXP(52 + 0) | U8TOU32(&inp[0]);
291 in1.u = EXP(52 + 32) | U8TOU32(&inp[4]);
292 in2.u = EXP(52 + 64) | U8TOU32(&inp[8]);
293 in3.u = EXP(52 + 96) | U8TOU32(&inp[12]) | pad;
294
295 x0 = in0.d - TWO(52) * TWO0;
296 x1 = in1.d - TWO(52) * TWO32;
297 x2 = in2.d - TWO(52) * TWO64;
298 x3 = in3.d - TWO(52) * TWO96;
299
300 x0 += h0lo;
301 x1 += h1lo;
302 x2 += h2lo;
303 x3 += h3lo;
304
305 goto fast_entry;
306 #endif
307
308 do {
309 in0.u = EXP(52 + 0) | U8TOU32(&inp[0]);
310 in1.u = EXP(52 + 32) | U8TOU32(&inp[4]);
311 in2.u = EXP(52 + 64) | U8TOU32(&inp[8]);
312 in3.u = EXP(52 + 96) | U8TOU32(&inp[12]) | pad;
313
314 x0 = in0.d - TWO(52) * TWO0;
315 x1 = in1.d - TWO(52) * TWO32;
316 x2 = in2.d - TWO(52) * TWO64;
317 x3 = in3.d - TWO(52) * TWO96;
318
319 /*
320 * note that there are multiple ways to accumulate input, e.g.
321 * one can as well accumulate to h0lo-h1lo-h1hi-h2hi...
322 */
323 h0lo += x0;
324 h0hi += x1;
325 h2lo += x2;
326 h2hi += x3;
327
328 /*
329 * carries that cross 32n-bit (and 130-bit) boundaries
330 */
331 c0lo = (h0lo + TWO(52) * TWO32) - TWO(52) * TWO32;
332 c1lo = (h1lo + TWO(52) * TWO64) - TWO(52) * TWO64;
333 c2lo = (h2lo + TWO(52) * TWO96) - TWO(52) * TWO96;
334 c3lo = (h3lo + TWO(52) * TWO130) - TWO(52) * TWO130;
335
336 c0hi = (h0hi + TWO(52) * TWO32) - TWO(52) * TWO32;
337 c1hi = (h1hi + TWO(52) * TWO64) - TWO(52) * TWO64;
338 c2hi = (h2hi + TWO(52) * TWO96) - TWO(52) * TWO96;
339 c3hi = (h3hi + TWO(52) * TWO130) - TWO(52) * TWO130;
340
341 /*
342 * base 2^48 -> base 2^32 with last reduction step
343 */
344 x1 = (h1lo - c1lo) + c0lo;
345 x2 = (h2lo - c2lo) + c1lo;
346 x3 = (h3lo - c3lo) + c2lo;
347 x0 = (h0lo - c0lo) + c3lo * (5.0 / TWO130);
348
349 x1 += (h1hi - c1hi) + c0hi;
350 x2 += (h2hi - c2hi) + c1hi;
351 x3 += (h3hi - c3hi) + c2hi;
352 x0 += (h0hi - c0hi) + c3hi * (5.0 / TWO130);
353
354 #ifndef __clang__
355 fast_entry:
356 #endif
357 /*
358 * base 2^32 * base 2^16 = base 2^48
359 */
360 h0lo = s3lo * x1 + s2lo * x2 + s1lo * x3 + r0lo * x0;
361 h1lo = r0lo * x1 + s3lo * x2 + s2lo * x3 + r1lo * x0;
362 h2lo = r1lo * x1 + r0lo * x2 + s3lo * x3 + r2lo * x0;
363 h3lo = r2lo * x1 + r1lo * x2 + r0lo * x3 + r3lo * x0;
364
365 h0hi = s3hi * x1 + s2hi * x2 + s1hi * x3 + r0hi * x0;
366 h1hi = r0hi * x1 + s3hi * x2 + s2hi * x3 + r1hi * x0;
367 h2hi = r1hi * x1 + r0hi * x2 + s3hi * x3 + r2hi * x0;
368 h3hi = r2hi * x1 + r1hi * x2 + r0hi * x3 + r3hi * x0;
369
370 inp += 16;
371 len -= 16;
372
373 } while (len >= 16);
374
375 /*
376 * carries that cross 32n-bit (and 130-bit) boundaries
377 */
378 c0lo = (h0lo + TWO(52) * TWO32) - TWO(52) * TWO32;
379 c1lo = (h1lo + TWO(52) * TWO64) - TWO(52) * TWO64;
380 c2lo = (h2lo + TWO(52) * TWO96) - TWO(52) * TWO96;
381 c3lo = (h3lo + TWO(52) * TWO130) - TWO(52) * TWO130;
382
383 c0hi = (h0hi + TWO(52) * TWO32) - TWO(52) * TWO32;
384 c1hi = (h1hi + TWO(52) * TWO64) - TWO(52) * TWO64;
385 c2hi = (h2hi + TWO(52) * TWO96) - TWO(52) * TWO96;
386 c3hi = (h3hi + TWO(52) * TWO130) - TWO(52) * TWO130;
387
388 /*
389 * base 2^48 -> base 2^32 with last reduction step
390 */
391 x1 = (h1lo - c1lo) + c0lo;
392 x2 = (h2lo - c2lo) + c1lo;
393 x3 = (h3lo - c3lo) + c2lo;
394 x0 = (h0lo - c0lo) + c3lo * (5.0 / TWO130);
395
396 x1 += (h1hi - c1hi) + c0hi;
397 x2 += (h2hi - c2hi) + c1hi;
398 x3 += (h3hi - c3hi) + c2hi;
399 x0 += (h0hi - c0hi) + c3hi * (5.0 / TWO130);
400
401 /*
402 * store base 2^32, with bias
403 */
404 st->h[1].d = x1 + TWO(52) * TWO32;
405 st->h[2].d = x2 + TWO(52) * TWO64;
406 st->h[3].d = x3 + TWO(52) * TWO96;
407 st->h[0].d = x0 + TWO(52) * TWO0;
408
409 /*
410 * restore original FPU control register
411 */
412 #if defined(__x86_64__)
413 asm volatile("ldmxcsr %0" ::"m"(mxcsr_orig));
414 #elif defined(__PPC__) || defined(__POWERPC__)
415 asm volatile("mtfsf 255,%0" ::"f"(fpscr_orig));
416 #elif defined(__s390x__)
417 asm volatile("lfpc %0" ::"m"(fpc_orig));
418 #elif defined(__sparc__)
419 asm volatile("ldx %0,%%fsr" ::"m"(fsr_orig));
420 #elif defined(__mips__)
421 asm volatile("ctc1 %0,$31" ::"r"(fcsr_orig));
422 #endif
423 }
424
poly1305_emit(void * ctx,unsigned char mac[16],const u32 nonce[4])425 void poly1305_emit(void *ctx, unsigned char mac[16], const u32 nonce[4])
426 {
427 poly1305_internal *st = (poly1305_internal *)ctx;
428 u64 h0, h1, h2, h3, h4;
429 u32 g0, g1, g2, g3, g4;
430 u64 t;
431 u32 mask;
432
433 /*
434 * thanks to bias masking exponent gives integer result
435 */
436 h0 = st->h[0].u & 0x000fffffffffffffULL;
437 h1 = st->h[1].u & 0x000fffffffffffffULL;
438 h2 = st->h[2].u & 0x000fffffffffffffULL;
439 h3 = st->h[3].u & 0x000fffffffffffffULL;
440
441 /*
442 * can be partially reduced, so reduce...
443 */
444 h4 = h3 >> 32;
445 h3 &= 0xffffffffU;
446 g4 = h4 & -4;
447 h4 &= 3;
448 g4 += g4 >> 2;
449
450 h0 += g4;
451 h1 += h0 >> 32;
452 h0 &= 0xffffffffU;
453 h2 += h1 >> 32;
454 h1 &= 0xffffffffU;
455 h3 += h2 >> 32;
456 h2 &= 0xffffffffU;
457
458 /* compute h + -p */
459 g0 = (u32)(t = h0 + 5);
460 g1 = (u32)(t = h1 + (t >> 32));
461 g2 = (u32)(t = h2 + (t >> 32));
462 g3 = (u32)(t = h3 + (t >> 32));
463 g4 = h4 + (u32)(t >> 32);
464
465 /* if there was carry, select g0-g3 */
466 mask = 0 - (g4 >> 2);
467 g0 &= mask;
468 g1 &= mask;
469 g2 &= mask;
470 g3 &= mask;
471 mask = ~mask;
472 g0 |= (h0 & mask);
473 g1 |= (h1 & mask);
474 g2 |= (h2 & mask);
475 g3 |= (h3 & mask);
476
477 /* mac = (h + nonce) % (2^128) */
478 g0 = (u32)(t = (u64)g0 + nonce[0]);
479 g1 = (u32)(t = (u64)g1 + (t >> 32) + nonce[1]);
480 g2 = (u32)(t = (u64)g2 + (t >> 32) + nonce[2]);
481 g3 = (u32)(t = (u64)g3 + (t >> 32) + nonce[3]);
482
483 U32TO8(mac + 0, g0);
484 U32TO8(mac + 4, g1);
485 U32TO8(mac + 8, g2);
486 U32TO8(mac + 12, g3);
487 }
488