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
, 1);
130 DISPATCH_PICK_FALLBACK(mpmont_reduce
, simple_redccore
);
133 /* --- @mulcore@ --- *
135 * Arguments: @mpw *dv, *dvl@ = base and limit of source/destination
136 * @const mpw *av, *avl@ = base and limit of first multiplicand
137 * @const mpw *bv, *bvl@ = base and limit of second multiplicand
138 * @const mpw *mv@ = base of modulus %$m$%
139 * @size_t n@ = length of modulus
140 * @const mpw *mi@ = base of REDC coefficient %$m'$%
144 * Use: Let %$a$% and %$b$% be the multiplicands. Let %$w = a b$%.
145 * Store in %$d$% the value %$a b + (m' a b \bmod R) m$%.
148 CPU_DISPATCH(static, (void), void, mulcore
,
149 (mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
150 const mpw
*bv
, const mpw
*bvl
, const mpw
*mv
,
151 size_t n
, const mpw
*mi
),
152 (dv
, dvl
, av
, avl
, bv
, bvl
, mv
, n
, mi
),
153 pick_mulcore
, simple_mulcore
);
155 static void simple_mulcore(mpw
*dv
, mpw
*dvl
,
156 const mpw
*av
, const mpw
*avl
,
157 const mpw
*bv
, const mpw
*bvl
,
158 const mpw
*mv
, size_t n
, const mpw
*mi
)
160 mpw ai
, b0
, y
, mi0
= *mi
;
162 const mpw
*mvl
= mv
+ n
;
165 /* --- Initial setup --- */
168 if (avl
- av
> bvl
- bv
) {
169 tv
= av
; av
= bv
; bv
= tv
;
170 tvl
= avl
; avl
= bvl
; bvl
= tvl
;
174 /* --- Multiply, until we run out of multiplicand --- */
176 while (i
< n
&& av
< avl
) {
178 y
= MPW((*dv
+ ai
*b0
)*mi0
);
179 MPX_UMLAN(dv
, dvl
, bv
, bvl
, ai
);
180 MPX_UMLAN(dv
, dvl
, mv
, mvl
, y
);
184 /* --- Continue reducing until we run out of modulus --- */
188 MPX_UMLAN(dv
, dvl
, mv
, mvl
, y
);
193 #define MAYBE_MUL4(impl) \
194 extern void mpxmont_mul4_##impl(mpw *dv, \
195 const mpw *av, const mpw *bv, \
197 size_t n, const mpw *mi); \
198 static void maybe_mul4_##impl(mpw *dv, mpw *dvl, \
199 const mpw *av, const mpw *avl, \
200 const mpw *bv, const mpw *bvl, \
201 const mpw *mv, size_t n, const mpw *mi) \
203 size_t an = avl - av, bn = bvl - bv; \
204 if (n%4 || an != n || bn != n) \
205 simple_mulcore(dv, dvl, av, avl, bv, bvl, mv, n, mi); \
207 mpxmont_mul4_##impl(dv, av, bv, mv, n, mi); \
208 MPX_ZERO(dv + 2*n + 1, dvl); \
218 MAYBE_MUL4(amd64_sse2
)
219 MAYBE_MUL4(amd64_avx
)
227 MAYBE_MUL4(arm64_simd
)
230 static mulcore__functype
*pick_mulcore(void)
233 DISPATCH_PICK_COND(mpmont_mul
, maybe_mul4_x86_avx
,
234 cpu_feature_p(CPUFEAT_X86_AVX
));
235 DISPATCH_PICK_COND(mpmont_mul
, maybe_mul4_x86_sse2
,
236 cpu_feature_p(CPUFEAT_X86_SSE2
));
239 DISPATCH_PICK_COND(mpmont_mul
, maybe_mul4_amd64_avx
,
240 cpu_feature_p(CPUFEAT_X86_AVX
));
241 DISPATCH_PICK_COND(mpmont_mul
, maybe_mul4_amd64_sse2
,
242 cpu_feature_p(CPUFEAT_X86_SSE2
));
245 DISPATCH_PICK_COND(mpmont_mul
, maybe_mul4_arm_neon
,
246 cpu_feature_p(CPUFEAT_ARM_NEON
));
249 DISPATCH_PICK_COND(mpmont_mul
, maybe_mul4_arm64_simd
, 1);
251 DISPATCH_PICK_FALLBACK(mpmont_mul
, simple_mulcore
);
254 /* --- @finish@ --- *
256 * Arguments: @const mpmont *mm@ = pointer to a Montgomery reduction
258 * *mp *d@ = pointer to mostly-reduced operand
262 * Use: Applies the finishing touches to Montgomery reduction. The
263 * operand @d@ is a multiple of %$R%$ at this point, so it needs
264 * to be shifted down; the result might need a further
265 * subtraction to get it into the right interval; and we may
266 * need to do an additional subtraction if %$d$% is negative.
269 static void finish(const mpmont
*mm
, mp
*d
)
271 mpw
*dv
= d
->v
, *dvl
= d
->vl
;
274 memmove(dv
, dv
+ n
, MPWS(dvl
- (dv
+ n
)));
277 if (MPX_UCMP(dv
, dvl
, >=, mm
->m
->v
, mm
->m
->vl
))
278 mpx_usub(dv
, dvl
, dv
, dvl
, mm
->m
->v
, mm
->m
->vl
);
281 mpx_usub(dv
, dvl
, mm
->m
->v
, mm
->m
->vl
, dv
, dvl
);
291 /*----- Reduction and multiplication --------------------------------------*/
293 /* --- @mpmont_create@ --- *
295 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
296 * @mp *m@ = modulus to use
298 * Returns: Zero on success, nonzero on error.
300 * Use: Initializes a Montgomery reduction context ready for use.
301 * The argument @m@ must be a positive odd integer.
304 #ifdef MPMONT_DISABLE
306 int mpmont_create(mpmont
*mm
, mp
*m
)
318 int mpmont_create(mpmont
*mm
, mp
*m
)
320 size_t n
= MP_LEN(m
);
321 mp
*r2
= mp_new(2 * n
+ 1, 0);
324 /* --- Take a copy of the modulus --- */
326 if (!MP_POSP(m
) || !MP_ODDP(m
))
330 /* --- Determine %$R^2$% --- */
333 MPX_ZERO(r2
->v
, r2
->vl
- 1);
336 /* --- Find the magic value @mi@ --- */
338 mp_build(&r
, r2
->v
+ n
, r2
->vl
);
339 mm
->mi
= mp_modinv(MP_NEW
, m
, &r
);
340 mm
->mi
= mp_sub(mm
->mi
, &r
, mm
->mi
);
341 MP_ENSURE(mm
->mi
, n
);
343 /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */
346 mp_div(0, &mm
->r2
, r2
, m
);
347 mm
->r
= mpmont_reduce(mm
, MP_NEW
, mm
->r2
);
354 /* --- @mpmont_destroy@ --- *
356 * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
360 * Use: Disposes of a context when it's no longer of any use to
364 void mpmont_destroy(mpmont
*mm
)
372 /* --- @mpmont_reduce@ --- *
374 * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context
375 * @mp *d@ = destination
376 * @mp *a@ = source, assumed positive
378 * Returns: Result, %$a R^{-1} \bmod m$%.
381 #ifdef MPMONT_DISABLE
383 mp
*mpmont_reduce(const mpmont
*mm
, mp
*d
, mp
*a
)
385 mp_div(0, &d
, a
, mm
->m
);
391 mp
*mpmont_reduce(const mpmont
*mm
, mp
*d
, mp
*a
)
395 /* --- Check for serious Karatsuba reduction --- */
397 if (n
> MPMONT_KTHRESH
) {
402 if (MP_LEN(a
) >= n
) vl
= a
->v
+ n
;
404 mp_build(&al
, a
->v
, vl
);
405 u
= mp_mul(MP_NEW
, &al
, mm
->mi
);
406 if (MP_LEN(u
) > n
) u
->vl
= u
->v
+ n
;
407 u
= mp_mul(u
, u
, mm
->m
);
413 /* --- Otherwise do it the hard way --- */
419 MP_DEST(d
, 2*mm
->n
+ 1, a
->f
);
420 redccore(d
->v
, d
->vl
, mm
->m
->v
, mm
->n
, mm
->mi
->v
);
423 /* --- Wrap everything up --- */
431 /* --- @mpmont_mul@ --- *
433 * Arguments: @const mpmont *mm@ = pointer to Montgomery reduction context
434 * @mp *d@ = destination
435 * @mp *a, *b@ = sources, assumed positive
437 * Returns: Result, %$a b R^{-1} \bmod m$%.
440 #ifdef MPMONT_DISABLE
442 mp
*mpmont_mul(const mpmont
*mm
, mp
*d
, mp
*a
, mp
*b
)
445 mp_div(0, &d
, d
, mm
->m
);
451 mp
*mpmont_mul(const mpmont
*mm
, mp
*d
, mp
*a
, mp
*b
)
455 if (n
> MPMONT_KTHRESH
) {
457 d
= mpmont_reduce(mm
, d
, d
);
459 a
= MP_COPY(a
); b
= MP_COPY(b
);
460 MP_DEST(d
, 2*n
+ 1, a
->f
| b
->f
| MP_UNDEF
);
461 mulcore(d
->v
, d
->vl
, a
->v
, a
->vl
, b
->v
, b
->vl
,
462 mm
->m
->v
, mm
->n
, mm
->mi
->v
);
463 d
->f
= ((a
->f
| b
->f
) & MP_BURN
) | ((a
->f
^ b
->f
) & MP_NEG
);
465 MP_DROP(a
); MP_DROP(b
);
473 /*----- Test rig ----------------------------------------------------------*/
477 #ifdef ENABLE_ASM_DEBUG
478 # include "regdump.h"
481 static int tcreate(dstr
*v
)
483 mp
*m
= *(mp
**)v
[0].buf
;
484 mp
*mi
= *(mp
**)v
[1].buf
;
485 mp
*r
= *(mp
**)v
[2].buf
;
486 mp
*r2
= *(mp
**)v
[3].buf
;
491 mpmont_create(&mm
, m
);
493 if (mm
.mi
->v
[0] != mi
->v
[0]) {
494 fprintf(stderr
, "\n*** bad mi: found %lu, expected %lu",
495 (unsigned long)mm
.mi
->v
[0], (unsigned long)mi
->v
[0]);
496 fputs("\nm = ", stderr
); mp_writefile(m
, stderr
, 10);
501 if (!MP_EQ(mm
.r
, r
)) {
502 fputs("\n*** bad r", stderr
);
503 fputs("\nm = ", stderr
); mp_writefile(m
, stderr
, 10);
504 fputs("\nexpected ", stderr
); mp_writefile(r
, stderr
, 10);
505 fputs("\n found ", stderr
); mp_writefile(mm
.r
, stderr
, 10);
510 if (!MP_EQ(mm
.r2
, r2
)) {
511 fputs("\n*** bad r2", stderr
);
512 fputs("\nm = ", stderr
); mp_writefile(m
, stderr
, 10);
513 fputs("\nexpected ", stderr
); mp_writefile(r2
, stderr
, 10);
514 fputs("\n found ", stderr
); mp_writefile(mm
.r2
, stderr
, 10);
524 assert(mparena_count(MPARENA_GLOBAL
) == 0);
528 static int tmul(dstr
*v
)
530 mp
*m
= *(mp
**)v
[0].buf
;
531 mp
*a
= *(mp
**)v
[1].buf
;
532 mp
*b
= *(mp
**)v
[2].buf
;
533 mp
*r
= *(mp
**)v
[3].buf
;
537 mpmont_create(&mm
, m
);
540 mp
*qr
= mp_mul(MP_NEW
, a
, b
);
541 mp_div(0, &qr
, qr
, m
);
544 fputs("\n*** classical modmul failed", stderr
);
545 fputs("\n m = ", stderr
); mp_writefile(m
, stderr
, 10);
546 fputs("\n a = ", stderr
); mp_writefile(a
, stderr
, 10);
547 fputs("\n b = ", stderr
); mp_writefile(b
, stderr
, 10);
548 fputs("\n r = ", stderr
); mp_writefile(r
, stderr
, 10);
549 fputs("\nqr = ", stderr
); mp_writefile(qr
, stderr
, 10);
558 mp
*ar
= mpmont_mul(&mm
, MP_NEW
, a
, mm
.r2
);
559 mp
*br
= mpmont_mul(&mm
, MP_NEW
, b
, mm
.r2
);
560 mp
*mr
= mpmont_mul(&mm
, MP_NEW
, ar
, br
);
561 mr
= mpmont_reduce(&mm
, mr
, mr
);
563 fputs("\n*** montgomery modmul failed", stderr
);
564 fputs("\n m = ", stderr
); mp_writefile(m
, stderr
, 10);
565 fputs("\n a = ", stderr
); mp_writefile(a
, stderr
, 10);
566 fputs("\n b = ", stderr
); mp_writefile(b
, stderr
, 10);
567 fputs("\n r = ", stderr
); mp_writefile(r
, stderr
, 10);
568 fputs("\nmr = ", stderr
); mp_writefile(mr
, stderr
, 10);
572 MP_DROP(ar
); MP_DROP(br
);
581 assert(mparena_count(MPARENA_GLOBAL
) == 0);
585 static test_chunk tests
[] = {
586 { "create", tcreate
, { &type_mp
, &type_mp
, &type_mp
, &type_mp
, 0 } },
587 { "mul", tmul
, { &type_mp
, &type_mp
, &type_mp
, &type_mp
, 0 } },
591 int main(int argc
, char *argv
[])
594 #ifdef ENABLE_ASM_DEBUG
597 test_run(argc
, argv
, tests
, SRCDIR
"/t/mpmont");
603 /*----- That's all, folks -------------------------------------------------*/