Dual 1024-bit exponentiation optimization for Intel IceLake CPU

with AVX512_IFMA + AVX512_VL instructions, primarily for RSA CRT private key
operations. It uses 256-bit registers to avoid CPU frequency scaling issues.
The performance speedup for RSA2k signature on ICL is ~2x.

Reviewed-by: Paul Dale <pauli@openssl.org>
Reviewed-by: Matt Caswell <matt@openssl.org>
(Merged from https://github.com/openssl/openssl/pull/13750)
This commit is contained in:
Andrey Matyukov 2020-12-08 22:53:39 +03:00 committed by Matt Caswell
parent db89d8f04b
commit c781eb1c63
14 changed files with 1576 additions and 13 deletions

View File

@ -178,6 +178,11 @@ OpenSSL 3.0
*Tomáš Mráz*
* Parallel dual-prime 1024-bit modular exponentiation for AVX512_IFMA
capable processors.
*Ilya Albrekht, Sergey Kirillov, Andrey Matyukov (Intel Corp)*
* Combining the Configure options no-ec and no-dh no longer disables TLSv1.3.
Typically if OpenSSL has no EC or DH algorithms then it cannot support
connections with TLSv1.3. However OpenSSL now supports "pluggable" groups

View File

@ -0,0 +1,743 @@
# Copyright 2020 The OpenSSL Project Authors. All Rights Reserved.
# Copyright (c) 2020, Intel Corporation. All Rights Reserved.
#
# Licensed under the Apache License 2.0 (the "License"). You may not use
# this file except in compliance with the License. You can obtain a copy
# in the file LICENSE in the source distribution or at
# https://www.openssl.org/source/license.html
#
#
# Originally written by Ilya Albrekht, Sergey Kirillov and Andrey Matyukov
# Intel Corporation
#
# December 2020
#
# Initial release.
#
# Implementation utilizes 256-bit (ymm) registers to avoid frequency scaling issues.
#
# IceLake-Client @ 1.3GHz
# |---------+----------------------+--------------+-------------|
# | | OpenSSL 3.0.0-alpha9 | this | Unit |
# |---------+----------------------+--------------+-------------|
# | rsa2048 | 2 127 659 | 1 015 625 | cycles/sign |
# | | 611 | 1280 / +109% | sign/s |
# |---------+----------------------+--------------+-------------|
#
# $output is the last argument if it looks like a file (it has an extension)
# $flavour is the first argument if it doesn't look like a file
$output = $#ARGV >= 0 && $ARGV[$#ARGV] =~ m|\.\w+$| ? pop : undef;
$flavour = $#ARGV >= 0 && $ARGV[0] !~ m|\.| ? shift : undef;
$win64=0; $win64=1 if ($flavour =~ /[nm]asm|mingw64/ || $output =~ /\.asm$/);
$avx512ifma=0;
$0 =~ m/(.*[\/\\])[^\/\\]+$/; $dir=$1;
( $xlate="${dir}x86_64-xlate.pl" and -f $xlate ) or
( $xlate="${dir}../../perlasm/x86_64-xlate.pl" and -f $xlate) or
die "can't locate x86_64-xlate.pl";
if (`$ENV{CC} -Wa,-v -c -o /dev/null -x assembler /dev/null 2>&1`
=~ /GNU assembler version ([2-9]\.[0-9]+)/) {
$avx512ifma = ($1>=2.26);
}
if (!$avx512 && $win64 && ($flavour =~ /nasm/ || $ENV{ASM} =~ /nasm/) &&
`nasm -v 2>&1` =~ /NASM version ([2-9]\.[0-9]+)(?:\.([0-9]+))?/) {
$avx512ifma = ($1==2.11 && $2>=8) + ($1>=2.12);
}
if (!$avx512 && `$ENV{CC} -v 2>&1` =~ /((?:clang|LLVM) version|.*based on LLVM) ([0-9]+\.[0-9]+)/) {
$avx512ifma = ($2>=6.0);
}
open OUT,"| \"$^X\" \"$xlate\" $flavour \"$output\""
or die "can't call $xlate: $!";
*STDOUT=*OUT;
if ($avx512ifma>0) {{{
@_6_args_universal_ABI = ("%rdi","%rsi","%rdx","%rcx","%r8","%r9");
$code.=<<___;
.extern OPENSSL_ia32cap_P
.globl rsaz_avx512ifma_eligible
.type rsaz_avx512ifma_eligible,\@abi-omnipotent
.align 32
rsaz_avx512ifma_eligible:
mov OPENSSL_ia32cap_P+8(%rip), %ecx
xor %eax,%eax
and \$`1<<31|1<<21|1<<17|1<<16`, %ecx # avx512vl + avx512ifma + avx512dq + avx512f
cmp \$`1<<31|1<<21|1<<17|1<<16`, %ecx
cmove %ecx,%eax
ret
.size rsaz_avx512ifma_eligible, .-rsaz_avx512ifma_eligible
___
###############################################################################
# Almost Montgomery Multiplication (AMM) for 20-digit number in radix 2^52.
#
# AMM is defined as presented in the paper
# "Efficient Software Implementations of Modular Exponentiation" by Shay Gueron.
#
# The input and output are presented in 2^52 radix domain, i.e.
# |res|, |a|, |b|, |m| are arrays of 20 64-bit qwords with 12 high bits zeroed.
# |k0| is a Montgomery coefficient, which is here k0 = -1/m mod 2^64
# (note, the implementation counts only 52 bits from it).
#
# NB: the AMM implementation does not perform "conditional" subtraction step as
# specified in the original algorithm as according to the paper "Enhanced Montgomery
# Multiplication" by Shay Gueron (see Lemma 1), the result will be always < 2*2^1024
# and can be used as a direct input to the next AMM iteration.
# This post-condition is true, provided the correct parameter |s| is choosen, i.e.
# s >= n + 2 * k, which matches our case: 1040 > 1024 + 2 * 1.
#
# void RSAZ_amm52x20_x1_256(BN_ULONG *res,
# const BN_ULONG *a,
# const BN_ULONG *b,
# const BN_ULONG *m,
# BN_ULONG k0);
###############################################################################
{
# input parameters ("%rdi","%rsi","%rdx","%rcx","%r8")
my ($res,$a,$b,$m,$k0) = @_6_args_universal_ABI;
my $mask52 = "%rax";
my $acc0_0 = "%r9";
my $acc0_0_low = "%r9d";
my $acc0_1 = "%r15";
my $acc0_1_low = "%r15d";
my $b_ptr = "%r11";
my $iter = "%ebx";
my $zero = "%ymm0";
my ($R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0) = ("%ymm1", map("%ymm$_",(16..19)));
my ($R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1) = ("%ymm2", map("%ymm$_",(20..23)));
my $Bi = "%ymm3";
my $Yi = "%ymm4";
# Registers mapping for normalization.
# We can reuse Bi, Yi registers here.
my $TMP = $Bi;
my $mask52x4 = $Yi;
my ($T0,$T0h,$T1,$T1h,$T2) = map("%ymm$_", (24..28));
sub amm52x20_x1() {
# _data_offset - offset in the |a| or |m| arrays pointing to the beginning
# of data for corresponding AMM operation;
# _b_offset - offset in the |b| array pointing to the next qword digit;
my ($_data_offset,$_b_offset,$_acc,$_R0,$_R0h,$_R1,$_R1h,$_R2,$_k0) = @_;
my $_R0_xmm = $_R0 =~ s/%y/%x/r;
$code.=<<___;
movq $_b_offset($b_ptr), %r13 # b[i]
vpbroadcastq %r13, $Bi # broadcast b[i]
movq $_data_offset($a), %rdx
mulx %r13, %r13, %r12 # a[0]*b[i] = (t0,t2)
addq %r13, $_acc # acc += t0
movq %r12, %r10
adcq \$0, %r10 # t2 += CF
movq $_k0, %r13
imulq $_acc, %r13 # acc * k0
andq $mask52, %r13 # yi = (acc * k0) & mask52
vpbroadcastq %r13, $Yi # broadcast y[i]
movq $_data_offset($m), %rdx
mulx %r13, %r13, %r12 # yi * m[0] = (t0,t1)
addq %r13, $_acc # acc += t0
adcq %r12, %r10 # t2 += (t1 + CF)
shrq \$52, $_acc
salq \$12, %r10
or %r10, $_acc # acc = ((acc >> 52) | (t2 << 12))
vpmadd52luq `$_data_offset+64*0`($a), $Bi, $_R0
vpmadd52luq `$_data_offset+64*0+32`($a), $Bi, $_R0h
vpmadd52luq `$_data_offset+64*1`($a), $Bi, $_R1
vpmadd52luq `$_data_offset+64*1+32`($a), $Bi, $_R1h
vpmadd52luq `$_data_offset+64*2`($a), $Bi, $_R2
vpmadd52luq `$_data_offset+64*0`($m), $Yi, $_R0
vpmadd52luq `$_data_offset+64*0+32`($m), $Yi, $_R0h
vpmadd52luq `$_data_offset+64*1`($m), $Yi, $_R1
vpmadd52luq `$_data_offset+64*1+32`($m), $Yi, $_R1h
vpmadd52luq `$_data_offset+64*2`($m), $Yi, $_R2
# Shift accumulators right by 1 qword, zero extending the highest one
valignq \$1, $_R0, $_R0h, $_R0
valignq \$1, $_R0h, $_R1, $_R0h
valignq \$1, $_R1, $_R1h, $_R1
valignq \$1, $_R1h, $_R2, $_R1h
valignq \$1, $_R2, $zero, $_R2
vmovq $_R0_xmm, %r13
addq %r13, $_acc # acc += R0[0]
vpmadd52huq `$_data_offset+64*0`($a), $Bi, $_R0
vpmadd52huq `$_data_offset+64*0+32`($a), $Bi, $_R0h
vpmadd52huq `$_data_offset+64*1`($a), $Bi, $_R1
vpmadd52huq `$_data_offset+64*1+32`($a), $Bi, $_R1h
vpmadd52huq `$_data_offset+64*2`($a), $Bi, $_R2
vpmadd52huq `$_data_offset+64*0`($m), $Yi, $_R0
vpmadd52huq `$_data_offset+64*0+32`($m), $Yi, $_R0h
vpmadd52huq `$_data_offset+64*1`($m), $Yi, $_R1
vpmadd52huq `$_data_offset+64*1+32`($m), $Yi, $_R1h
vpmadd52huq `$_data_offset+64*2`($m), $Yi, $_R2
___
}
# Normalization routine: handles carry bits in R0..R2 QWs and
# gets R0..R2 back to normalized 2^52 representation.
#
# Uses %r8-14,%e[bcd]x
sub amm52x20_x1_norm {
my ($_acc,$_R0,$_R0h,$_R1,$_R1h,$_R2) = @_;
$code.=<<___;
# Put accumulator to low qword in R0
vpbroadcastq $_acc, $TMP
vpblendd \$3, $TMP, $_R0, $_R0
# Extract "carries" (12 high bits) from each QW of R0..R2
# Save them to LSB of QWs in T0..T2
vpsrlq \$52, $_R0, $T0
vpsrlq \$52, $_R0h, $T0h
vpsrlq \$52, $_R1, $T1
vpsrlq \$52, $_R1h, $T1h
vpsrlq \$52, $_R2, $T2
# "Shift left" T0..T2 by 1 QW
valignq \$3, $T1h, $T2, $T2
valignq \$3, $T1, $T1h, $T1h
valignq \$3, $T0h, $T1, $T1
valignq \$3, $T0, $T0h, $T0h
valignq \$3, $zero, $T0, $T0
# Drop "carries" from R0..R2 QWs
vpandq $mask52x4, $_R0, $_R0
vpandq $mask52x4, $_R0h, $_R0h
vpandq $mask52x4, $_R1, $_R1
vpandq $mask52x4, $_R1h, $_R1h
vpandq $mask52x4, $_R2, $_R2
# Sum R0..R2 with corresponding adjusted carries
vpaddq $T0, $_R0, $_R0
vpaddq $T0h, $_R0h, $_R0h
vpaddq $T1, $_R1, $_R1
vpaddq $T1h, $_R1h, $_R1h
vpaddq $T2, $_R2, $_R2
# Now handle carry bits from this addition
# Get mask of QWs which 52-bit parts overflow...
vpcmpuq \$1, $_R0, $mask52x4, %k1 # OP=lt
vpcmpuq \$1, $_R0h, $mask52x4, %k2
vpcmpuq \$1, $_R1, $mask52x4, %k3
vpcmpuq \$1, $_R1h, $mask52x4, %k4
vpcmpuq \$1, $_R2, $mask52x4, %k5
kmovb %k1, %r14d # k1
kmovb %k2, %r13d # k1h
kmovb %k3, %r12d # k2
kmovb %k4, %r11d # k2h
kmovb %k5, %r10d # k3
# ...or saturated
vpcmpuq \$0, $_R0, $mask52x4, %k1 # OP=eq
vpcmpuq \$0, $_R0h, $mask52x4, %k2
vpcmpuq \$0, $_R1, $mask52x4, %k3
vpcmpuq \$0, $_R1h, $mask52x4, %k4
vpcmpuq \$0, $_R2, $mask52x4, %k5
kmovb %k1, %r9d # k4
kmovb %k2, %r8d # k4h
kmovb %k3, %ebx # k5
kmovb %k4, %ecx # k5h
kmovb %k5, %edx # k6
# Get mask of QWs where carries shall be propagated to.
# Merge 4-bit masks to 8-bit values to use add with carry.
shl \$4, %r13b
or %r13b, %r14b
shl \$4, %r11b
or %r11b, %r12b
add %r14b, %r14b
adc %r12b, %r12b
adc %r10b, %r10b
shl \$4, %r8b
or %r8b,%r9b
shl \$4, %cl
or %cl, %bl
add %r9b, %r14b
adc %bl, %r12b
adc %dl, %r10b
xor %r9b, %r14b
xor %bl, %r12b
xor %dl, %r10b
kmovb %r14d, %k1
shr \$4, %r14b
kmovb %r14d, %k2
kmovb %r12d, %k3
shr \$4, %r12b
kmovb %r12d, %k4
kmovb %r10d, %k5
# Add carries according to the obtained mask
vpsubq $mask52x4, $_R0, ${_R0}{%k1}
vpsubq $mask52x4, $_R0h, ${_R0h}{%k2}
vpsubq $mask52x4, $_R1, ${_R1}{%k3}
vpsubq $mask52x4, $_R1h, ${_R1h}{%k4}
vpsubq $mask52x4, $_R2, ${_R2}{%k5}
vpandq $mask52x4, $_R0, $_R0
vpandq $mask52x4, $_R0h, $_R0h
vpandq $mask52x4, $_R1, $_R1
vpandq $mask52x4, $_R1h, $_R1h
vpandq $mask52x4, $_R2, $_R2
___
}
$code.=<<___;
.text
.globl RSAZ_amm52x20_x1_256
.type RSAZ_amm52x20_x1_256,\@function,5
.align 32
RSAZ_amm52x20_x1_256:
.cfi_startproc
endbranch
push %rbx
.cfi_push %rbx
push %rbp
.cfi_push %rbp
push %r12
.cfi_push %r12
push %r13
.cfi_push %r13
push %r14
.cfi_push %r14
push %r15
.cfi_push %r15
.Lrsaz_amm52x20_x1_256_body:
# Zeroing accumulators
vpxord $zero, $zero, $zero
vmovdqa64 $zero, $R0_0
vmovdqa64 $zero, $R0_0h
vmovdqa64 $zero, $R1_0
vmovdqa64 $zero, $R1_0h
vmovdqa64 $zero, $R2_0
xorl $acc0_0_low, $acc0_0_low
movq $b, $b_ptr # backup address of b
movq \$0xfffffffffffff, $mask52 # 52-bit mask
# Loop over 20 digits unrolled by 4
mov \$5, $iter
.align 32
.Lloop5:
___
foreach my $idx (0..3) {
&amm52x20_x1(0,8*$idx,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$k0);
}
$code.=<<___;
lea `4*8`($b_ptr), $b_ptr
dec $iter
jne .Lloop5
vmovdqa64 .Lmask52x4(%rip), $mask52x4
___
&amm52x20_x1_norm($acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0);
$code.=<<___;
vmovdqu64 $R0_0, ($res)
vmovdqu64 $R0_0h, 32($res)
vmovdqu64 $R1_0, 64($res)
vmovdqu64 $R1_0h, 96($res)
vmovdqu64 $R2_0, 128($res)
vzeroupper
mov 0(%rsp),%r15
.cfi_restore %r15
mov 8(%rsp),%r14
.cfi_restore %r14
mov 16(%rsp),%r13
.cfi_restore %r13
mov 24(%rsp),%r12
.cfi_restore %r12
mov 32(%rsp),%rbp
.cfi_restore %rbp
mov 40(%rsp),%rbx
.cfi_restore %rbx
lea 48(%rsp),%rsp
.cfi_adjust_cfa_offset -48
.Lrsaz_amm52x20_x1_256_epilogue:
ret
.cfi_endproc
.size RSAZ_amm52x20_x1_256, .-RSAZ_amm52x20_x1_256
___
$code.=<<___;
.data
.align 32
.Lmask52x4:
.quad 0xfffffffffffff
.quad 0xfffffffffffff
.quad 0xfffffffffffff
.quad 0xfffffffffffff
___
###############################################################################
# Dual Almost Montgomery Multiplication for 20-digit number in radix 2^52
#
# See description of RSAZ_amm52x20_x1_256() above for details about Almost
# Montgomery Multiplication algorithm and function input parameters description.
#
# This function does two AMMs for two independent inputs, hence dual.
#
# void RSAZ_amm52x20_x2_256(BN_ULONG out[2][20],
# const BN_ULONG a[2][20],
# const BN_ULONG b[2][20],
# const BN_ULONG m[2][20],
# const BN_ULONG k0[2]);
###############################################################################
$code.=<<___;
.text
.globl RSAZ_amm52x20_x2_256
.type RSAZ_amm52x20_x2_256,\@function,5
.align 32
RSAZ_amm52x20_x2_256:
.cfi_startproc
endbranch
push %rbx
.cfi_push %rbx
push %rbp
.cfi_push %rbp
push %r12
.cfi_push %r12
push %r13
.cfi_push %r13
push %r14
.cfi_push %r14
push %r15
.cfi_push %r15
.Lrsaz_amm52x20_x2_256_body:
# Zeroing accumulators
vpxord $zero, $zero, $zero
vmovdqa64 $zero, $R0_0
vmovdqa64 $zero, $R0_0h
vmovdqa64 $zero, $R1_0
vmovdqa64 $zero, $R1_0h
vmovdqa64 $zero, $R2_0
vmovdqa64 $zero, $R0_1
vmovdqa64 $zero, $R0_1h
vmovdqa64 $zero, $R1_1
vmovdqa64 $zero, $R1_1h
vmovdqa64 $zero, $R2_1
xorl $acc0_0_low, $acc0_0_low
xorl $acc0_1_low, $acc0_1_low
movq $b, $b_ptr # backup address of b
movq \$0xfffffffffffff, $mask52 # 52-bit mask
mov \$20, $iter
.align 32
.Lloop20:
___
&amm52x20_x1( 0, 0,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,"($k0)");
# 20*8 = offset of the next dimension in two-dimension array
&amm52x20_x1(20*8,20*8,$acc0_1,$R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1,"8($k0)");
$code.=<<___;
lea 8($b_ptr), $b_ptr
dec $iter
jne .Lloop20
vmovdqa64 .Lmask52x4(%rip), $mask52x4
___
&amm52x20_x1_norm($acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0);
&amm52x20_x1_norm($acc0_1,$R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1);
$code.=<<___;
vmovdqu64 $R0_0, ($res)
vmovdqu64 $R0_0h, 32($res)
vmovdqu64 $R1_0, 64($res)
vmovdqu64 $R1_0h, 96($res)
vmovdqu64 $R2_0, 128($res)
vmovdqu64 $R0_1, 160($res)
vmovdqu64 $R0_1h, 192($res)
vmovdqu64 $R1_1, 224($res)
vmovdqu64 $R1_1h, 256($res)
vmovdqu64 $R2_1, 288($res)
vzeroupper
mov 0(%rsp),%r15
.cfi_restore %r15
mov 8(%rsp),%r14
.cfi_restore %r14
mov 16(%rsp),%r13
.cfi_restore %r13
mov 24(%rsp),%r12
.cfi_restore %r12
mov 32(%rsp),%rbp
.cfi_restore %rbp
mov 40(%rsp),%rbx
.cfi_restore %rbx
lea 48(%rsp),%rsp
.cfi_adjust_cfa_offset -48
.Lrsaz_amm52x20_x2_256_epilogue:
ret
.cfi_endproc
.size RSAZ_amm52x20_x2_256, .-RSAZ_amm52x20_x2_256
___
}
###############################################################################
# Constant time extraction from the precomputed table of powers base^i, where
# i = 0..2^EXP_WIN_SIZE-1
#
# The input |red_table| contains precomputations for two independent base values,
# so the |tbl_idx| indicates for which base shall we extract the value.
# |red_table_idx| is a power index.
#
# Extracted value (output) is 20 digit number in 2^52 radix.
#
# void extract_multiplier_2x20_win5(BN_ULONG *red_Y,
# const BN_ULONG red_table[1 << EXP_WIN_SIZE][2][20],
# int red_table_idx,
# int tbl_idx); # 0 or 1
#
# EXP_WIN_SIZE = 5
###############################################################################
{
# input parameters
my ($out,$red_tbl,$red_tbl_idx,$tbl_idx) = @_6_args_universal_ABI;
my ($t0,$t1,$t2,$t3,$t4) = map("%ymm$_", (0..4));
my $t4xmm = $t4 =~ s/%y/%x/r;
my ($tmp0,$tmp1,$tmp2,$tmp3,$tmp4) = map("%ymm$_", (16..20));
my ($cur_idx,$idx,$ones) = map("%ymm$_", (21..23));
$code.=<<___;
.text
.align 32
.globl extract_multiplier_2x20_win5
.type extract_multiplier_2x20_win5,\@function,4
extract_multiplier_2x20_win5:
.cfi_startproc
endbranch
leaq ($tbl_idx,$tbl_idx,4), %rax
salq \$5, %rax
addq %rax, $red_tbl
vmovdqa64 .Lones(%rip), $ones # broadcast ones
vpbroadcastq $red_tbl_idx, $idx
leaq `(1<<5)*2*20*8`($red_tbl), %rax # holds end of the tbl
vpxor $t4xmm, $t4xmm, $t4xmm
vmovdqa64 $t4, $t3 # zeroing t0..4, cur_idx
vmovdqa64 $t4, $t2
vmovdqa64 $t4, $t1
vmovdqa64 $t4, $t0
vmovdqa64 $t4, $cur_idx
.align 32
.Lloop:
vpcmpq \$0, $cur_idx, $idx, %k1 # mask of (idx == cur_idx)
addq \$320, $red_tbl # 320 = 2 * 20 digits * 8 bytes
vpaddq $ones, $cur_idx, $cur_idx # increment cur_idx
vmovdqu64 -320($red_tbl), $tmp0 # load data from red_tbl
vmovdqu64 -288($red_tbl), $tmp1
vmovdqu64 -256($red_tbl), $tmp2
vmovdqu64 -224($red_tbl), $tmp3
vmovdqu64 -192($red_tbl), $tmp4
vpblendmq $tmp0, $t0, ${t0}{%k1} # extract data when mask is not zero
vpblendmq $tmp1, $t1, ${t1}{%k1}
vpblendmq $tmp2, $t2, ${t2}{%k1}
vpblendmq $tmp3, $t3, ${t3}{%k1}
vpblendmq $tmp4, $t4, ${t4}{%k1}
cmpq $red_tbl, %rax
jne .Lloop
vmovdqu64 $t0, ($out) # store t0..4
vmovdqu64 $t1, 32($out)
vmovdqu64 $t2, 64($out)
vmovdqu64 $t3, 96($out)
vmovdqu64 $t4, 128($out)
ret
.cfi_endproc
.size extract_multiplier_2x20_win5, .-extract_multiplier_2x20_win5
___
$code.=<<___;
.data
.align 32
.Lones:
.quad 1,1,1,1
___
}
if ($win64) {
$rec="%rcx";
$frame="%rdx";
$context="%r8";
$disp="%r9";
$code.=<<___
.extern __imp_RtlVirtualUnwind
.type rsaz_def_handler,\@abi-omnipotent
.align 16
rsaz_def_handler:
push %rsi
push %rdi
push %rbx
push %rbp
push %r12
push %r13
push %r14
push %r15
pushfq
sub \$64,%rsp
mov 120($context),%rax # pull context->Rax
mov 248($context),%rbx # pull context->Rip
mov 8($disp),%rsi # disp->ImageBase
mov 56($disp),%r11 # disp->HandlerData
mov 0(%r11),%r10d # HandlerData[0]
lea (%rsi,%r10),%r10 # prologue label
cmp %r10,%rbx # context->Rip<.Lprologue
jb .Lcommon_seh_tail
mov 152($context),%rax # pull context->Rsp
mov 4(%r11),%r10d # HandlerData[1]
lea (%rsi,%r10),%r10 # epilogue label
cmp %r10,%rbx # context->Rip>=.Lepilogue
jae .Lcommon_seh_tail
lea 48(%rax),%rax
mov -8(%rax),%rbx
mov -16(%rax),%rbp
mov -24(%rax),%r12
mov -32(%rax),%r13
mov -40(%rax),%r14
mov -48(%rax),%r15
mov %rbx,144($context) # restore context->Rbx
mov %rbp,160($context) # restore context->Rbp
mov %r12,216($context) # restore context->R12
mov %r13,224($context) # restore context->R13
mov %r14,232($context) # restore context->R14
mov %r15,240($context) # restore context->R14
.Lcommon_seh_tail:
mov 8(%rax),%rdi
mov 16(%rax),%rsi
mov %rax,152($context) # restore context->Rsp
mov %rsi,168($context) # restore context->Rsi
mov %rdi,176($context) # restore context->Rdi
mov 40($disp),%rdi # disp->ContextRecord
mov $context,%rsi # context
mov \$154,%ecx # sizeof(CONTEXT)
.long 0xa548f3fc # cld; rep movsq
mov $disp,%rsi
xor %rcx,%rcx # arg1, UNW_FLAG_NHANDLER
mov 8(%rsi),%rdx # arg2, disp->ImageBase
mov 0(%rsi),%r8 # arg3, disp->ControlPc
mov 16(%rsi),%r9 # arg4, disp->FunctionEntry
mov 40(%rsi),%r10 # disp->ContextRecord
lea 56(%rsi),%r11 # &disp->HandlerData
lea 24(%rsi),%r12 # &disp->EstablisherFrame
mov %r10,32(%rsp) # arg5
mov %r11,40(%rsp) # arg6
mov %r12,48(%rsp) # arg7
mov %rcx,56(%rsp) # arg8, (NULL)
call *__imp_RtlVirtualUnwind(%rip)
mov \$1,%eax # ExceptionContinueSearch
add \$64,%rsp
popfq
pop %r15
pop %r14
pop %r13
pop %r12
pop %rbp
pop %rbx
pop %rdi
pop %rsi
ret
.size rsaz_def_handler,.-rsaz_def_handler
.section .pdata
.align 4
.rva .LSEH_begin_RSAZ_amm52x20_x1_256
.rva .LSEH_end_RSAZ_amm52x20_x1_256
.rva .LSEH_info_RSAZ_amm52x20_x1_256
.rva .LSEH_begin_extract_multiplier_2x20_win5
.rva .LSEH_end_extract_multiplier_2x20_win5
.rva .LSEH_info_extract_multiplier_2x20_win5
.rva .LSEH_begin_RSAZ_amm52x20_x2_256
.rva .LSEH_end_RSAZ_amm52x20_x2_256
.rva .LSEH_info_RSAZ_amm52x20_x2_256
.section .xdata
.align 8
.LSEH_info_RSAZ_amm52x20_x1_256:
.byte 9,0,0,0
.rva rsaz_def_handler
.rva .Lrsaz_amm52x20_x1_256_body,.Lrsaz_amm52x20_x1_256_epilogue
.LSEH_info_extract_multiplier_2x20_win5:
.byte 9,0,0,0
.rva rsaz_def_handler
.rva .LSEH_begin_extract_multiplier_2x20_win5,.LSEH_begin_extract_multiplier_2x20_win5
.LSEH_info_RSAZ_amm52x20_x2_256:
.byte 9,0,0,0
.rva rsaz_def_handler
.rva .Lrsaz_amm52x20_x2_256_body,.Lrsaz_amm52x20_x2_256_epilogue
___
}
}}} else {{{ # fallback for old assembler
$code.=<<___;
.text
.globl rsaz_avx512ifma_eligible
.type rsaz_avx512ifma_eligible,\@abi-omnipotent
rsaz_avx512ifma_eligible:
xor %eax,%eax
ret
.size rsaz_avx512ifma_eligible, .-rsaz_avx512ifma_eligible
.globl RSAZ_amm52x20_x1_256
.globl RSAZ_amm52x20_x2_256
.globl extract_multiplier_2x20_win5
.type RSAZ_amm52x20_x1_256,\@abi-omnipotent
RSAZ_amm52x20_x1_256:
RSAZ_amm52x20_x2_256:
extract_multiplier_2x20_win5:
.byte 0x0f,0x0b # ud2
ret
.size RSAZ_amm52x20_x1_256, .-RSAZ_amm52x20_x1_256
___
}}}
$code =~ s/\`([^\`]*)\`/eval $1/gem;
print $code;
close STDOUT or die "error closing STDOUT: $!";

View File

@ -1390,3 +1390,85 @@ int BN_mod_exp_simple(BIGNUM *r, const BIGNUM *a, const BIGNUM *p,
bn_check_top(r);
return ret;
}
/*
* This is a variant of modular exponentiation optimization that does
* parallel 2-primes exponentiation using 256-bit (AVX512VL) AVX512_IFMA ISA
* in 52-bit binary redundant representation.
* If such instructions are not available, or input data size is not supported,
* it falls back to two BN_mod_exp_mont_consttime() calls.
*/
int BN_mod_exp_mont_consttime_x2(BIGNUM *rr1, const BIGNUM *a1, const BIGNUM *p1,
const BIGNUM *m1, BN_MONT_CTX *in_mont1,
BIGNUM *rr2, const BIGNUM *a2, const BIGNUM *p2,
const BIGNUM *m2, BN_MONT_CTX *in_mont2,
BN_CTX *ctx)
{
int ret = 0;
#ifdef RSAZ_ENABLED
BN_MONT_CTX *mont1 = NULL;
BN_MONT_CTX *mont2 = NULL;
if (rsaz_avx512ifma_eligible() &&
((a1->top == 16) && (p1->top == 16) && (BN_num_bits(m1) == 1024) &&
(a2->top == 16) && (p2->top == 16) && (BN_num_bits(m2) == 1024))) {
if (bn_wexpand(rr1, 16) == NULL)
goto err;
if (bn_wexpand(rr2, 16) == NULL)
goto err;
/* Ensure that montgomery contexts are initialized */
if (in_mont1 != NULL) {
mont1 = in_mont1;
} else {
if ((mont1 = BN_MONT_CTX_new()) == NULL)
goto err;
if (!BN_MONT_CTX_set(mont1, m1, ctx))
goto err;
}
if (in_mont2 != NULL) {
mont2 = in_mont2;
} else {
if ((mont2 = BN_MONT_CTX_new()) == NULL)
goto err;
if (!BN_MONT_CTX_set(mont2, m2, ctx))
goto err;
}
ret = RSAZ_mod_exp_avx512_x2(rr1->d, a1->d, p1->d, m1->d, mont1->RR.d,
mont1->n0[0],
rr2->d, a2->d, p2->d, m2->d, mont2->RR.d,
mont2->n0[0],
1024 /* factor bit size */);
rr1->top = 16;
rr1->neg = 0;
bn_correct_top(rr1);
bn_check_top(rr1);
rr2->top = 16;
rr2->neg = 0;
bn_correct_top(rr2);
bn_check_top(rr2);
goto err;
}
#endif
/* rr1 = a1^p1 mod m1 */
ret = BN_mod_exp_mont_consttime(rr1, a1, p1, m1, ctx, in_mont1);
/* rr2 = a2^p2 mod m2 */
ret &= BN_mod_exp_mont_consttime(rr2, a2, p2, m2, ctx, in_mont2);
#ifdef RSAZ_ENABLED
err:
if (in_mont2 == NULL)
BN_MONT_CTX_free(mont2);
if (in_mont1 == NULL)
BN_MONT_CTX_free(mont1);
#endif
return ret;
}

View File

@ -24,7 +24,7 @@ IF[{- !$disabled{asm} -}]
$BNASM_x86_64=\
x86_64-mont.s x86_64-mont5.s x86_64-gf2m.s rsaz_exp.c rsaz-x86_64.s \
rsaz-avx2.s
rsaz-avx2.s rsaz_exp_x2.c rsaz-avx512.s
IF[{- $config{target} !~ /^VC/ -}]
$BNASM_x86_64=asm/x86_64-gcc.c $BNASM_x86_64
ELSE
@ -154,6 +154,7 @@ GENERATE[x86_64-mont5.s]=asm/x86_64-mont5.pl
GENERATE[x86_64-gf2m.s]=asm/x86_64-gf2m.pl
GENERATE[rsaz-x86_64.s]=asm/rsaz-x86_64.pl
GENERATE[rsaz-avx2.s]=asm/rsaz-avx2.pl
GENERATE[rsaz-avx512.s]=asm/rsaz-avx512.pl
GENERATE[bn-ia64.s]=asm/ia64.S
GENERATE[ia64-mont.s]=asm/ia64-mont.pl

View File

@ -1,6 +1,6 @@
/*
* Copyright 2013-2018 The OpenSSL Project Authors. All Rights Reserved.
* Copyright (c) 2012, Intel Corporation. All Rights Reserved.
* Copyright 2013-2020 The OpenSSL Project Authors. All Rights Reserved.
* Copyright (c) 2020, Intel Corporation. All Rights Reserved.
*
* Licensed under the Apache License 2.0 (the "License"). You may not use
* this file except in compliance with the License. You can obtain a copy
@ -35,6 +35,23 @@ void RSAZ_512_mod_exp(BN_ULONG result[8],
const BN_ULONG m_norm[8], BN_ULONG k0,
const BN_ULONG RR[8]);
int rsaz_avx512ifma_eligible(void);
int RSAZ_mod_exp_avx512_x2(BN_ULONG *res1,
const BN_ULONG *base1,
const BN_ULONG *exponent1,
const BN_ULONG *m1,
const BN_ULONG *RR1,
BN_ULONG k0_1,
BN_ULONG *res2,
const BN_ULONG *base2,
const BN_ULONG *exponent2,
const BN_ULONG *m2,
const BN_ULONG *RR2,
BN_ULONG k0_2,
int factor_size);
# endif
#endif

542
crypto/bn/rsaz_exp_x2.c Normal file
View File

@ -0,0 +1,542 @@
/*
* Copyright 2020 The OpenSSL Project Authors. All Rights Reserved.
* Copyright (c) 2020, Intel Corporation. All Rights Reserved.
*
* Licensed under the Apache License 2.0 (the "License"). You may not use
* this file except in compliance with the License. You can obtain a copy
* in the file LICENSE in the source distribution or at
* https://www.openssl.org/source/license.html
*
*
* Originally written by Ilya Albrekht, Sergey Kirillov and Andrey Matyukov
* Intel Corporation
*
*/
#include <openssl/opensslconf.h>
#include "rsaz_exp.h"
#ifndef RSAZ_ENABLED
NON_EMPTY_TRANSLATION_UNIT
#else
# include <assert.h>
# include <string.h>
# if defined(__GNUC__)
# define ALIGN64 __attribute__((aligned(64)))
# elif defined(_MSC_VER)
# define ALIGN64 __declspec(align(64))
# else
# define ALIGN64
# endif
# define ALIGN_OF(ptr, boundary) \
((unsigned char *)(ptr) + (boundary - (((size_t)(ptr)) & (boundary - 1))))
/* Internal radix */
# define DIGIT_SIZE (52)
/* 52-bit mask */
# define DIGIT_MASK ((uint64_t)0xFFFFFFFFFFFFF)
# define BITS2WORD8_SIZE(x) (((x) + 7) >> 3)
# define BITS2WORD64_SIZE(x) (((x) + 63) >> 6)
static ossl_inline uint64_t get_digit52(const uint8_t *in, int in_len);
static ossl_inline void put_digit52(uint8_t *out, int out_len, uint64_t digit);
static void to_words52(BN_ULONG *out, int out_len, const BN_ULONG *in,
int in_bitsize);
static void from_words52(BN_ULONG *bn_out, int out_bitsize, const BN_ULONG *in);
static ossl_inline void set_bit(BN_ULONG *a, int idx);
/* Number of |digit_size|-bit digits in |bitsize|-bit value */
static ossl_inline int number_of_digits(int bitsize, int digit_size)
{
return (bitsize + digit_size - 1) / digit_size;
}
typedef void (*AMM52)(BN_ULONG *res, const BN_ULONG *base,
const BN_ULONG *exp, const BN_ULONG *m, BN_ULONG k0);
typedef void (*EXP52_x2)(BN_ULONG *res, const BN_ULONG *base,
const BN_ULONG *exp[2], const BN_ULONG *m,
const BN_ULONG *rr, const BN_ULONG k0[2]);
/*
* For details of the methods declared below please refer to
* crypto/bn/asm/rsaz-avx512.pl
*
* Naming notes:
* amm = Almost Montgomery Multiplication
* ams = Almost Montgomery Squaring
* 52x20 - data represented as array of 20 digits in 52-bit radix
* _x1_/_x2_ - 1 or 2 independent inputs/outputs
* _256 suffix - uses 256-bit (AVX512VL) registers
*/
/*AMM = Almost Montgomery Multiplication. */
void RSAZ_amm52x20_x1_256(BN_ULONG *res, const BN_ULONG *base,
const BN_ULONG *exp, const BN_ULONG *m,
BN_ULONG k0);
void RSAZ_exp52x20_x2_256(BN_ULONG *res, const BN_ULONG *base,
const BN_ULONG *exp[2], const BN_ULONG *m,
const BN_ULONG *rr, const BN_ULONG k0[2]);
void RSAZ_amm52x20_x2_256(BN_ULONG *out, const BN_ULONG *a,
const BN_ULONG *b, const BN_ULONG *m,
const BN_ULONG k0[2]);
void extract_multiplier_2x20_win5(BN_ULONG *red_Y,
const BN_ULONG *red_table,
int red_table_idx, int tbl_idx);
/*
* Dual Montgomery modular exponentiation using prime moduli of the
* same bit size, optimized with AVX512 ISA.
*
* Input and output parameters for each exponentiation are independent and
* denoted here by index |i|, i = 1..2.
*
* Input and output are all in regular 2^64 radix.
*
* Each moduli shall be |factor_size| bit size.
*
* NOTE: currently only 2x1024 case is supported.
*
* [out] res|i| - result of modular exponentiation: array of qword values
* in regular (2^64) radix. Size of array shall be enough
* to hold |factor_size| bits.
* [in] base|i| - base
* [in] exp|i| - exponent
* [in] m|i| - moduli
* [in] rr|i| - Montgomery parameter RR = R^2 mod m|i|
* [in] k0_|i| - Montgomery parameter k0 = -1/m|i| mod 2^64
* [in] factor_size - moduli bit size
*
* \return 0 in case of failure,
* 1 in case of success.
*/
int RSAZ_mod_exp_avx512_x2(BN_ULONG *res1,
const BN_ULONG *base1,
const BN_ULONG *exp1,
const BN_ULONG *m1,
const BN_ULONG *rr1,
BN_ULONG k0_1,
BN_ULONG *res2,
const BN_ULONG *base2,
const BN_ULONG *exp2,
const BN_ULONG *m2,
const BN_ULONG *rr2,
BN_ULONG k0_2,
int factor_size)
{
int ret = 0;
/*
* Number of word-size (BN_ULONG) digits to store exponent in redundant
* representation.
*/
int exp_digits = number_of_digits(factor_size + 2, DIGIT_SIZE);
int coeff_pow = 4 * (DIGIT_SIZE * exp_digits - factor_size);
BN_ULONG *base1_red, *m1_red, *rr1_red;
BN_ULONG *base2_red, *m2_red, *rr2_red;
BN_ULONG *coeff_red;
BN_ULONG *storage = NULL;
BN_ULONG *storage_aligned = NULL;
BN_ULONG storage_len_bytes = 7 * exp_digits * sizeof(BN_ULONG);
/* AMM = Almost Montgomery Multiplication */
AMM52 amm = NULL;
/* Dual (2-exps in parallel) exponentiation */
EXP52_x2 exp_x2 = NULL;
const BN_ULONG *exp[2] = {0};
BN_ULONG k0[2] = {0};
/* Only 1024-bit factor size is supported now */
switch (factor_size) {
case 1024:
amm = RSAZ_amm52x20_x1_256;
exp_x2 = RSAZ_exp52x20_x2_256;
break;
default:
goto err;
}
storage = (BN_ULONG *)OPENSSL_malloc(storage_len_bytes + 64);
if (storage == NULL)
goto err;
storage_aligned = (BN_ULONG *)ALIGN_OF(storage, 64);
/* Memory layout for red(undant) representations */
base1_red = storage_aligned;
base2_red = storage_aligned + 1 * exp_digits;
m1_red = storage_aligned + 2 * exp_digits;
m2_red = storage_aligned + 3 * exp_digits;
rr1_red = storage_aligned + 4 * exp_digits;
rr2_red = storage_aligned + 5 * exp_digits;
coeff_red = storage_aligned + 6 * exp_digits;
/* Convert base_i, m_i, rr_i, from regular to 52-bit radix */
to_words52(base1_red, exp_digits, base1, factor_size);
to_words52(base2_red, exp_digits, base2, factor_size);
to_words52(m1_red, exp_digits, m1, factor_size);
to_words52(m2_red, exp_digits, m2, factor_size);
to_words52(rr1_red, exp_digits, rr1, factor_size);
to_words52(rr2_red, exp_digits, rr2, factor_size);
/*
* Compute target domain Montgomery converters RR' for each modulus
* based on precomputed original domain's RR.
*
* RR -> RR' transformation steps:
* (1) coeff = 2^k
* (2) t = AMM(RR,RR) = RR^2 / R' mod m
* (3) RR' = AMM(t, coeff) = RR^2 * 2^k / R'^2 mod m
* where
* k = 4 * (52 * digits52 - modlen)
* R = 2^(64 * ceil(modlen/64)) mod m
* RR = R^2 mod M
* R' = 2^(52 * ceil(modlen/52)) mod m
*
* modlen = 1024: k = 64, RR = 2^2048 mod m, RR' = 2^2080 mod m
*/
memset(coeff_red, 0, exp_digits * sizeof(BN_ULONG));
/* (1) in reduced domain representation */
set_bit(coeff_red, 64 * (int)(coeff_pow / 52) + coeff_pow % 52);
amm(rr1_red, rr1_red, rr1_red, m1_red, k0_1); /* (2) for m1 */
amm(rr1_red, rr1_red, coeff_red, m1_red, k0_1); /* (3) for m1 */
amm(rr2_red, rr2_red, rr2_red, m2_red, k0_2); /* (2) for m2 */
amm(rr2_red, rr2_red, coeff_red, m2_red, k0_2); /* (3) for m2 */
exp[0] = exp1;
exp[1] = exp2;
k0[0] = k0_1;
k0[1] = k0_2;
exp_x2(rr1_red, base1_red, exp, m1_red, rr1_red, k0);
/* Convert rr_i back to regular radix */
from_words52(res1, factor_size, rr1_red);
from_words52(res2, factor_size, rr2_red);
ret = 1;
err:
if (storage != NULL) {
OPENSSL_cleanse(storage, storage_len_bytes);
OPENSSL_free(storage);
}
return ret;
}
/*
* Dual 1024-bit w-ary modular exponentiation using prime moduli of the same
* bit size using Almost Montgomery Multiplication, optimized with AVX512_IFMA
* ISA.
*
* The parameter w (window size) = 5.
*
* [out] res - result of modular exponentiation: 2x20 qword
* values in 2^52 radix.
* [in] base - base (2x20 qword values in 2^52 radix)
* [in] exp - array of 2 pointers to 16 qword values in 2^64 radix.
* Exponent is not converted to redundant representation.
* [in] m - moduli (2x20 qword values in 2^52 radix)
* [in] rr - Montgomery parameter for 2 moduli: RR = 2^2080 mod m.
* (2x20 qword values in 2^52 radix)
* [in] k0 - Montgomery parameter for 2 moduli: k0 = -1/m mod 2^64
*
* \return (void).
*/
void RSAZ_exp52x20_x2_256(BN_ULONG *out, /* [2][20] */
const BN_ULONG *base, /* [2][20] */
const BN_ULONG *exp[2], /* 2x16 */
const BN_ULONG *m, /* [2][20] */
const BN_ULONG *rr, /* [2][20] */
const BN_ULONG k0[2])
{
# define BITSIZE_MODULUS (1024)
# define EXP_WIN_SIZE (5)
# define EXP_WIN_MASK ((1U << EXP_WIN_SIZE) - 1)
/*
* Number of digits (64-bit words) in redundant representation to handle
* modulus bits
*/
# define RED_DIGITS (20)
# define EXP_DIGITS (16)
# define DAMM RSAZ_amm52x20_x2_256
/*
* Squaring is done using multiplication now. That can be a subject of
* optimization in future.
*/
# define DAMS(r,a,m,k0) \
RSAZ_amm52x20_x2_256((r),(a),(a),(m),(k0))
/* Allocate stack for red(undant) result Y and multiplier X */
ALIGN64 BN_ULONG red_Y[2][RED_DIGITS];
ALIGN64 BN_ULONG red_X[2][RED_DIGITS];
/* Allocate expanded exponent */
ALIGN64 BN_ULONG expz[2][EXP_DIGITS + 1];
/* Pre-computed table of base powers */
ALIGN64 BN_ULONG red_table[1U << EXP_WIN_SIZE][2][RED_DIGITS];
int idx;
memset(red_Y, 0, sizeof(red_Y));
memset(red_table, 0, sizeof(red_table));
memset(red_X, 0, sizeof(red_X));
/*
* Compute table of powers base^i, i = 0, ..., (2^EXP_WIN_SIZE) - 1
* table[0] = mont(x^0) = mont(1)
* table[1] = mont(x^1) = mont(x)
*/
red_X[0][0] = 1;
red_X[1][0] = 1;
DAMM(red_table[0][0], (const BN_ULONG*)red_X, rr, m, k0);
DAMM(red_table[1][0], base, rr, m, k0);
for (idx = 1; idx < (int)((1U << EXP_WIN_SIZE) / 2); idx++) {
DAMS(red_table[2 * idx + 0][0], red_table[1 * idx][0], m, k0);
DAMM(red_table[2 * idx + 1][0], red_table[2 * idx][0], red_table[1][0], m, k0);
}
/* Copy and expand exponents */
memcpy(expz[0], exp[0], EXP_DIGITS * sizeof(BN_ULONG));
expz[0][EXP_DIGITS] = 0;
memcpy(expz[1], exp[1], EXP_DIGITS * sizeof(BN_ULONG));
expz[1][EXP_DIGITS] = 0;
/* Exponentiation */
{
int rem = BITSIZE_MODULUS % EXP_WIN_SIZE;
int delta = rem ? rem : EXP_WIN_SIZE;
BN_ULONG table_idx_mask = EXP_WIN_MASK;
int exp_bit_no = BITSIZE_MODULUS - delta;
int exp_chunk_no = exp_bit_no / 64;
int exp_chunk_shift = exp_bit_no % 64;
/* Process 1-st exp window - just init result */
BN_ULONG red_table_idx_0 = expz[0][exp_chunk_no];
BN_ULONG red_table_idx_1 = expz[1][exp_chunk_no];
/*
* The function operates with fixed moduli sizes divisible by 64,
* thus table index here is always in supported range [0, EXP_WIN_SIZE).
*/
red_table_idx_0 >>= exp_chunk_shift;
red_table_idx_1 >>= exp_chunk_shift;
extract_multiplier_2x20_win5(red_Y[0], (const BN_ULONG*)red_table, (int)red_table_idx_0, 0);
extract_multiplier_2x20_win5(red_Y[1], (const BN_ULONG*)red_table, (int)red_table_idx_1, 1);
/* Process other exp windows */
for (exp_bit_no -= EXP_WIN_SIZE; exp_bit_no >= 0; exp_bit_no -= EXP_WIN_SIZE) {
/* Extract pre-computed multiplier from the table */
{
BN_ULONG T;
exp_chunk_no = exp_bit_no / 64;
exp_chunk_shift = exp_bit_no % 64;
{
red_table_idx_0 = expz[0][exp_chunk_no];
T = expz[0][exp_chunk_no + 1];
red_table_idx_0 >>= exp_chunk_shift;
/*
* Get additional bits from then next quadword
* when 64-bit boundaries are crossed.
*/
if (exp_chunk_shift > 64 - EXP_WIN_SIZE) {
T <<= (64 - exp_chunk_shift);
red_table_idx_0 ^= T;
}
red_table_idx_0 &= table_idx_mask;
extract_multiplier_2x20_win5(red_X[0], (const BN_ULONG*)red_table, (int)red_table_idx_0, 0);
}
{
red_table_idx_1 = expz[1][exp_chunk_no];
T = expz[1][exp_chunk_no + 1];
red_table_idx_1 >>= exp_chunk_shift;
/*
* Get additional bits from then next quadword
* when 64-bit boundaries are crossed.
*/
if (exp_chunk_shift > 64 - EXP_WIN_SIZE) {
T <<= (64 - exp_chunk_shift);
red_table_idx_1 ^= T;
}
red_table_idx_1 &= table_idx_mask;
extract_multiplier_2x20_win5(red_X[1], (const BN_ULONG*)red_table, (int)red_table_idx_1, 1);
}
}
/* Series of squaring */
DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
DAMM((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, (const BN_ULONG*)red_X, m, k0);
}
}
/*
*
* NB: After the last AMM of exponentiation in Montgomery domain, the result
* may be 1025-bit, but the conversion out of Montgomery domain performs an
* AMM(x,1) which guarantees that the final result is less than |m|, so no
* conditional subtraction is needed here. See "Efficient Software
* Implementations of Modular Exponentiation" (by Shay Gueron) paper for details.
*/
/* Convert result back in regular 2^52 domain */
memset(red_X, 0, sizeof(red_X));
red_X[0][0] = 1;
red_X[1][0] = 1;
DAMM(out, (const BN_ULONG*)red_Y, (const BN_ULONG*)red_X, m, k0);
/* Clear exponents */
OPENSSL_cleanse(expz, sizeof(expz));
OPENSSL_cleanse(red_Y, sizeof(red_Y));
# undef DAMS
# undef DAMM
# undef EXP_DIGITS
# undef RED_DIGITS
# undef EXP_WIN_MASK
# undef EXP_WIN_SIZE
# undef BITSIZE_MODULUS
}
static ossl_inline uint64_t get_digit52(const uint8_t *in, int in_len)
{
uint64_t digit = 0;
assert(in != NULL);
for (; in_len > 0; in_len--) {
digit <<= 8;
digit += (uint64_t)(in[in_len - 1]);
}
return digit;
}
/*
* Convert array of words in regular (base=2^64) representation to array of
* words in redundant (base=2^52) one.
*/
static void to_words52(BN_ULONG *out, int out_len,
const BN_ULONG *in, int in_bitsize)
{
uint8_t *in_str = NULL;
assert(out != NULL);
assert(in != NULL);
/* Check destination buffer capacity */
assert(out_len >= number_of_digits(in_bitsize, DIGIT_SIZE));
in_str = (uint8_t *)in;
for (; in_bitsize >= (2 * DIGIT_SIZE); in_bitsize -= (2 * DIGIT_SIZE), out += 2) {
out[0] = (*(uint64_t *)in_str) & DIGIT_MASK;
in_str += 6;
out[1] = ((*(uint64_t *)in_str) >> 4) & DIGIT_MASK;
in_str += 7;
out_len -= 2;
}
if (in_bitsize > DIGIT_SIZE) {
uint64_t digit = get_digit52(in_str, 7);
out[0] = digit & DIGIT_MASK;
in_str += 6;
in_bitsize -= DIGIT_SIZE;
digit = get_digit52(in_str, BITS2WORD8_SIZE(in_bitsize));
out[1] = digit >> 4;
out += 2;
out_len -= 2;
} else if (in_bitsize > 0) {
out[0] = get_digit52(in_str, BITS2WORD8_SIZE(in_bitsize));
out++;
out_len--;
}
while (out_len > 0) {
*out = 0;
out_len--;
out++;
}
}
static ossl_inline void put_digit52(uint8_t *pStr, int strLen, uint64_t digit)
{
assert(pStr != NULL);
for (; strLen > 0; strLen--) {
*pStr++ = (uint8_t)(digit & 0xFF);
digit >>= 8;
}
}
/*
* Convert array of words in redundant (base=2^52) representation to array of
* words in regular (base=2^64) one.
*/
static void from_words52(BN_ULONG *out, int out_bitsize, const BN_ULONG *in)
{
int i;
int out_len = BITS2WORD64_SIZE(out_bitsize);
assert(out != NULL);
assert(in != NULL);
for (i = 0; i < out_len; i++)
out[i] = 0;
{
uint8_t *out_str = (uint8_t *)out;
for (; out_bitsize >= (2 * DIGIT_SIZE); out_bitsize -= (2 * DIGIT_SIZE), in += 2) {
(*(uint64_t *)out_str) = in[0];
out_str += 6;
(*(uint64_t *)out_str) ^= in[1] << 4;
out_str += 7;
}
if (out_bitsize > DIGIT_SIZE) {
put_digit52(out_str, 7, in[0]);
out_str += 6;
out_bitsize -= DIGIT_SIZE;
put_digit52(out_str, BITS2WORD8_SIZE(out_bitsize),
(in[1] << 4 | in[0] >> 48));
} else if (out_bitsize) {
put_digit52(out_str, BITS2WORD8_SIZE(out_bitsize), in[0]);
}
}
}
/*
* Set bit at index |idx| in the words array |a|.
* It does not do any boundaries checks, make sure the index is valid before
* calling the function.
*/
static ossl_inline void set_bit(BN_ULONG *a, int idx)
{
assert(a != NULL);
{
int i, j;
i = idx / BN_BITS2;
j = idx % BN_BITS2;
a[i] |= (((BN_ULONG)1) << j);
}
}
#endif

View File

@ -688,15 +688,20 @@ static int rsa_ossl_mod_exp(BIGNUM *r0, const BIGNUM *I, RSA *rsa, BN_CTX *ctx)
if (/* m1 = I moq q */
!bn_from_mont_fixed_top(m1, I, rsa->_method_mod_q, ctx)
|| !bn_to_mont_fixed_top(m1, m1, rsa->_method_mod_q, ctx)
/* m1 = m1^dmq1 mod q */
|| !BN_mod_exp_mont_consttime(m1, m1, rsa->dmq1, rsa->q, ctx,
rsa->_method_mod_q)
/* r1 = I mod p */
|| !bn_from_mont_fixed_top(r1, I, rsa->_method_mod_p, ctx)
|| !bn_to_mont_fixed_top(r1, r1, rsa->_method_mod_p, ctx)
/* r1 = r1^dmp1 mod p */
|| !BN_mod_exp_mont_consttime(r1, r1, rsa->dmp1, rsa->p, ctx,
rsa->_method_mod_p)
/*
* Use parallel exponentiations optimization if possible,
* otherwise fallback to two sequential exponentiations:
* m1 = m1^dmq1 mod q
* r1 = r1^dmp1 mod p
*/
|| !BN_mod_exp_mont_consttime_x2(m1, m1, rsa->dmq1, rsa->q,
rsa->_method_mod_q,
r1, r1, rsa->dmp1, rsa->p,
rsa->_method_mod_p,
ctx)
/* r1 = (r1 - m1) mod p */
/*
* bn_mod_sub_fixed_top is not regular modular subtraction,

View File

@ -215,7 +215,7 @@ OPENSSL_ia32_cpuid:
cmp \$0xe6,%eax
je .Ldone
andl \$0x3fdeffff,8(%rdi) # ~(1<<31|1<<30|1<<21|1<<16)
# clear AVX512F+BW+VL+FIMA, all of
# clear AVX512F+BW+VL+IFMA, all of
# them are EVEX-encoded, which requires
# ZMM state support even if one uses
# only XMM and YMM :-(

View File

@ -698,6 +698,10 @@ DEPEND[html/man3/BN_generate_prime.html]=man3/BN_generate_prime.pod
GENERATE[html/man3/BN_generate_prime.html]=man3/BN_generate_prime.pod
DEPEND[man/man3/BN_generate_prime.3]=man3/BN_generate_prime.pod
GENERATE[man/man3/BN_generate_prime.3]=man3/BN_generate_prime.pod
DEPEND[html/man3/BN_mod_exp_mont.html]=man3/BN_mod_exp_mont.pod
GENERATE[html/man3/BN_mod_exp_mont.html]=man3/BN_mod_exp_mont.pod
DEPEND[man/man3/BN_mod_exp_mont.3]=man3/BN_mod_exp_mont.pod
GENERATE[man/man3/BN_mod_exp_mont.3]=man3/BN_mod_exp_mont.pod
DEPEND[html/man3/BN_mod_inverse.html]=man3/BN_mod_inverse.pod
GENERATE[html/man3/BN_mod_inverse.html]=man3/BN_mod_inverse.pod
DEPEND[man/man3/BN_mod_inverse.3]=man3/BN_mod_inverse.pod
@ -2808,6 +2812,7 @@ html/man3/BN_bn2bin.html \
html/man3/BN_cmp.html \
html/man3/BN_copy.html \
html/man3/BN_generate_prime.html \
html/man3/BN_mod_exp_mont.html \
html/man3/BN_mod_inverse.html \
html/man3/BN_mod_mul_montgomery.html \
html/man3/BN_mod_mul_reciprocal.html \
@ -3379,6 +3384,7 @@ man/man3/BN_bn2bin.3 \
man/man3/BN_cmp.3 \
man/man3/BN_copy.3 \
man/man3/BN_generate_prime.3 \
man/man3/BN_mod_exp_mont.3 \
man/man3/BN_mod_inverse.3 \
man/man3/BN_mod_mul_montgomery.3 \
man/man3/BN_mod_mul_reciprocal.3 \

View File

@ -0,0 +1,65 @@
=pod
=head1 NAME
BN_mod_exp_mont, BN_mod_exp_mont_consttime, BN_mod_exp_mont_consttime_x2 -
Montgomery exponentiation
=head1 SYNOPSIS
#include <openssl/bn.h>
int BN_mod_exp_mont(BIGNUM *rr, const BIGNUM *a, const BIGNUM *p,
const BIGNUM *m, BN_CTX *ctx, BN_MONT_CTX *in_mont);
int BN_mod_exp_mont_consttime(BIGNUM *rr, const BIGNUM *a, const BIGNUM *p,
const BIGNUM *m, BN_CTX *ctx,
BN_MONT_CTX *in_mont);
int BN_mod_exp_mont_consttime_x2(BIGNUM *rr1, const BIGNUM *a1,
const BIGNUM *p1, const BIGNUM *m1,
BN_MONT_CTX *in_mont1, BIGNUM *rr2,
const BIGNUM *a2, const BIGNUM *p2,
const BIGNUM *m2, BN_MONT_CTX *in_mont2,
BN_CTX *ctx);
=head1 DESCRIPTION
BN_mod_exp_mont() computes I<a> to the I<p>-th power modulo I<m> (C<rr=a^p % m>)
using Montgomery multiplication. I<in_mont> is a Montgomery context and can be
NULL. In the case I<in_mont> is NULL, it will be initialized within the
function, so you can save time on initialization if you provide it in advance.
BN_mod_exp_mont_consttime() computes I<a> to the I<p>-th power modulo I<m>
(C<rr=a^p % m>) using Montgomery multiplication. It is a variant of
L<BN_mod_exp_mont(3)> that uses fixed windows and the special precomputation
memory layout to limit data-dependency to a minimum to protect secret exponents.
It is called automatically when L<BN_mod_exp_mont(3)> is called with parameters
I<a>, I<p>, I<m>, any of which have B<BN_FLG_CONSTTIME> flag.
BN_mod_exp_mont_consttime_x2() computes two independent exponentiations I<a1> to
the I<p1>-th power modulo I<m1> (C<rr1=a1^p1 % m1>) and I<a2> to the I<p2>-th
power modulo I<m2> (C<rr2=a2^p2 % m2>) using Montgomery multiplication. For some
fixed and equal modulus sizes I<m1> and I<m2> it uses optimizations that allow
to speedup two exponentiations. In all other cases the function reduces to two
calls of L<BN_mod_exp_mont_consttime(3)>.
=head1 RETURN VALUES
For all functions 1 is returned for success, 0 on error.
The error codes can be obtained by L<ERR_get_error(3)>.
=head1 SEE ALSO
L<ERR_get_error(3)>, L<BN_mod_exp_mont(3)>
=head1 COPYRIGHT
Copyright 2000-2020 The OpenSSL Project Authors. All Rights Reserved.
Licensed under the Apache License 2.0 (the "License"). You may not use
this file except in compliance with the License. You can obtain a copy
in the file LICENSE in the source distribution or at
L<https://www.openssl.org/source/license.html>.
=cut

View File

@ -312,6 +312,11 @@ int BN_mod_exp2_mont(BIGNUM *r, const BIGNUM *a1, const BIGNUM *p1,
BN_CTX *ctx, BN_MONT_CTX *m_ctx);
int BN_mod_exp_simple(BIGNUM *r, const BIGNUM *a, const BIGNUM *p,
const BIGNUM *m, BN_CTX *ctx);
int BN_mod_exp_mont_consttime_x2(BIGNUM *rr1, const BIGNUM *a1, const BIGNUM *p1,
const BIGNUM *m1, BN_MONT_CTX *in_mont1,
BIGNUM *rr2, const BIGNUM *a2, const BIGNUM *p2,
const BIGNUM *m2, BN_MONT_CTX *in_mont2,
BN_CTX *ctx);
int BN_mask_bits(BIGNUM *a, int n);
# ifndef OPENSSL_NO_STDIO

View File

@ -1,5 +1,5 @@
/*
* Copyright 1995-2017 The OpenSSL Project Authors. All Rights Reserved.
* Copyright 1995-2020 The OpenSSL Project Authors. All Rights Reserved.
*
* Licensed under the Apache License 2.0 (the "License"). You may not use
* this file except in compliance with the License. You can obtain a copy
@ -198,9 +198,102 @@ static int test_mod_exp(int round)
return ret;
}
static int test_mod_exp_x2(int idx)
{
BN_CTX *ctx;
int ret = 0;
BIGNUM *r_mont_const_x2_1 = NULL;
BIGNUM *r_mont_const_x2_2 = NULL;
BIGNUM *r_simple1 = NULL;
BIGNUM *r_simple2 = NULL;
BIGNUM *a1 = NULL;
BIGNUM *b1 = NULL;
BIGNUM *m1 = NULL;
BIGNUM *a2 = NULL;
BIGNUM *b2 = NULL;
BIGNUM *m2 = NULL;
int factor_size = 0;
/*
* Currently only 1024-bit factor size is supported.
*/
if (idx <= 100)
factor_size = 1024;
if (!TEST_ptr(ctx = BN_CTX_new()))
goto err;
if (!TEST_ptr(r_mont_const_x2_1 = BN_new())
|| !TEST_ptr(r_mont_const_x2_2 = BN_new())
|| !TEST_ptr(r_simple1 = BN_new())
|| !TEST_ptr(r_simple2 = BN_new())
|| !TEST_ptr(a1 = BN_new())
|| !TEST_ptr(b1 = BN_new())
|| !TEST_ptr(m1 = BN_new())
|| !TEST_ptr(a2 = BN_new())
|| !TEST_ptr(b2 = BN_new())
|| !TEST_ptr(m2 = BN_new()))
goto err;
BN_rand(a1, factor_size, BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ANY);
BN_rand(b1, factor_size, BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ANY);
BN_rand(m1, factor_size, BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ODD);
BN_rand(a2, factor_size, BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ANY);
BN_rand(b2, factor_size, BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ANY);
BN_rand(m2, factor_size, BN_RAND_TOP_ONE, BN_RAND_BOTTOM_ODD);
if (!TEST_true(BN_mod(a1, a1, m1, ctx))
|| !TEST_true(BN_mod(b1, b1, m1, ctx))
|| !TEST_true(BN_mod(a2, a2, m2, ctx))
|| !TEST_true(BN_mod(b2, b2, m2, ctx))
|| !TEST_true(BN_mod_exp_simple(r_simple1, a1, b1, m1, ctx))
|| !TEST_true(BN_mod_exp_simple(r_simple2, a2, b2, m2, ctx))
|| !TEST_true(BN_mod_exp_mont_consttime_x2(r_mont_const_x2_1, a1, b1, m1, NULL,
r_mont_const_x2_2, a2, b2, m2, NULL,
ctx)))
goto err;
if (!TEST_BN_eq(r_simple1, r_mont_const_x2_1)
|| !TEST_BN_eq(r_simple2, r_mont_const_x2_2)) {
if (BN_cmp(r_simple1, r_mont_const_x2_1) != 0)
TEST_info("simple and mont const time x2 (#1) results differ");
if (BN_cmp(r_simple2, r_mont_const_x2_2) != 0)
TEST_info("simple and mont const time x2 (#2) results differ");
BN_print_var(a1);
BN_print_var(b1);
BN_print_var(m1);
BN_print_var(a2);
BN_print_var(b2);
BN_print_var(m2);
BN_print_var(r_simple1);
BN_print_var(r_simple2);
BN_print_var(r_mont_const_x2_1);
BN_print_var(r_mont_const_x2_2);
goto err;
}
ret = 1;
err:
BN_free(r_mont_const_x2_1);
BN_free(r_mont_const_x2_2);
BN_free(r_simple1);
BN_free(r_simple2);
BN_free(a1);
BN_free(b1);
BN_free(m1);
BN_free(a2);
BN_free(b2);
BN_free(m2);
BN_CTX_free(ctx);
return ret;
}
int setup_tests(void)
{
ADD_TEST(test_mod_exp_zero);
ADD_ALL_TESTS(test_mod_exp, 200);
ADD_ALL_TESTS(test_mod_exp_x2, 100);
return 1;
}

View File

@ -5313,6 +5313,7 @@ EVP_RAND_CTX_gettable_params ? 3_0_0 EXIST::FUNCTION:
EVP_RAND_CTX_settable_params ? 3_0_0 EXIST::FUNCTION:
RAND_set_DRBG_type ? 3_0_0 EXIST::FUNCTION:
RAND_set_seed_source_type ? 3_0_0 EXIST::FUNCTION:
BN_mod_exp_mont_consttime_x2 ? 3_0_0 EXIST::FUNCTION:
BIO_f_readbuffer ? 3_0_0 EXIST::FUNCTION:
EVP_DigestInit_ex2 ? 3_0_0 EXIST::FUNCTION:
EVP_EncryptInit_ex2 ? 3_0_0 EXIST::FUNCTION:

View File

@ -261,8 +261,6 @@ BN_is_negative(3)
BN_kronecker(3)
BN_mod_add_quick(3)
BN_mod_exp2_mont(3)
BN_mod_exp_mont(3)
BN_mod_exp_mont_consttime(3)
BN_mod_exp_mont_word(3)
BN_mod_exp_recp(3)
BN_mod_exp_simple(3)