5 * (c) 1999 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
12 * Catacomb is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU Library General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
17 * Catacomb is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Library General Public License for more details.
22 * You should have received a copy of the GNU Library General Public
23 * License along with Catacomb; if not, write to the Free
24 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
28 /*----- Header files ------------------------------------------------------*/
35 /*----- Tweakables --------------------------------------------------------*/
37 /* --- @MPMONT_DISABLE@ --- *
39 * Replace all the clever Montgomery reduction with good old-fashioned long
43 /* #define MPMONT_DISABLE */
45 #define MPMONT_KTHRESH (16*MPK_THRESH)
47 /*----- Low-level implementation ------------------------------------------*/
49 #ifndef MPMONT_DISABLE
51 /* --- @redccore@ --- *
53 * Arguments: @mpw *dv, *dvl@ = base and limit of source/destination
54 * @const mpw *mv@ = base of modulus %$m$%
55 * @size_t n@ = length of modulus
56 * @const mpw *mi@ = base of REDC coefficient %$m'$%
60 * Use: Let %$a$% be the input operand. Store in %$d$% the value
61 * %$a + (m' a \bmod R) m$%. The destination has space for at
62 * least %$2 n + 1$% words of result.
65 CPU_DISPATCH(static, (void), void, redccore,
66 (mpw *dv, mpw *dvl, const mpw *mv, size_t n, const mpw *mi),
67 (dv, dvl, mv, n, mi), pick_redccore, simple_redccore);
69 static void simple_redccore(mpw *dv, mpw *dvl, const mpw *mv,
70 size_t n, const mpw *mi)
75 for (i = 0; i < n; i++) {
76 MPX_UMLAN(dv, dvl, mv, mv + n, MPW(*dv*mi0));
81 #define MAYBE_REDC4(impl) \
82 extern void mpxmont_redc4_##impl(mpw *dv, mpw *dvl, const mpw *mv, \
83 size_t n, const mpw *mi); \
84 static void maybe_redc4_##impl(mpw *dv, mpw *dvl, const mpw *mv, \
85 size_t n, const mpw *mi) \
87 if (n%4) simple_redccore(dv, dvl, mv, n, mi); \
88 else mpxmont_redc4_##impl(dv, dvl, mv, n, mi); \
97 MAYBE_REDC4(amd64_sse2)
98 MAYBE_REDC4(amd64_avx)
102 MAYBE_REDC4(arm_neon)
106 MAYBE_REDC4(arm64_simd)
109 static redccore__functype *pick_redccore(void)
112 DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_x86_avx,
113 cpu_feature_p(CPUFEAT_X86_AVX));
114 DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_x86_sse2,
115 cpu_feature_p(CPUFEAT_X86_SSE2));
118 DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_amd64_avx,
119 cpu_feature_p(CPUFEAT_X86_AVX));
120 DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_amd64_sse2,
121 cpu_feature_p(CPUFEAT_X86_SSE2));
124 DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_arm_neon,
125 cpu_feature_p(CPUFEAT_ARM_NEON));
128 DISPATCH_PICK_COND(mpmont_reduce, maybe_redc4_arm64_simd,
129 cpu_feature_p(CPUFEAT_ARM_NEON));
131 DISPATCH_PICK_FALLBACK(mpmont_reduce, simple_redccore);
134 /* --- @mulcore@ --- *
136 * Arguments: @mpw *dv, *dvl@ = base and limit of source/destination
137 * @const mpw *av, *avl@ = base and limit of first multiplicand
138 * @const mpw *bv, *bvl@ = base and limit of second multiplicand
139 * @const mpw *mv@ = base of modulus %$m$%
140 * @size_t n@ = length of modulus
141 * @const mpw *mi@ = base of REDC coefficient %$m'$%
145 * Use: Let %$a$% and %$b$% be the multiplicands. Let %$w = a b$%.
146 * Store in %$d$% the value %$a b + (m' a b \bmod R) m$%.
149 CPU_DISPATCH(static, (void), void, mulcore,
150 (mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
151 const mpw *bv, const mpw *bvl, const mpw *mv,
152 size_t n, const mpw *mi),
153 (dv, dvl, av, avl, bv, bvl, mv, n, mi),
154 pick_mulcore, simple_mulcore);
156 static void simple_mulcore(mpw *dv, mpw *dvl,
157 const mpw *av, const mpw *avl,
158 const mpw *bv, const mpw *bvl,
159 const mpw *mv, size_t n, const mpw *mi)
161 mpw ai, b0, y, mi0 = *mi;
163 const mpw *mvl = mv + n;
166 /* --- Initial setup --- */
169 if (avl - av > bvl - bv) {
170 tv = av; av = bv; bv = tv;
171 tvl = avl; avl = bvl; bvl = tvl;
175 /* --- Multiply, until we run out of multiplicand --- */
177 while (i < n && av < avl) {
179 y = MPW((*dv + ai*b0)*mi0);
180 MPX_UMLAN(dv, dvl, bv, bvl, ai);
181 MPX_UMLAN(dv, dvl, mv, mvl, y);
185 /* --- Continue reducing until we run out of modulus --- */
189 MPX_UMLAN(dv, dvl, mv, mvl, y);
194 #define MAYBE_MUL4(impl) \
195 extern void mpxmont_mul4_##impl(mpw *dv, \
196 const mpw *av, const mpw *bv, \
198 size_t n, const mpw *mi); \
199 static void maybe_mul4_##impl(mpw *dv, mpw *dvl, \
200 const mpw *av, const mpw *avl, \
201 const mpw *bv, const mpw *bvl, \
202 const mpw *mv, size_t n, const mpw *mi) \
204 size_t an = avl - av, bn = bvl - bv; \
205 if (n%4 || an != n || bn != n) \
206 simple_mulcore(dv, dvl, av, avl, bv, bvl, mv, n, mi); \
208 mpxmont_mul4_##impl(dv, av, bv, mv, n, mi); \
209 MPX_ZERO(dv + 2*n + 1, dvl); \
219 MAYBE_MUL4(amd64_sse2)
220 MAYBE_MUL4(amd64_avx)
228 MAYBE_MUL4(arm64_simd)
231 static mulcore__functype *pick_mulcore(void)
234 DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_x86_avx,
235 cpu_feature_p(CPUFEAT_X86_AVX));
236 DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_x86_sse2,
237 cpu_feature_p(CPUFEAT_X86_SSE2));
240 DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_amd64_avx,
241 cpu_feature_p(CPUFEAT_X86_AVX));
242 DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_amd64_sse2,
243 cpu_feature_p(CPUFEAT_X86_SSE2));
246 DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_arm_neon,
247 cpu_feature_p(CPUFEAT_ARM_NEON));
250 DISPATCH_PICK_COND(mpmont_mul, maybe_mul4_arm64_simd,
251 cpu_feature_p(CPUFEAT_ARM_NEON));
253 DISPATCH_PICK_FALLBACK(mpmont_mul, simple_mulcore);
256 /* --- @finish@ --- *
258 * Arguments: @const mpmont *mm@ = pointer to a Montgomery reduction
260 * *mp *d@ = pointer to mostly-reduced operand
264 * Use: Applies the finishing touches to Montgomery reduction. The
265 * operand @d@ is a multiple of %$R%$ at this point, so it needs
266 * to be shifted down; the result might need a further
267 * subtraction to get it into the right interval; and we may
268 * need to do an additional subtraction if %$d$% is negative.
271 static void finish(const mpmont *mm, mp *d)
273 mpw *dv = d->v, *dvl = d->vl;
276 memmove(dv, dv + n, MPWS(dvl - (dv + n)));
279 if (MPX_UCMP(dv, dvl, >=, mm->m->v, mm->m->vl))
280 mpx_usub(dv, dvl, dv, dvl, mm->m->v, mm->m->vl);
283 mpx_usub(dv, dvl, mm->m->v, mm->m->vl, dv, dvl);
293 /*----- Reduction and multiplication --------------------------------------*/
295 /* --- @mpmont_create@ --- *
297 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
298 * @mp *m@ = modulus to use
300 * Returns: Zero on success, nonzero on error.
302 * Use: Initializes a Montgomery reduction context ready for use.
303 * The argument @m@ must be a positive odd integer.
306 #ifdef MPMONT_DISABLE
308 int mpmont_create(mpmont *mm, mp *m)
320 int mpmont_create(mpmont *mm, mp *m)
322 size_t n = MP_LEN(m);
323 mp *r2 = mp_new(2 * n + 1, 0);
326 /* --- Take a copy of the modulus --- */
328 if (!MP_POSP(m) || !MP_ODDP(m))
332 /* --- Determine %$R^2$% --- */
335 MPX_ZERO(r2->v, r2->vl - 1);
338 /* --- Find the magic value @mi@ --- */
340 mp_build(&r, r2->v + n, r2->vl);
341 mm->mi = mp_modinv(MP_NEW, m, &r);
342 mm->mi = mp_sub(mm->mi, &r, mm->mi);
343 MP_ENSURE(mm->mi, n);
345 /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */
348 mp_div(0, &mm->r2, r2, m);
349 mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
356 /* --- @mpmont_destroy@ --- *
358 * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
362 * Use: Disposes of a context when it's no longer of any use to
366 void mpmont_destroy(mpmont *mm)
374 /* --- @mpmont_reduce@ --- *
376 * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context
377 * @mp *d@ = destination
378 * @mp *a@ = source, assumed positive
380 * Returns: Result, %$a R^{-1} \bmod m$%.
383 #ifdef MPMONT_DISABLE
385 mp *mpmont_reduce(const mpmont *mm, mp *d, mp *a)
387 mp_div(0, &d, a, mm->m);
393 mp *mpmont_reduce(const mpmont *mm, mp *d, mp *a)
397 /* --- Check for serious Karatsuba reduction --- */
399 if (n > MPMONT_KTHRESH) {
404 if (MP_LEN(a) >= n) vl = a->v + n;
406 mp_build(&al, a->v, vl);
407 u = mp_mul(MP_NEW, &al, mm->mi);
408 if (MP_LEN(u) > n) u->vl = u->v + n;
409 u = mp_mul(u, u, mm->m);
415 /* --- Otherwise do it the hard way --- */
421 MP_DEST(d, 2*mm->n + 1, a->f);
422 redccore(d->v, d->vl, mm->m->v, mm->n, mm->mi->v);
425 /* --- Wrap everything up --- */
433 /* --- @mpmont_mul@ --- *
435 * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context
436 * @mp *d@ = destination
437 * @mp *a, *b@ = sources, assumed positive
439 * Returns: Result, %$a b R^{-1} \bmod m$%.
442 #ifdef MPMONT_DISABLE
444 mp *mpmont_mul(const mpmont *mm, mp *d, mp *a, mp *b)
447 mp_div(0, &d, d, mm->m);
453 mp *mpmont_mul(const mpmont *mm, mp *d, mp *a, mp *b)
457 if (n > MPMONT_KTHRESH) {
459 d = mpmont_reduce(mm, d, d);
461 a = MP_COPY(a); b = MP_COPY(b);
462 MP_DEST(d, 2*n + 1, a->f | b->f | MP_UNDEF);
463 mulcore(d->v, d->vl, a->v, a->vl, b->v, b->vl,
464 mm->m->v, mm->n, mm->mi->v);
465 d->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ b->f) & MP_NEG);
467 MP_DROP(a); MP_DROP(b);
475 /*----- Test rig ----------------------------------------------------------*/
479 #ifdef ENABLE_ASM_DEBUG
480 # include "regdump.h"
483 static int tcreate(dstr *v)
485 mp *m = *(mp **)v[0].buf;
486 mp *mi = *(mp **)v[1].buf;
487 mp *r = *(mp **)v[2].buf;
488 mp *r2 = *(mp **)v[3].buf;
493 mpmont_create(&mm, m);
495 if (mm.mi->v[0] != mi->v[0]) {
496 fprintf(stderr, "\n*** bad mi: found %lu, expected %lu",
497 (unsigned long)mm.mi->v[0], (unsigned long)mi->v[0]);
498 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
503 if (!MP_EQ(mm.r, r)) {
504 fputs("\n*** bad r", stderr);
505 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
506 fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10);
507 fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10);
512 if (!MP_EQ(mm.r2, r2)) {
513 fputs("\n*** bad r2", stderr);
514 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
515 fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10);
516 fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10);
526 assert(mparena_count(MPARENA_GLOBAL) == 0);
530 static int tmul(dstr *v)
532 mp *m = *(mp **)v[0].buf;
533 mp *a = *(mp **)v[1].buf;
534 mp *b = *(mp **)v[2].buf;
535 mp *r = *(mp **)v[3].buf;
539 mpmont_create(&mm, m);
542 mp *qr = mp_mul(MP_NEW, a, b);
543 mp_div(0, &qr, qr, m);
546 fputs("\n*** classical modmul failed", stderr);
547 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
548 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
549 fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
550 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
551 fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10);
560 mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
561 mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
562 mp *mr = mpmont_mul(&mm, MP_NEW, ar, br);
563 mr = mpmont_reduce(&mm, mr, mr);
565 fputs("\n*** montgomery modmul failed", stderr);
566 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
567 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
568 fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
569 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
570 fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
574 MP_DROP(ar); MP_DROP(br);
583 assert(mparena_count(MPARENA_GLOBAL) == 0);
587 static test_chunk tests[] = {
588 { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
589 { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
593 int main(int argc, char *argv[])
596 #ifdef ENABLE_ASM_DEBUG
599 test_run(argc, argv, tests, SRCDIR "/t/mpmont");
605 /*----- That's all, folks -------------------------------------------------*/