openssl
/
openssl
mirror of https://github.com/openssl/openssl
1
0
Fork 0
Code Issues Releases Wiki Activity

Dual 1536/2048-bit exponentiation optimization for Intel IceLake CPU 

It uses AVX512_IFMA + AVX512_VL (with 256-bit wide registers) ISA to
keep lower power license.

Reviewed-by: Matt Caswell <matt@openssl.org>
Reviewed-by: Paul Dale <pauli@openssl.org>
(Merged from https://github.com/openssl/openssl/pull/14908)
This commit is contained in:
Andrey Matyukov 2020-12-08 22:53:39 +03:00 committed by Pauli
parent e67edf60f2
commit f87b4c4ea6
8 changed files with 2229 additions and 329 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
CHANGES.md
View File

@ -65,6 +65,11 @@ OpenSSL 3.1
*Dmitry Belyavskiy*
* Parallel dual-prime 1536/2048-bit modular exponentiation for
AVX512_IFMA capable processors.
*Sergey Kirillov, Andrey Matyukov (Intel Corp)*
OpenSSL 3.0
-----------

308
crypto/bn/asm/rsaz-avx512.pl → crypto/bn/asm/rsaz-2k-avx512.pl
View File

@ -7,7 +7,8 @@
# https://www.openssl.org/source/license.html
#
#
# Originally written by Ilya Albrekht, Sergey Kirillov and Andrey Matyukov
# Originally written by Sergey Kirillov and Andrey Matyukov.
# Special thanks to Ilya Albrekht for his valuable hints.
# Intel Corporation
#
# December 2020
@ -77,26 +78,29 @@ ___
###############################################################################
# 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.
# AMM is defined as presented in the paper [1].
#
# 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.
# NB: the AMM implementation does not perform "conditional" subtraction step
# specified in the original algorithm as according to the Lemma 1 from the paper
# [2], the result will be always < 2*m and can be used as a direct input to
# the next AMM iteration. This post-condition is true, provided the correct
# parameter |s| (notion of the Lemma 1 from [2]) is choosen, i.e. s >= n + 2 * k,
# which matches our case: 1040 > 1024 + 2 * 1.
#
# void ossl_rsaz_amm52x20_x1_256(BN_ULONG *res,
# const BN_ULONG *a,
# const BN_ULONG *b,
# const BN_ULONG *m,
# BN_ULONG k0);
# [1] Gueron, S. Efficient software implementations of modular exponentiation.
# DOI: 10.1007/s13389-012-0031-5
# [2] Gueron, S. Enhanced Montgomery Multiplication.
# DOI: 10.1007/3-540-36400-5_5
#
# void ossl_rsaz_amm52x20_x1_ifma256(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")
@ -112,16 +116,13 @@ 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";
my $Bi = "%ymm1";
my $Yi = "%ymm2";
my ($R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0) = ("%ymm3",map("%ymm$_",(16..19)));
my ($R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1) = ("%ymm4",map("%ymm$_",(20..23)));
# 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));
my ($T0,$T0h,$T1,$T1h,$T2) = ("$zero", "$Bi", "$Yi", map("%ymm$_", (25..26)));
sub amm52x20_x1() {
# _data_offset - offset in the |a| or |m| arrays pointing to the beginning
@ -190,16 +191,16 @@ $code.=<<___;
___
}
# Normalization routine: handles carry bits in R0..R2 QWs and
# gets R0..R2 back to normalized 2^52 representation.
# Normalization routine: handles carry bits and gets bignum qwords 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
vpbroadcastq $_acc, $T0
vpblendd \$3, $T0, $_R0, $_R0
# Extract "carries" (12 high bits) from each QW of R0..R2
# Save them to LSB of QWs in T0..T2
@ -214,14 +215,14 @@ $code.=<<___;
valignq \$3, $T1, $T1h, $T1h
valignq \$3, $T0h, $T1, $T1
valignq \$3, $T0, $T0h, $T0h
valignq \$3, $zero, $T0, $T0
valignq \$3, .Lzeros(%rip), $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
vpandq .Lmask52x4(%rip), $_R0, $_R0
vpandq .Lmask52x4(%rip), $_R0h, $_R0h
vpandq .Lmask52x4(%rip), $_R1, $_R1
vpandq .Lmask52x4(%rip), $_R1h, $_R1h
vpandq .Lmask52x4(%rip), $_R2, $_R2
# Sum R0..R2 with corresponding adjusted carries
vpaddq $T0, $_R0, $_R0
@ -232,11 +233,11 @@ $code.=<<___;
# 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
vpcmpuq \$6, .Lmask52x4(%rip), $_R0, %k1 # OP=nle (i.e. gt)
vpcmpuq \$6, .Lmask52x4(%rip), $_R0h, %k2
vpcmpuq \$6, .Lmask52x4(%rip), $_R1, %k3
vpcmpuq \$6, .Lmask52x4(%rip), $_R1h, %k4
vpcmpuq \$6, .Lmask52x4(%rip), $_R2, %k5
kmovb %k1, %r14d # k1
kmovb %k2, %r13d # k1h
kmovb %k3, %r12d # k2
@ -244,11 +245,11 @@ $code.=<<___;
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
vpcmpuq \$0, .Lmask52x4(%rip), $_R0, %k1 # OP=eq
vpcmpuq \$0, .Lmask52x4(%rip), $_R0h, %k2
vpcmpuq \$0, .Lmask52x4(%rip), $_R1, %k3
vpcmpuq \$0, .Lmask52x4(%rip), $_R1h, %k4
vpcmpuq \$0, .Lmask52x4(%rip), $_R2, %k5
kmovb %k1, %r9d # k4
kmovb %k2, %r8d # k4h
kmovb %k3, %ebx # k5
@ -288,27 +289,27 @@ $code.=<<___;
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}
vpsubq .Lmask52x4(%rip), $_R0, ${_R0}{%k1}
vpsubq .Lmask52x4(%rip), $_R0h, ${_R0h}{%k2}
vpsubq .Lmask52x4(%rip), $_R1, ${_R1}{%k3}
vpsubq .Lmask52x4(%rip), $_R1h, ${_R1h}{%k4}
vpsubq .Lmask52x4(%rip), $_R2, ${_R2}{%k5}
vpandq $mask52x4, $_R0, $_R0
vpandq $mask52x4, $_R0h, $_R0h
vpandq $mask52x4, $_R1, $_R1
vpandq $mask52x4, $_R1h, $_R1h
vpandq $mask52x4, $_R2, $_R2
vpandq .Lmask52x4(%rip), $_R0, $_R0
vpandq .Lmask52x4(%rip), $_R0h, $_R0h
vpandq .Lmask52x4(%rip), $_R1, $_R1
vpandq .Lmask52x4(%rip), $_R1h, $_R1h
vpandq .Lmask52x4(%rip), $_R2, $_R2
___
}
$code.=<<___;
.text
.globl ossl_rsaz_amm52x20_x1_256
.type ossl_rsaz_amm52x20_x1_256,\@function,5
.globl ossl_rsaz_amm52x20_x1_ifma256
.type ossl_rsaz_amm52x20_x1_ifma256,\@function,5
.align 32
ossl_rsaz_amm52x20_x1_256:
ossl_rsaz_amm52x20_x1_ifma256:
.cfi_startproc
endbranch
push %rbx
@ -323,7 +324,7 @@ ossl_rsaz_amm52x20_x1_256:
.cfi_push %r14
push %r15
.cfi_push %r15
.Lrsaz_amm52x20_x1_256_body:
.Lossl_rsaz_amm52x20_x1_ifma256_body:
# Zeroing accumulators
vpxord $zero, $zero, $zero
@ -351,17 +352,15 @@ $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)
vmovdqu64 $R0_0, `0*32`($res)
vmovdqu64 $R0_0h, `1*32`($res)
vmovdqu64 $R1_0, `2*32`($res)
vmovdqu64 $R1_0h, `3*32`($res)
vmovdqu64 $R2_0, `4*32`($res)
vzeroupper
mov 0(%rsp),%r15
@ -378,10 +377,10 @@ $code.=<<___;
.cfi_restore %rbx
lea 48(%rsp),%rsp
.cfi_adjust_cfa_offset -48
.Lrsaz_amm52x20_x1_256_epilogue:
.Lossl_rsaz_amm52x20_x1_ifma256_epilogue:
ret
.cfi_endproc
.size ossl_rsaz_amm52x20_x1_256, .-ossl_rsaz_amm52x20_x1_256
.size ossl_rsaz_amm52x20_x1_ifma256, .-ossl_rsaz_amm52x20_x1_ifma256
___
$code.=<<___;
@ -397,25 +396,25 @@ ___
###############################################################################
# Dual Almost Montgomery Multiplication for 20-digit number in radix 2^52
#
# See description of ossl_rsaz_amm52x20_x1_256() above for details about Almost
# See description of ossl_rsaz_amm52x20_x1_ifma256() 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 ossl_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]);
# void ossl_rsaz_amm52x20_x2_ifma256(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 ossl_rsaz_amm52x20_x2_256
.type ossl_rsaz_amm52x20_x2_256,\@function,5
.globl ossl_rsaz_amm52x20_x2_ifma256
.type ossl_rsaz_amm52x20_x2_ifma256,\@function,5
.align 32
ossl_rsaz_amm52x20_x2_256:
ossl_rsaz_amm52x20_x2_ifma256:
.cfi_startproc
endbranch
push %rbx
@ -430,7 +429,7 @@ ossl_rsaz_amm52x20_x2_256:
.cfi_push %r14
push %r15
.cfi_push %r15
.Lrsaz_amm52x20_x2_256_body:
.Lossl_rsaz_amm52x20_x2_ifma256_body:
# Zeroing accumulators
vpxord $zero, $zero, $zero
@ -463,24 +462,22 @@ $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_0, `0*32`($res)
vmovdqu64 $R0_0h, `1*32`($res)
vmovdqu64 $R1_0, `2*32`($res)
vmovdqu64 $R1_0h, `3*32`($res)
vmovdqu64 $R2_0, `4*32`($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)
vmovdqu64 $R0_1, `5*32`($res)
vmovdqu64 $R0_1h, `6*32`($res)
vmovdqu64 $R1_1, `7*32`($res)
vmovdqu64 $R1_1h, `8*32`($res)
vmovdqu64 $R2_1, `9*32`($res)
vzeroupper
mov 0(%rsp),%r15
@ -497,10 +494,10 @@ $code.=<<___;
.cfi_restore %rbx
lea 48(%rsp),%rsp
.cfi_adjust_cfa_offset -48
.Lrsaz_amm52x20_x2_256_epilogue:
.Lossl_rsaz_amm52x20_x2_ifma256_epilogue:
ret
.cfi_endproc
.size ossl_rsaz_amm52x20_x2_256, .-ossl_rsaz_amm52x20_x2_256
.size ossl_rsaz_amm52x20_x2_ifma256, .-ossl_rsaz_amm52x20_x2_ifma256
___
}
@ -508,77 +505,76 @@ ___
# 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.
# The input |red_table| contains precomputations for two independent base values.
# |red_table_idx1| and |red_table_idx2| are corresponding power indexes.
#
# Extracted value (output) is 20 digit number in 2^52 radix.
# Extracted value (output) is 2 20 digit numbers in 2^52 radix.
#
# void ossl_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
# int red_table_idx1, int red_table_idx2);
#
# EXP_WIN_SIZE = 5
###############################################################################
{
# input parameters
my ($out,$red_tbl,$red_tbl_idx,$tbl_idx) = @_6_args_universal_ABI;
my ($out,$red_tbl,$red_tbl_idx1,$red_tbl_idx2)=$win64 ? ("%rcx","%rdx","%r8", "%r9") : # Win64 order
("%rdi","%rsi","%rdx","%rcx"); # Unix order
my ($t0,$t1,$t2,$t3,$t4) = map("%ymm$_", (0..4));
my $t4xmm = $t4;
$t4xmm =~ s/%y/%x/;
my ($tmp0,$tmp1,$tmp2,$tmp3,$tmp4) = map("%ymm$_", (16..20));
my ($cur_idx,$idx,$ones) = map("%ymm$_", (21..23));
my ($t0,$t1,$t2,$t3,$t4,$t5) = map("%ymm$_", (0..5));
my ($t6,$t7,$t8,$t9) = map("%ymm$_", (16..19));
my ($tmp,$cur_idx,$idx1,$idx2,$ones) = map("%ymm$_", (20..24));
my @t = ($t0,$t1,$t2,$t3,$t4,$t5,$t6,$t7,$t8,$t9);
my $t0xmm = $t0;
$t0xmm =~ s/%y/%x/;
$code.=<<___;
.text
.align 32
.globl ossl_extract_multiplier_2x20_win5
.type ossl_extract_multiplier_2x20_win5,\@function,4
.type ossl_extract_multiplier_2x20_win5,\@abi-omnipotent
ossl_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
vpbroadcastq $red_tbl_idx1, $idx1
vpbroadcastq $red_tbl_idx2, $idx2
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
# zeroing t0..n, cur_idx
vpxor $t0xmm, $t0xmm, $t0xmm
vmovdqa64 $t0, $cur_idx
___
foreach (1..9) {
$code.="vmovdqa64 $t0, $t[$_] \n";
}
$code.=<<___;
.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}
vpcmpq \$0, $cur_idx, $idx1, %k1 # mask of (idx1 == cur_idx)
vpcmpq \$0, $cur_idx, $idx2, %k2 # mask of (idx2 == cur_idx)
___
foreach (0..9) {
my $mask = $_<5?"%k1":"%k2";
$code.=<<___;
vmovdqu64 `${_}*32`($red_tbl), $tmp # load data from red_tbl
vpblendmq $tmp, $t[$_], ${t[$_]}{$mask} # extract data when mask is not zero
___
}
$code.=<<___;
vpaddq $ones, $cur_idx, $cur_idx # increment cur_idx
addq \$`2*20*8`, $red_tbl
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)
___
# store t0..n
foreach (0..9) {
$code.="vmovdqu64 $t[$_], `${_}*32`($out) \n";
}
$code.=<<___;
ret
.cfi_endproc
.size ossl_extract_multiplier_2x20_win5, .-ossl_extract_multiplier_2x20_win5
@ -588,6 +584,8 @@ $code.=<<___;
.align 32
.Lones:
.quad 1,1,1,1
.Lzeros:
.quad 0,0,0,0
___
}
@ -597,7 +595,7 @@ $frame="%rdx";
$context="%r8";
$disp="%r9";
$code.=<<___
$code.=<<___;
.extern __imp_RtlVirtualUnwind
.type rsaz_def_handler,\@abi-omnipotent
.align 16
@ -688,32 +686,24 @@ rsaz_def_handler:
.section .pdata
.align 4
.rva .LSEH_begin_ossl_rsaz_amm52x20_x1_256
.rva .LSEH_end_ossl_rsaz_amm52x20_x1_256
.rva .LSEH_info_ossl_rsaz_amm52x20_x1_256
.rva .LSEH_begin_ossl_rsaz_amm52x20_x1_ifma256
.rva .LSEH_end_ossl_rsaz_amm52x20_x1_ifma256
.rva .LSEH_info_ossl_rsaz_amm52x20_x1_ifma256
.rva .LSEH_begin_ossl_rsaz_amm52x20_x2_256
.rva .LSEH_end_ossl_rsaz_amm52x20_x2_256
.rva .LSEH_info_ossl_rsaz_amm52x20_x2_256
.rva .LSEH_begin_ossl_extract_multiplier_2x20_win5
.rva .LSEH_end_ossl_extract_multiplier_2x20_win5
.rva .LSEH_info_ossl_extract_multiplier_2x20_win5
.rva .LSEH_begin_ossl_rsaz_amm52x20_x2_ifma256
.rva .LSEH_end_ossl_rsaz_amm52x20_x2_ifma256
.rva .LSEH_info_ossl_rsaz_amm52x20_x2_ifma256
.section .xdata
.align 8
.LSEH_info_ossl_rsaz_amm52x20_x1_256:
.LSEH_info_ossl_rsaz_amm52x20_x1_ifma256:
.byte 9,0,0,0
.rva rsaz_def_handler
.rva .Lrsaz_amm52x20_x1_256_body,.Lrsaz_amm52x20_x1_256_epilogue
.LSEH_info_ossl_rsaz_amm52x20_x2_256:
.rva .Lossl_rsaz_amm52x20_x1_ifma256_body,.Lossl_rsaz_amm52x20_x1_ifma256_epilogue
.LSEH_info_ossl_rsaz_amm52x20_x2_ifma256:
.byte 9,0,0,0
.rva rsaz_def_handler
.rva .Lrsaz_amm52x20_x2_256_body,.Lrsaz_amm52x20_x2_256_epilogue
.LSEH_info_ossl_extract_multiplier_2x20_win5:
.byte 9,0,0,0
.rva rsaz_def_handler
.rva .LSEH_begin_ossl_extract_multiplier_2x20_win5,.LSEH_begin_ossl_extract_multiplier_2x20_win5
.rva .Lossl_rsaz_amm52x20_x2_ifma256_body,.Lossl_rsaz_amm52x20_x2_ifma256_epilogue
___
}
}}} else {{{ # fallback for old assembler
@ -727,16 +717,16 @@ ossl_rsaz_avx512ifma_eligible:
ret
.size ossl_rsaz_avx512ifma_eligible, .-ossl_rsaz_avx512ifma_eligible
.globl ossl_rsaz_amm52x20_x1_256
.globl ossl_rsaz_amm52x20_x2_256
.globl ossl_rsaz_amm52x20_x1_ifma256
.globl ossl_rsaz_amm52x20_x2_ifma256
.globl ossl_extract_multiplier_2x20_win5
.type ossl_rsaz_amm52x20_x1_256,\@abi-omnipotent
ossl_rsaz_amm52x20_x1_256:
ossl_rsaz_amm52x20_x2_256:
.type ossl_rsaz_amm52x20_x1_ifma256,\@abi-omnipotent
ossl_rsaz_amm52x20_x1_ifma256:
ossl_rsaz_amm52x20_x2_ifma256:
ossl_extract_multiplier_2x20_win5:
.byte 0x0f,0x0b # ud2
ret
.size ossl_rsaz_amm52x20_x1_256, .-ossl_rsaz_amm52x20_x1_256
.size ossl_rsaz_amm52x20_x1_ifma256, .-ossl_rsaz_amm52x20_x1_ifma256
___
}}}

874
crypto/bn/asm/rsaz-3k-avx512.pl Normal file
View File

@ -0,0 +1,874 @@
# Copyright 2021 The OpenSSL Project Authors. All Rights Reserved.
# Copyright (c) 2021, 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 Sergey Kirillov and Andrey Matyukov
# Intel Corporation
#
# March 2021
#
# Initial release.
#
# Implementation utilizes 256-bit (ymm) registers to avoid frequency scaling issues.
#
# IceLake-Client @ 1.3GHz
# |---------+-----------------------+---------------+-------------|
# | | OpenSSL 3.0.0-alpha15 | this | Unit |
# |---------+-----------------------+---------------+-------------|
# | rsa3072 | 6 397 637 | 2 866 593 | cycles/sign |
# | | 203.2 | 453.5 / +123% | 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>=7.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");
###############################################################################
# Almost Montgomery Multiplication (AMM) for 30-digit number in radix 2^52.
#
# AMM is defined as presented in the paper [1].
#
# The input and output are presented in 2^52 radix domain, i.e.
# |res|, |a|, |b|, |m| are arrays of 32 64-bit qwords with 12 high bits zeroed
#
# NOTE: the function uses zero-padded data - 2 high QWs is a padding.
#
# |k0| is a Montgomery coefficient, which is here k0 = -1/m mod 2^64
#
# NB: the AMM implementation does not perform "conditional" subtraction step
# specified in the original algorithm as according to the Lemma 1 from the paper
# [2], the result will be always < 2*m and can be used as a direct input to
# the next AMM iteration. This post-condition is true, provided the correct
# parameter |s| (notion of the Lemma 1 from [2]) is choosen, i.e. s >= n + 2 * k,
# which matches our case: 1560 > 1536 + 2 * 1.
#
# [1] Gueron, S. Efficient software implementations of modular exponentiation.
# DOI: 10.1007/s13389-012-0031-5
# [2] Gueron, S. Enhanced Montgomery Multiplication.
# DOI: 10.1007/3-540-36400-5_5
#
# void ossl_rsaz_amm52x30_x1_ifma256(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 $Bi = "%ymm1";
my $Yi = "%ymm2";
my ($R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h) = map("%ymm$_",(3..10));
my ($R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1,$R2_1h,$R3_1,$R3_1h) = map("%ymm$_",(11..18));
# Registers mapping for normalization
my ($T0,$T0h,$T1,$T1h,$T2,$T2h,$T3,$T3h) = ("$zero", "$Bi", "$Yi", map("%ymm$_", (19..23)));
sub amm52x30_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,$_R2h,$_R3,$_R3h,$_k0) = @_;
my $_R0_xmm = $_R0;
$_R0_xmm =~ s/%y/%x/;
$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*2+32`($a), $Bi, $_R2h
vpmadd52luq `$_data_offset+64*3`($a), $Bi, $_R3
vpmadd52luq `$_data_offset+64*3+32`($a), $Bi, $_R3h
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
vpmadd52luq `$_data_offset+64*2+32`($m), $Yi, $_R2h
vpmadd52luq `$_data_offset+64*3`($m), $Yi, $_R3
vpmadd52luq `$_data_offset+64*3+32`($m), $Yi, $_R3h
# 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, $_R2h, $_R2
valignq \$1, $_R2h, $_R3, $_R2h
valignq \$1, $_R3, $_R3h, $_R3
valignq \$1, $_R3h, $zero, $_R3h
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*2+32`($a), $Bi, $_R2h
vpmadd52huq `$_data_offset+64*3`($a), $Bi, $_R3
vpmadd52huq `$_data_offset+64*3+32`($a), $Bi, $_R3h
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
vpmadd52huq `$_data_offset+64*2+32`($m), $Yi, $_R2h
vpmadd52huq `$_data_offset+64*3`($m), $Yi, $_R3
vpmadd52huq `$_data_offset+64*3+32`($m), $Yi, $_R3h
___
}
# Normalization routine: handles carry bits and gets bignum qwords to normalized
# 2^52 representation.
#
# Uses %r8-14,%e[abcd]x
sub amm52x30_x1_norm {
my ($_acc,$_R0,$_R0h,$_R1,$_R1h,$_R2,$_R2h,$_R3,$_R3h) = @_;
$code.=<<___;
# Put accumulator to low qword in R0
vpbroadcastq $_acc, $T0
vpblendd \$3, $T0, $_R0, $_R0
# Extract "carries" (12 high bits) from each QW of the bignum
# Save them to LSB of QWs in T0..Tn
vpsrlq \$52, $_R0, $T0
vpsrlq \$52, $_R0h, $T0h
vpsrlq \$52, $_R1, $T1
vpsrlq \$52, $_R1h, $T1h
vpsrlq \$52, $_R2, $T2
vpsrlq \$52, $_R2h, $T2h
vpsrlq \$52, $_R3, $T3
vpsrlq \$52, $_R3h, $T3h
# "Shift left" T0..Tn by 1 QW
valignq \$3, $T3, $T3h, $T3h
valignq \$3, $T2h, $T3, $T3
valignq \$3, $T2, $T2h, $T2h
valignq \$3, $T1h, $T2, $T2
valignq \$3, $T1, $T1h, $T1h
valignq \$3, $T0h, $T1, $T1
valignq \$3, $T0, $T0h, $T0h
valignq \$3, .Lzeros(%rip), $T0, $T0
# Drop "carries" from R0..Rn QWs
vpandq .Lmask52x4(%rip), $_R0, $_R0
vpandq .Lmask52x4(%rip), $_R0h, $_R0h
vpandq .Lmask52x4(%rip), $_R1, $_R1
vpandq .Lmask52x4(%rip), $_R1h, $_R1h
vpandq .Lmask52x4(%rip), $_R2, $_R2
vpandq .Lmask52x4(%rip), $_R2h, $_R2h
vpandq .Lmask52x4(%rip), $_R3, $_R3
vpandq .Lmask52x4(%rip), $_R3h, $_R3h
# Sum R0..Rn with corresponding adjusted carries
vpaddq $T0, $_R0, $_R0
vpaddq $T0h, $_R0h, $_R0h
vpaddq $T1, $_R1, $_R1
vpaddq $T1h, $_R1h, $_R1h
vpaddq $T2, $_R2, $_R2
vpaddq $T2h, $_R2h, $_R2h
vpaddq $T3, $_R3, $_R3
vpaddq $T3h, $_R3h, $_R3h
# Now handle carry bits from this addition
# Get mask of QWs whose 52-bit parts overflow
vpcmpuq \$6,.Lmask52x4(%rip),${_R0},%k1 # OP=nle (i.e. gt)
vpcmpuq \$6,.Lmask52x4(%rip),${_R0h},%k2
kmovb %k1,%r14d
kmovb %k2,%r13d
shl \$4,%r13b
or %r13b,%r14b
vpcmpuq \$6,.Lmask52x4(%rip),${_R1},%k1
vpcmpuq \$6,.Lmask52x4(%rip),${_R1h},%k2
kmovb %k1,%r13d
kmovb %k2,%r12d
shl \$4,%r12b
or %r12b,%r13b
vpcmpuq \$6,.Lmask52x4(%rip),${_R2},%k1
vpcmpuq \$6,.Lmask52x4(%rip),${_R2h},%k2
kmovb %k1,%r12d
kmovb %k2,%r11d
shl \$4,%r11b
or %r11b,%r12b
vpcmpuq \$6,.Lmask52x4(%rip),${_R3},%k1
vpcmpuq \$6,.Lmask52x4(%rip),${_R3h},%k2
kmovb %k1,%r11d
kmovb %k2,%r10d
shl \$4,%r10b
or %r10b,%r11b
addb %r14b,%r14b
adcb %r13b,%r13b
adcb %r12b,%r12b
adcb %r11b,%r11b
# Get mask of QWs whose 52-bit parts saturated
vpcmpuq \$0,.Lmask52x4(%rip),${_R0},%k1 # OP=eq
vpcmpuq \$0,.Lmask52x4(%rip),${_R0h},%k2
kmovb %k1,%r9d
kmovb %k2,%r8d
shl \$4,%r8b
or %r8b,%r9b
vpcmpuq \$0,.Lmask52x4(%rip),${_R1},%k1
vpcmpuq \$0,.Lmask52x4(%rip),${_R1h},%k2
kmovb %k1,%r8d
kmovb %k2,%edx
shl \$4,%dl
or %dl,%r8b
vpcmpuq \$0,.Lmask52x4(%rip),${_R2},%k1
vpcmpuq \$0,.Lmask52x4(%rip),${_R2h},%k2
kmovb %k1,%edx
kmovb %k2,%ecx
shl \$4,%cl
or %cl,%dl
vpcmpuq \$0,.Lmask52x4(%rip),${_R3},%k1
vpcmpuq \$0,.Lmask52x4(%rip),${_R3h},%k2
kmovb %k1,%ecx
kmovb %k2,%ebx
shl \$4,%bl
or %bl,%cl
addb %r9b,%r14b
adcb %r8b,%r13b
adcb %dl,%r12b
adcb %cl,%r11b
xor %r9b,%r14b
xor %r8b,%r13b
xor %dl,%r12b
xor %cl,%r11b
kmovb %r14d,%k1
shr \$4,%r14b
kmovb %r14d,%k2
kmovb %r13d,%k3
shr \$4,%r13b
kmovb %r13d,%k4
kmovb %r12d,%k5
shr \$4,%r12b
kmovb %r12d,%k6
kmovb %r11d,%k7
vpsubq .Lmask52x4(%rip), $_R0, ${_R0}{%k1}
vpsubq .Lmask52x4(%rip), $_R0h, ${_R0h}{%k2}
vpsubq .Lmask52x4(%rip), $_R1, ${_R1}{%k3}
vpsubq .Lmask52x4(%rip), $_R1h, ${_R1h}{%k4}
vpsubq .Lmask52x4(%rip), $_R2, ${_R2}{%k5}
vpsubq .Lmask52x4(%rip), $_R2h, ${_R2h}{%k6}
vpsubq .Lmask52x4(%rip), $_R3, ${_R3}{%k7}
vpandq .Lmask52x4(%rip), $_R0, $_R0
vpandq .Lmask52x4(%rip), $_R0h, $_R0h
vpandq .Lmask52x4(%rip), $_R1, $_R1
vpandq .Lmask52x4(%rip), $_R1h, $_R1h
vpandq .Lmask52x4(%rip), $_R2, $_R2
vpandq .Lmask52x4(%rip), $_R2h, $_R2h
vpandq .Lmask52x4(%rip), $_R3, $_R3
shr \$4,%r11b
kmovb %r11d,%k1
vpsubq .Lmask52x4(%rip), $_R3h, ${_R3h}{%k1}
vpandq .Lmask52x4(%rip), $_R3h, $_R3h
___
}
$code.=<<___;
.text
.globl ossl_rsaz_amm52x30_x1_ifma256
.type ossl_rsaz_amm52x30_x1_ifma256,\@function,5
.align 32
ossl_rsaz_amm52x30_x1_ifma256:
.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
___
$code.=<<___ if ($win64);
lea -168(%rsp),%rsp # 16*10 + (8 bytes to get correct 16-byte SIMD alignment)
vmovdqa64 %xmm6, `0*16`(%rsp) # save non-volatile registers
vmovdqa64 %xmm7, `1*16`(%rsp)
vmovdqa64 %xmm8, `2*16`(%rsp)
vmovdqa64 %xmm9, `3*16`(%rsp)
vmovdqa64 %xmm10,`4*16`(%rsp)
vmovdqa64 %xmm11,`5*16`(%rsp)
vmovdqa64 %xmm12,`6*16`(%rsp)
vmovdqa64 %xmm13,`7*16`(%rsp)
vmovdqa64 %xmm14,`8*16`(%rsp)
vmovdqa64 %xmm15,`9*16`(%rsp)
.Lossl_rsaz_amm52x30_x1_ifma256_body:
___
$code.=<<___;
# 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, $R2_0h
vmovdqa64 $zero, $R3_0
vmovdqa64 $zero, $R3_0h
xorl $acc0_0_low, $acc0_0_low
movq $b, $b_ptr # backup address of b
movq \$0xfffffffffffff, $mask52 # 52-bit mask
# Loop over 30 digits unrolled by 4
mov \$7, $iter
.align 32
.Lloop7:
___
foreach my $idx (0..3) {
&amm52x30_x1(0,8*$idx,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h,$k0);
}
$code.=<<___;
lea `4*8`($b_ptr), $b_ptr
dec $iter
jne .Lloop7
___
&amm52x30_x1(0,8*0,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h,$k0);
&amm52x30_x1(0,8*1,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h,$k0);
&amm52x30_x1_norm($acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h);
$code.=<<___;
vmovdqu64 $R0_0, `0*32`($res)
vmovdqu64 $R0_0h, `1*32`($res)
vmovdqu64 $R1_0, `2*32`($res)
vmovdqu64 $R1_0h, `3*32`($res)
vmovdqu64 $R2_0, `4*32`($res)
vmovdqu64 $R2_0h, `5*32`($res)
vmovdqu64 $R3_0, `6*32`($res)
vmovdqu64 $R3_0h, `7*32`($res)
vzeroupper
lea (%rsp),%rax
.cfi_def_cfa_register %rax
___
$code.=<<___ if ($win64);
vmovdqa64 `0*16`(%rax),%xmm6
vmovdqa64 `1*16`(%rax),%xmm7
vmovdqa64 `2*16`(%rax),%xmm8
vmovdqa64 `3*16`(%rax),%xmm9
vmovdqa64 `4*16`(%rax),%xmm10
vmovdqa64 `5*16`(%rax),%xmm11
vmovdqa64 `6*16`(%rax),%xmm12
vmovdqa64 `7*16`(%rax),%xmm13
vmovdqa64 `8*16`(%rax),%xmm14
vmovdqa64 `9*16`(%rax),%xmm15
lea 168(%rsp),%rax
___
$code.=<<___;
mov 0(%rax),%r15
.cfi_restore %r15
mov 8(%rax),%r14
.cfi_restore %r14
mov 16(%rax),%r13
.cfi_restore %r13
mov 24(%rax),%r12
.cfi_restore %r12
mov 32(%rax),%rbp
.cfi_restore %rbp
mov 40(%rax),%rbx
.cfi_restore %rbx
lea 48(%rax),%rsp # restore rsp
.cfi_def_cfa %rsp,8
.Lossl_rsaz_amm52x30_x1_ifma256_epilogue:
ret
.cfi_endproc
.size ossl_rsaz_amm52x30_x1_ifma256, .-ossl_rsaz_amm52x30_x1_ifma256
___
$code.=<<___;
.data
.align 32
.Lmask52x4:
.quad 0xfffffffffffff
.quad 0xfffffffffffff
.quad 0xfffffffffffff
.quad 0xfffffffffffff
___
###############################################################################
# Dual Almost Montgomery Multiplication for 30-digit number in radix 2^52
#
# See description of ossl_rsaz_amm52x30_x1_ifma256() above for details about Almost
# Montgomery Multiplication algorithm and function input parameters description.
#
# This function does two AMMs for two independent inputs, hence dual.
#
# NOTE: the function uses zero-padded data - 2 high QWs is a padding.
#
# void ossl_rsaz_amm52x30_x2_ifma256(BN_ULONG out[2][32],
# const BN_ULONG a[2][32],
# const BN_ULONG b[2][32],
# const BN_ULONG m[2][32],
# const BN_ULONG k0[2]);
###############################################################################
$code.=<<___;
.text
.globl ossl_rsaz_amm52x30_x2_ifma256
.type ossl_rsaz_amm52x30_x2_ifma256,\@function,5
.align 32
ossl_rsaz_amm52x30_x2_ifma256:
.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
___
$code.=<<___ if ($win64);
lea -168(%rsp),%rsp
vmovdqa64 %xmm6, `0*16`(%rsp) # save non-volatile registers
vmovdqa64 %xmm7, `1*16`(%rsp)
vmovdqa64 %xmm8, `2*16`(%rsp)
vmovdqa64 %xmm9, `3*16`(%rsp)
vmovdqa64 %xmm10,`4*16`(%rsp)
vmovdqa64 %xmm11,`5*16`(%rsp)
vmovdqa64 %xmm12,`6*16`(%rsp)
vmovdqa64 %xmm13,`7*16`(%rsp)
vmovdqa64 %xmm14,`8*16`(%rsp)
vmovdqa64 %xmm15,`9*16`(%rsp)
.Lossl_rsaz_amm52x30_x2_ifma256_body:
___
$code.=<<___;
# 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, $R2_0h
vmovdqa64 $zero, $R3_0
vmovdqa64 $zero, $R3_0h
vmovdqa64 $zero, $R0_1
vmovdqa64 $zero, $R0_1h
vmovdqa64 $zero, $R1_1
vmovdqa64 $zero, $R1_1h
vmovdqa64 $zero, $R2_1
vmovdqa64 $zero, $R2_1h
vmovdqa64 $zero, $R3_1
vmovdqa64 $zero, $R3_1h
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 \$30, $iter
.align 32
.Lloop30:
___
&amm52x30_x1( 0, 0,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h,"($k0)");
# 32*8 = offset of the next dimension in two-dimension array
&amm52x30_x1(32*8,32*8,$acc0_1,$R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1,$R2_1h,$R3_1,$R3_1h,"8($k0)");
$code.=<<___;
lea 8($b_ptr), $b_ptr
dec $iter
jne .Lloop30
___
&amm52x30_x1_norm($acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h);
&amm52x30_x1_norm($acc0_1,$R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1,$R2_1h,$R3_1,$R3_1h);
$code.=<<___;
vmovdqu64 $R0_0, `0*32`($res)
vmovdqu64 $R0_0h, `1*32`($res)
vmovdqu64 $R1_0, `2*32`($res)
vmovdqu64 $R1_0h, `3*32`($res)
vmovdqu64 $R2_0, `4*32`($res)
vmovdqu64 $R2_0h, `5*32`($res)
vmovdqu64 $R3_0, `6*32`($res)
vmovdqu64 $R3_0h, `7*32`($res)
vmovdqu64 $R0_1, `8*32`($res)
vmovdqu64 $R0_1h, `9*32`($res)
vmovdqu64 $R1_1, `10*32`($res)
vmovdqu64 $R1_1h, `11*32`($res)
vmovdqu64 $R2_1, `12*32`($res)
vmovdqu64 $R2_1h, `13*32`($res)
vmovdqu64 $R3_1, `14*32`($res)
vmovdqu64 $R3_1h, `15*32`($res)
vzeroupper
lea (%rsp),%rax
.cfi_def_cfa_register %rax
___
$code.=<<___ if ($win64);
vmovdqa64 `0*16`(%rax),%xmm6
vmovdqa64 `1*16`(%rax),%xmm7
vmovdqa64 `2*16`(%rax),%xmm8
vmovdqa64 `3*16`(%rax),%xmm9
vmovdqa64 `4*16`(%rax),%xmm10
vmovdqa64 `5*16`(%rax),%xmm11
vmovdqa64 `6*16`(%rax),%xmm12
vmovdqa64 `7*16`(%rax),%xmm13
vmovdqa64 `8*16`(%rax),%xmm14
vmovdqa64 `9*16`(%rax),%xmm15
lea 168(%rsp),%rax
___
$code.=<<___;
mov 0(%rax),%r15
.cfi_restore %r15
mov 8(%rax),%r14
.cfi_restore %r14
mov 16(%rax),%r13
.cfi_restore %r13
mov 24(%rax),%r12
.cfi_restore %r12
mov 32(%rax),%rbp
.cfi_restore %rbp
mov 40(%rax),%rbx
.cfi_restore %rbx
lea 48(%rax),%rsp
.cfi_def_cfa %rsp,8
.Lossl_rsaz_amm52x30_x2_ifma256_epilogue:
ret
.cfi_endproc
.size ossl_rsaz_amm52x30_x2_ifma256, .-ossl_rsaz_amm52x30_x2_ifma256
___
}
###############################################################################
# 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.
# |red_table_idx1| and |red_table_idx2| are corresponding power indexes.
#
# Extracted value (output) is 2 (30 + 2) digits numbers in 2^52 radix.
# (2 high QW is zero padding)
#
# void ossl_extract_multiplier_2x30_win5(BN_ULONG *red_Y,
# const BN_ULONG red_table[1 << EXP_WIN_SIZE][2][32],
# int red_table_idx1, int red_table_idx2);
#
# EXP_WIN_SIZE = 5
###############################################################################
{
# input parameters
my ($out,$red_tbl,$red_tbl_idx1,$red_tbl_idx2)=$win64 ? ("%rcx","%rdx","%r8", "%r9") : # Win64 order
("%rdi","%rsi","%rdx","%rcx"); # Unix order
my ($t0,$t1,$t2,$t3,$t4,$t5) = map("%ymm$_", (0..5));
my ($t6,$t7,$t8,$t9,$t10,$t11,$t12,$t13,$t14,$t15) = map("%ymm$_", (16..25));
my ($tmp,$cur_idx,$idx1,$idx2,$ones) = map("%ymm$_", (26..30));
my @t = ($t0,$t1,$t2,$t3,$t4,$t5,$t6,$t7,$t8,$t9,$t10,$t11,$t12,$t13,$t14,$t15);
my $t0xmm = $t0;
$t0xmm =~ s/%y/%x/;
$code.=<<___;
.text
.align 32
.globl ossl_extract_multiplier_2x30_win5
.type ossl_extract_multiplier_2x30_win5,\@abi-omnipotent
ossl_extract_multiplier_2x30_win5:
.cfi_startproc
endbranch
vmovdqa64 .Lones(%rip), $ones # broadcast ones
vpbroadcastq $red_tbl_idx1, $idx1
vpbroadcastq $red_tbl_idx2, $idx2
leaq `(1<<5)*2*32*8`($red_tbl), %rax # holds end of the tbl
# zeroing t0..n, cur_idx
vpxor $t0xmm, $t0xmm, $t0xmm
vmovdqa64 $t0, $cur_idx
___
foreach (1..15) {
$code.="vmovdqa64 $t0, $t[$_] \n";
}
$code.=<<___;
.align 32
.Lloop:
vpcmpq \$0, $cur_idx, $idx1, %k1 # mask of (idx1 == cur_idx)
vpcmpq \$0, $cur_idx, $idx2, %k2 # mask of (idx2 == cur_idx)
___
foreach (0..15) {
my $mask = $_<8?"%k1":"%k2";
$code.=<<___;
vmovdqu64 `${_}*32`($red_tbl), $tmp # load data from red_tbl
vpblendmq $tmp, $t[$_], ${t[$_]}{$mask} # extract data when mask is not zero
___
}
$code.=<<___;
vpaddq $ones, $cur_idx, $cur_idx # increment cur_idx
addq \$`2*32*8`, $red_tbl
cmpq $red_tbl, %rax
jne .Lloop
___
# store t0..n
foreach (0..15) {
$code.="vmovdqu64 $t[$_], `${_}*32`($out) \n";
}
$code.=<<___;
ret
.cfi_endproc
.size ossl_extract_multiplier_2x30_win5, .-ossl_extract_multiplier_2x30_win5
___
$code.=<<___;
.data
.align 32
.Lones:
.quad 1,1,1,1
.Lzeros:
.quad 0,0,0,0
___
}
if ($win64) {
$rec="%rcx";
$frame="%rdx";
$context="%r8";
$disp="%r9";
$code.=<<___;
.extern __imp_RtlVirtualUnwind
.type rsaz_avx_handler,\@abi-omnipotent
.align 16
rsaz_avx_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 4(%r11),%r10d # HandlerData[1]
lea (%rsi,%r10),%r10 # epilogue label
cmp %r10,%rbx # context->Rip>=.Lepilogue
jae .Lcommon_seh_tail
mov 152($context),%rax # pull context->Rsp
lea (%rax),%rsi # %xmm save area
lea 512($context),%rdi # & context.Xmm6
mov \$20,%ecx # 10*sizeof(%xmm0)/sizeof(%rax)
.long 0xa548f3fc # cld; rep movsq
lea `48+168`(%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_avx_handler,.-rsaz_avx_handler
.section .pdata
.align 4
.rva .LSEH_begin_ossl_rsaz_amm52x30_x1_ifma256
.rva .LSEH_end_ossl_rsaz_amm52x30_x1_ifma256
.rva .LSEH_info_ossl_rsaz_amm52x30_x1_ifma256
.rva .LSEH_begin_ossl_rsaz_amm52x30_x2_ifma256
.rva .LSEH_end_ossl_rsaz_amm52x30_x2_ifma256
.rva .LSEH_info_ossl_rsaz_amm52x30_x2_ifma256
.section .xdata
.align 8
.LSEH_info_ossl_rsaz_amm52x30_x1_ifma256:
.byte 9,0,0,0
.rva rsaz_avx_handler
.rva .Lossl_rsaz_amm52x30_x1_ifma256_body,.Lossl_rsaz_amm52x30_x1_ifma256_epilogue
.LSEH_info_ossl_rsaz_amm52x30_x2_ifma256:
.byte 9,0,0,0
.rva rsaz_avx_handler
.rva .Lossl_rsaz_amm52x30_x2_ifma256_body,.Lossl_rsaz_amm52x30_x2_ifma256_epilogue
___
}
}}} else {{{ # fallback for old assembler
$code.=<<___;
.text
.globl ossl_rsaz_amm52x30_x1_ifma256
.globl ossl_rsaz_amm52x30_x2_ifma256
.globl ossl_extract_multiplier_2x30_win5
.type ossl_rsaz_amm52x30_x1_ifma256,\@abi-omnipotent
ossl_rsaz_amm52x30_x1_ifma256:
ossl_rsaz_amm52x30_x2_ifma256:
ossl_extract_multiplier_2x30_win5:
.byte 0x0f,0x0b # ud2
ret
.size ossl_rsaz_amm52x30_x1_ifma256, .-ossl_rsaz_amm52x30_x1_ifma256
___
}}}
$code =~ s/\`([^\`]*)\`/eval $1/gem;
print $code;
close STDOUT or die "error closing STDOUT: $!";

930
crypto/bn/asm/rsaz-4k-avx512.pl Normal file
View File

@ -0,0 +1,930 @@
# Copyright 2021 The OpenSSL Project Authors. All Rights Reserved.
# Copyright (c) 2021, 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 Sergey Kirillov and Andrey Matyukov
# Intel Corporation
#
# March 2021
#
# Initial release.
#
# Implementation utilizes 256-bit (ymm) registers to avoid frequency scaling issues.
#
# IceLake-Client @ 1.3GHz
# |---------+-----------------------+---------------+-------------|
# | | OpenSSL 3.0.0-alpha15 | this | Unit |
# |---------+-----------------------+---------------+-------------|
# | rsa4096 | 14 301 4300 | 5 813 953 | cycles/sign |
# | | 90.9 | 223.6 / +146% | 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>=7.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");
###############################################################################
# Almost Montgomery Multiplication (AMM) for 40-digit number in radix 2^52.
#
# AMM is defined as presented in the paper [1].
#
# The input and output are presented in 2^52 radix domain, i.e.
# |res|, |a|, |b|, |m| are arrays of 40 64-bit qwords with 12 high bits zeroed.
# |k0| is a Montgomery coefficient, which is here k0 = -1/m mod 2^64
#
# NB: the AMM implementation does not perform "conditional" subtraction step
# specified in the original algorithm as according to the Lemma 1 from the paper
# [2], the result will be always < 2*m and can be used as a direct input to
# the next AMM iteration. This post-condition is true, provided the correct
# parameter |s| (notion of the Lemma 1 from [2]) is choosen, i.e. s >= n + 2 * k,
# which matches our case: 2080 > 2048 + 2 * 1.
#
# [1] Gueron, S. Efficient software implementations of modular exponentiation.
# DOI: 10.1007/s13389-012-0031-5
# [2] Gueron, S. Enhanced Montgomery Multiplication.
# DOI: 10.1007/3-540-36400-5_5
#
# void ossl_rsaz_amm52x40_x1_ifma256(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 $Bi = "%ymm1";
my $Yi = "%ymm2";
my ($R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h,$R4_0,$R4_0h) = map("%ymm$_",(3..12));
my ($R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1,$R2_1h,$R3_1,$R3_1h,$R4_1,$R4_1h) = map("%ymm$_",(13..22));
# Registers mapping for normalization
my ($T0,$T0h,$T1,$T1h,$T2,$T2h,$T3,$T3h,$T4,$T4h) = ("$zero", "$Bi", "$Yi", map("%ymm$_", (23..29)));
sub amm52x40_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,$_R2h,$_R3,$_R3h,$_R4,$_R4h,$_k0) = @_;
my $_R0_xmm = $_R0;
$_R0_xmm =~ s/%y/%x/;
$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*2+32`($a), $Bi, $_R2h
vpmadd52luq `$_data_offset+64*3`($a), $Bi, $_R3
vpmadd52luq `$_data_offset+64*3+32`($a), $Bi, $_R3h
vpmadd52luq `$_data_offset+64*4`($a), $Bi, $_R4
vpmadd52luq `$_data_offset+64*4+32`($a), $Bi, $_R4h
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
vpmadd52luq `$_data_offset+64*2+32`($m), $Yi, $_R2h
vpmadd52luq `$_data_offset+64*3`($m), $Yi, $_R3
vpmadd52luq `$_data_offset+64*3+32`($m), $Yi, $_R3h
vpmadd52luq `$_data_offset+64*4`($m), $Yi, $_R4
vpmadd52luq `$_data_offset+64*4+32`($m), $Yi, $_R4h
# 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, $_R2h, $_R2
valignq \$1, $_R2h, $_R3, $_R2h
valignq \$1, $_R3, $_R3h, $_R3
valignq \$1, $_R3h, $_R4, $_R3h
valignq \$1, $_R4, $_R4h, $_R4
valignq \$1, $_R4h, $zero, $_R4h
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*2+32`($a), $Bi, $_R2h
vpmadd52huq `$_data_offset+64*3`($a), $Bi, $_R3
vpmadd52huq `$_data_offset+64*3+32`($a), $Bi, $_R3h
vpmadd52huq `$_data_offset+64*4`($a), $Bi, $_R4
vpmadd52huq `$_data_offset+64*4+32`($a), $Bi, $_R4h
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
vpmadd52huq `$_data_offset+64*2+32`($m), $Yi, $_R2h
vpmadd52huq `$_data_offset+64*3`($m), $Yi, $_R3
vpmadd52huq `$_data_offset+64*3+32`($m), $Yi, $_R3h
vpmadd52huq `$_data_offset+64*4`($m), $Yi, $_R4
vpmadd52huq `$_data_offset+64*4+32`($m), $Yi, $_R4h
___
}
# Normalization routine: handles carry bits and gets bignum qwords to normalized
# 2^52 representation.
#
# Uses %r8-14,%e[abcd]x
sub amm52x40_x1_norm {
my ($_acc,$_R0,$_R0h,$_R1,$_R1h,$_R2,$_R2h,$_R3,$_R3h,$_R4,$_R4h) = @_;
$code.=<<___;
# Put accumulator to low qword in R0
vpbroadcastq $_acc, $T0
vpblendd \$3, $T0, $_R0, $_R0
# Extract "carries" (12 high bits) from each QW of the bignum
# Save them to LSB of QWs in T0..Tn
vpsrlq \$52, $_R0, $T0
vpsrlq \$52, $_R0h, $T0h
vpsrlq \$52, $_R1, $T1
vpsrlq \$52, $_R1h, $T1h
vpsrlq \$52, $_R2, $T2
vpsrlq \$52, $_R2h, $T2h
vpsrlq \$52, $_R3, $T3
vpsrlq \$52, $_R3h, $T3h
vpsrlq \$52, $_R4, $T4
vpsrlq \$52, $_R4h, $T4h
# "Shift left" T0..Tn by 1 QW
valignq \$3, $T4, $T4h, $T4h
valignq \$3, $T3h, $T4, $T4
valignq \$3, $T3, $T3h, $T3h
valignq \$3, $T2h, $T3, $T3
valignq \$3, $T2, $T2h, $T2h
valignq \$3, $T1h, $T2, $T2
valignq \$3, $T1, $T1h, $T1h
valignq \$3, $T0h, $T1, $T1
valignq \$3, $T0, $T0h, $T0h
valignq \$3, .Lzeros(%rip), $T0, $T0
# Drop "carries" from R0..Rn QWs
vpandq .Lmask52x4(%rip), $_R0, $_R0
vpandq .Lmask52x4(%rip), $_R0h, $_R0h
vpandq .Lmask52x4(%rip), $_R1, $_R1
vpandq .Lmask52x4(%rip), $_R1h, $_R1h
vpandq .Lmask52x4(%rip), $_R2, $_R2
vpandq .Lmask52x4(%rip), $_R2h, $_R2h
vpandq .Lmask52x4(%rip), $_R3, $_R3
vpandq .Lmask52x4(%rip), $_R3h, $_R3h
vpandq .Lmask52x4(%rip), $_R4, $_R4
vpandq .Lmask52x4(%rip), $_R4h, $_R4h
# Sum R0..Rn with corresponding adjusted carries
vpaddq $T0, $_R0, $_R0
vpaddq $T0h, $_R0h, $_R0h
vpaddq $T1, $_R1, $_R1
vpaddq $T1h, $_R1h, $_R1h
vpaddq $T2, $_R2, $_R2
vpaddq $T2h, $_R2h, $_R2h
vpaddq $T3, $_R3, $_R3
vpaddq $T3h, $_R3h, $_R3h
vpaddq $T4, $_R4, $_R4
vpaddq $T4h, $_R4h, $_R4h
# Now handle carry bits from this addition
# Get mask of QWs whose 52-bit parts overflow
vpcmpuq \$6,.Lmask52x4(%rip),${_R0},%k1 # OP=nle (i.e. gt)
vpcmpuq \$6,.Lmask52x4(%rip),${_R0h},%k2
kmovb %k1,%r14d
kmovb %k2,%r13d
shl \$4,%r13b
or %r13b,%r14b
vpcmpuq \$6,.Lmask52x4(%rip),${_R1},%k1
vpcmpuq \$6,.Lmask52x4(%rip),${_R1h},%k2
kmovb %k1,%r13d
kmovb %k2,%r12d
shl \$4,%r12b
or %r12b,%r13b
vpcmpuq \$6,.Lmask52x4(%rip),${_R2},%k1
vpcmpuq \$6,.Lmask52x4(%rip),${_R2h},%k2
kmovb %k1,%r12d
kmovb %k2,%r11d
shl \$4,%r11b
or %r11b,%r12b
vpcmpuq \$6,.Lmask52x4(%rip),${_R3},%k1
vpcmpuq \$6,.Lmask52x4(%rip),${_R3h},%k2
kmovb %k1,%r11d
kmovb %k2,%r10d
shl \$4,%r10b
or %r10b,%r11b
vpcmpuq \$6,.Lmask52x4(%rip),${_R4},%k1
vpcmpuq \$6,.Lmask52x4(%rip),${_R4h},%k2
kmovb %k1,%r10d
kmovb %k2,%r9d
shl \$4,%r9b
or %r9b,%r10b
addb %r14b,%r14b
adcb %r13b,%r13b
adcb %r12b,%r12b
adcb %r11b,%r11b
adcb %r10b,%r10b
# Get mask of QWs whose 52-bit parts saturated
vpcmpuq \$0,.Lmask52x4(%rip),${_R0},%k1 # OP=eq
vpcmpuq \$0,.Lmask52x4(%rip),${_R0h},%k2
kmovb %k1,%r9d
kmovb %k2,%r8d
shl \$4,%r8b
or %r8b,%r9b
vpcmpuq \$0,.Lmask52x4(%rip),${_R1},%k1
vpcmpuq \$0,.Lmask52x4(%rip),${_R1h},%k2
kmovb %k1,%r8d
kmovb %k2,%edx
shl \$4,%dl
or %dl,%r8b
vpcmpuq \$0,.Lmask52x4(%rip),${_R2},%k1
vpcmpuq \$0,.Lmask52x4(%rip),${_R2h},%k2
kmovb %k1,%edx
kmovb %k2,%ecx
shl \$4,%cl
or %cl,%dl
vpcmpuq \$0,.Lmask52x4(%rip),${_R3},%k1
vpcmpuq \$0,.Lmask52x4(%rip),${_R3h},%k2
kmovb %k1,%ecx
kmovb %k2,%ebx
shl \$4,%bl
or %bl,%cl
vpcmpuq \$0,.Lmask52x4(%rip),${_R4},%k1
vpcmpuq \$0,.Lmask52x4(%rip),${_R4h},%k2
kmovb %k1,%ebx
kmovb %k2,%eax
shl \$4,%al
or %al,%bl
addb %r9b,%r14b
adcb %r8b,%r13b
adcb %dl,%r12b
adcb %cl,%r11b
adcb %bl,%r10b
xor %r9b,%r14b
xor %r8b,%r13b
xor %dl,%r12b
xor %cl,%r11b
xor %bl,%r10b
kmovb %r14d,%k1
shr \$4,%r14b
kmovb %r14d,%k2
kmovb %r13d,%k3
shr \$4,%r13b
kmovb %r13d,%k4
kmovb %r12d,%k5
shr \$4,%r12b
kmovb %r12d,%k6
kmovb %r11d,%k7
vpsubq .Lmask52x4(%rip), $_R0, ${_R0}{%k1}
vpsubq .Lmask52x4(%rip), $_R0h, ${_R0h}{%k2}
vpsubq .Lmask52x4(%rip), $_R1, ${_R1}{%k3}
vpsubq .Lmask52x4(%rip), $_R1h, ${_R1h}{%k4}
vpsubq .Lmask52x4(%rip), $_R2, ${_R2}{%k5}
vpsubq .Lmask52x4(%rip), $_R2h, ${_R2h}{%k6}
vpsubq .Lmask52x4(%rip), $_R3, ${_R3}{%k7}
vpandq .Lmask52x4(%rip), $_R0, $_R0
vpandq .Lmask52x4(%rip), $_R0h, $_R0h
vpandq .Lmask52x4(%rip), $_R1, $_R1
vpandq .Lmask52x4(%rip), $_R1h, $_R1h
vpandq .Lmask52x4(%rip), $_R2, $_R2
vpandq .Lmask52x4(%rip), $_R2h, $_R2h
vpandq .Lmask52x4(%rip), $_R3, $_R3
shr \$4,%r11b
kmovb %r11d,%k1
kmovb %r10d,%k2
shr \$4,%r10b
kmovb %r10d,%k3
vpsubq .Lmask52x4(%rip), $_R3h, ${_R3h}{%k1}
vpsubq .Lmask52x4(%rip), $_R4, ${_R4}{%k2}
vpsubq .Lmask52x4(%rip), $_R4h, ${_R4h}{%k3}
vpandq .Lmask52x4(%rip), $_R3h, $_R3h
vpandq .Lmask52x4(%rip), $_R4, $_R4
vpandq .Lmask52x4(%rip), $_R4h, $_R4h
___
}
$code.=<<___;
.text
.globl ossl_rsaz_amm52x40_x1_ifma256
.type ossl_rsaz_amm52x40_x1_ifma256,\@function,5
.align 32
ossl_rsaz_amm52x40_x1_ifma256:
.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
___
$code.=<<___ if ($win64);
lea -168(%rsp),%rsp # 16*10 + (8 bytes to get correct 16-byte SIMD alignment)
vmovdqa64 %xmm6, `0*16`(%rsp) # save non-volatile registers
vmovdqa64 %xmm7, `1*16`(%rsp)
vmovdqa64 %xmm8, `2*16`(%rsp)
vmovdqa64 %xmm9, `3*16`(%rsp)
vmovdqa64 %xmm10,`4*16`(%rsp)
vmovdqa64 %xmm11,`5*16`(%rsp)
vmovdqa64 %xmm12,`6*16`(%rsp)
vmovdqa64 %xmm13,`7*16`(%rsp)
vmovdqa64 %xmm14,`8*16`(%rsp)
vmovdqa64 %xmm15,`9*16`(%rsp)
.Lossl_rsaz_amm52x40_x1_ifma256_body:
___
$code.=<<___;
# 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, $R2_0h
vmovdqa64 $zero, $R3_0
vmovdqa64 $zero, $R3_0h
vmovdqa64 $zero, $R4_0
vmovdqa64 $zero, $R4_0h
xorl $acc0_0_low, $acc0_0_low
movq $b, $b_ptr # backup address of b
movq \$0xfffffffffffff, $mask52 # 52-bit mask
# Loop over 40 digits unrolled by 4
mov \$10, $iter
.align 32
.Lloop10:
___
foreach my $idx (0..3) {
&amm52x40_x1(0,8*$idx,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h,$R4_0,$R4_0h,$k0);
}
$code.=<<___;
lea `4*8`($b_ptr), $b_ptr
dec $iter
jne .Lloop10
___
&amm52x40_x1_norm($acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h,$R4_0,$R4_0h);
$code.=<<___;
vmovdqu64 $R0_0, `0*32`($res)
vmovdqu64 $R0_0h, `1*32`($res)
vmovdqu64 $R1_0, `2*32`($res)
vmovdqu64 $R1_0h, `3*32`($res)
vmovdqu64 $R2_0, `4*32`($res)
vmovdqu64 $R2_0h, `5*32`($res)
vmovdqu64 $R3_0, `6*32`($res)
vmovdqu64 $R3_0h, `7*32`($res)
vmovdqu64 $R4_0, `8*32`($res)
vmovdqu64 $R4_0h, `9*32`($res)
vzeroupper
lea (%rsp),%rax
.cfi_def_cfa_register %rax
___
$code.=<<___ if ($win64);
vmovdqa64 `0*16`(%rax),%xmm6
vmovdqa64 `1*16`(%rax),%xmm7
vmovdqa64 `2*16`(%rax),%xmm8
vmovdqa64 `3*16`(%rax),%xmm9
vmovdqa64 `4*16`(%rax),%xmm10
vmovdqa64 `5*16`(%rax),%xmm11
vmovdqa64 `6*16`(%rax),%xmm12
vmovdqa64 `7*16`(%rax),%xmm13
vmovdqa64 `8*16`(%rax),%xmm14
vmovdqa64 `9*16`(%rax),%xmm15
lea 168(%rsp),%rax
___
$code.=<<___;
mov 0(%rax),%r15
.cfi_restore %r15
mov 8(%rax),%r14
.cfi_restore %r14
mov 16(%rax),%r13
.cfi_restore %r13
mov 24(%rax),%r12
.cfi_restore %r12
mov 32(%rax),%rbp
.cfi_restore %rbp
mov 40(%rax),%rbx
.cfi_restore %rbx
lea 48(%rax),%rsp # restore rsp
.cfi_def_cfa %rsp,8
.Lossl_rsaz_amm52x40_x1_ifma256_epilogue:
ret
.cfi_endproc
.size ossl_rsaz_amm52x40_x1_ifma256, .-ossl_rsaz_amm52x40_x1_ifma256
___
$code.=<<___;
.data
.align 32
.Lmask52x4:
.quad 0xfffffffffffff
.quad 0xfffffffffffff
.quad 0xfffffffffffff
.quad 0xfffffffffffff
___
###############################################################################
# Dual Almost Montgomery Multiplication for 40-digit number in radix 2^52
#
# See description of ossl_rsaz_amm52x40_x1_ifma256() 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 ossl_rsaz_amm52x40_x2_ifma256(BN_ULONG out[2][40],
# const BN_ULONG a[2][40],
# const BN_ULONG b[2][40],
# const BN_ULONG m[2][40],
# const BN_ULONG k0[2]);
###############################################################################
$code.=<<___;
.text
.globl ossl_rsaz_amm52x40_x2_ifma256
.type ossl_rsaz_amm52x40_x2_ifma256,\@function,5
.align 32
ossl_rsaz_amm52x40_x2_ifma256:
.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
___
$code.=<<___ if ($win64);
lea -168(%rsp),%rsp
vmovdqa64 %xmm6, `0*16`(%rsp) # save non-volatile registers
vmovdqa64 %xmm7, `1*16`(%rsp)
vmovdqa64 %xmm8, `2*16`(%rsp)
vmovdqa64 %xmm9, `3*16`(%rsp)
vmovdqa64 %xmm10,`4*16`(%rsp)
vmovdqa64 %xmm11,`5*16`(%rsp)
vmovdqa64 %xmm12,`6*16`(%rsp)
vmovdqa64 %xmm13,`7*16`(%rsp)
vmovdqa64 %xmm14,`8*16`(%rsp)
vmovdqa64 %xmm15,`9*16`(%rsp)
.Lossl_rsaz_amm52x40_x2_ifma256_body:
___
$code.=<<___;
# 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, $R2_0h
vmovdqa64 $zero, $R3_0
vmovdqa64 $zero, $R3_0h
vmovdqa64 $zero, $R4_0
vmovdqa64 $zero, $R4_0h
vmovdqa64 $zero, $R0_1
vmovdqa64 $zero, $R0_1h
vmovdqa64 $zero, $R1_1
vmovdqa64 $zero, $R1_1h
vmovdqa64 $zero, $R2_1
vmovdqa64 $zero, $R2_1h
vmovdqa64 $zero, $R3_1
vmovdqa64 $zero, $R3_1h
vmovdqa64 $zero, $R4_1
vmovdqa64 $zero, $R4_1h
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 \$40, $iter
.align 32
.Lloop40:
___
&amm52x40_x1( 0, 0,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h,$R4_0,$R4_0h,"($k0)");
# 40*8 = offset of the next dimension in two-dimension array
&amm52x40_x1(40*8,40*8,$acc0_1,$R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1,$R2_1h,$R3_1,$R3_1h,$R4_1,$R4_1h,"8($k0)");
$code.=<<___;
lea 8($b_ptr), $b_ptr
dec $iter
jne .Lloop40
___
&amm52x40_x1_norm($acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$R2_0h,$R3_0,$R3_0h,$R4_0,$R4_0h);
&amm52x40_x1_norm($acc0_1,$R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1,$R2_1h,$R3_1,$R3_1h,$R4_1,$R4_1h);
$code.=<<___;
vmovdqu64 $R0_0, `0*32`($res)
vmovdqu64 $R0_0h, `1*32`($res)
vmovdqu64 $R1_0, `2*32`($res)
vmovdqu64 $R1_0h, `3*32`($res)
vmovdqu64 $R2_0, `4*32`($res)
vmovdqu64 $R2_0h, `5*32`($res)
vmovdqu64 $R3_0, `6*32`($res)
vmovdqu64 $R3_0h, `7*32`($res)
vmovdqu64 $R4_0, `8*32`($res)
vmovdqu64 $R4_0h, `9*32`($res)
vmovdqu64 $R0_1, `10*32`($res)
vmovdqu64 $R0_1h, `11*32`($res)
vmovdqu64 $R1_1, `12*32`($res)
vmovdqu64 $R1_1h, `13*32`($res)
vmovdqu64 $R2_1, `14*32`($res)
vmovdqu64 $R2_1h, `15*32`($res)
vmovdqu64 $R3_1, `16*32`($res)
vmovdqu64 $R3_1h, `17*32`($res)
vmovdqu64 $R4_1, `18*32`($res)
vmovdqu64 $R4_1h, `19*32`($res)
vzeroupper
lea (%rsp),%rax
.cfi_def_cfa_register %rax
___
$code.=<<___ if ($win64);
vmovdqa64 `0*16`(%rax),%xmm6
vmovdqa64 `1*16`(%rax),%xmm7
vmovdqa64 `2*16`(%rax),%xmm8
vmovdqa64 `3*16`(%rax),%xmm9
vmovdqa64 `4*16`(%rax),%xmm10
vmovdqa64 `5*16`(%rax),%xmm11
vmovdqa64 `6*16`(%rax),%xmm12
vmovdqa64 `7*16`(%rax),%xmm13
vmovdqa64 `8*16`(%rax),%xmm14
vmovdqa64 `9*16`(%rax),%xmm15
lea 168(%rsp),%rax
___
$code.=<<___;
mov 0(%rax),%r15
.cfi_restore %r15
mov 8(%rax),%r14
.cfi_restore %r14
mov 16(%rax),%r13
.cfi_restore %r13
mov 24(%rax),%r12
.cfi_restore %r12
mov 32(%rax),%rbp
.cfi_restore %rbp
mov 40(%rax),%rbx
.cfi_restore %rbx
lea 48(%rax),%rsp
.cfi_def_cfa %rsp,8
.Lossl_rsaz_amm52x40_x2_ifma256_epilogue:
ret
.cfi_endproc
.size ossl_rsaz_amm52x40_x2_ifma256, .-ossl_rsaz_amm52x40_x2_ifma256
___
}
###############################################################################
# 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.
# |red_table_idx1| and |red_table_idx2| are corresponding power indexes.
#
# Extracted value (output) is 2 40 digits numbers in 2^52 radix.
#
# void ossl_extract_multiplier_2x40_win5(BN_ULONG *red_Y,
# const BN_ULONG red_table[1 << EXP_WIN_SIZE][2][40],
# int red_table_idx1, int red_table_idx2);
#
# EXP_WIN_SIZE = 5
###############################################################################
{
# input parameters
my ($out,$red_tbl,$red_tbl_idx1,$red_tbl_idx2)=$win64 ? ("%rcx","%rdx","%r8", "%r9") : # Win64 order
("%rdi","%rsi","%rdx","%rcx"); # Unix order
my ($t0,$t1,$t2,$t3,$t4,$t5) = map("%ymm$_", (0..5));
my ($t6,$t7,$t8,$t9) = map("%ymm$_", (16..19));
my ($tmp,$cur_idx,$idx1,$idx2,$ones) = map("%ymm$_", (20..24));
my @t = ($t0,$t1,$t2,$t3,$t4,$t5,$t6,$t7,$t8,$t9);
my $t0xmm = $t0;
$t0xmm =~ s/%y/%x/;
sub get_table_value_consttime() {
my ($_idx,$_offset) = @_;
$code.=<<___;
vpxorq $cur_idx, $cur_idx, $cur_idx
.align 32
.Lloop_$_offset:
vpcmpq \$0, $cur_idx, $_idx, %k1 # mask of (idx == cur_idx)
___
foreach (0..9) {
$code.=<<___;
vmovdqu64 `$_offset+${_}*32`($red_tbl), $tmp # load data from red_tbl
vpblendmq $tmp, $t[$_], ${t[$_]}{%k1} # extract data when mask is not zero
___
}
$code.=<<___;
vpaddq $ones, $cur_idx, $cur_idx # increment cur_idx
addq \$`2*40*8`, $red_tbl
cmpq $red_tbl, %rax
jne .Lloop_$_offset
___
}
$code.=<<___;
.text
.align 32
.globl ossl_extract_multiplier_2x40_win5
.type ossl_extract_multiplier_2x40_win5,\@abi-omnipotent
ossl_extract_multiplier_2x40_win5:
.cfi_startproc
endbranch
vmovdqa64 .Lones(%rip), $ones # broadcast ones
vpbroadcastq $red_tbl_idx1, $idx1
vpbroadcastq $red_tbl_idx2, $idx2
leaq `(1<<5)*2*40*8`($red_tbl), %rax # holds end of the tbl
# backup red_tbl address
movq $red_tbl, %r10
# zeroing t0..n, cur_idx
vpxor $t0xmm, $t0xmm, $t0xmm
___
foreach (1..9) {
$code.="vmovdqa64 $t0, $t[$_] \n";
}
&get_table_value_consttime($idx1, 0);
foreach (0..9) {
$code.="vmovdqu64 $t[$_], `(0+$_)*32`($out) \n";
}
$code.="movq %r10, $red_tbl \n";
&get_table_value_consttime($idx2, 40*8);
foreach (0..9) {
$code.="vmovdqu64 $t[$_], `(10+$_)*32`($out) \n";
}
$code.=<<___;
ret
.cfi_endproc
.size ossl_extract_multiplier_2x40_win5, .-ossl_extract_multiplier_2x40_win5
___
$code.=<<___;
.data
.align 32
.Lones:
.quad 1,1,1,1
.Lzeros:
.quad 0,0,0,0
___
}
if ($win64) {
$rec="%rcx";
$frame="%rdx";
$context="%r8";
$disp="%r9";
$code.=<<___;
.extern __imp_RtlVirtualUnwind
.type rsaz_avx_handler,\@abi-omnipotent
.align 16
rsaz_avx_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 4(%r11),%r10d # HandlerData[1]
lea (%rsi,%r10),%r10 # epilogue label
cmp %r10,%rbx # context->Rip>=.Lepilogue
jae .Lcommon_seh_tail
mov 152($context),%rax # pull context->Rsp
lea (%rax),%rsi # %xmm save area
lea 512($context),%rdi # & context.Xmm6
mov \$20,%ecx # 10*sizeof(%xmm0)/sizeof(%rax)
.long 0xa548f3fc # cld; rep movsq
lea `48+168`(%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_avx_handler,.-rsaz_avx_handler
.section .pdata
.align 4
.rva .LSEH_begin_ossl_rsaz_amm52x40_x1_ifma256
.rva .LSEH_end_ossl_rsaz_amm52x40_x1_ifma256
.rva .LSEH_info_ossl_rsaz_amm52x40_x1_ifma256
.rva .LSEH_begin_ossl_rsaz_amm52x40_x2_ifma256
.rva .LSEH_end_ossl_rsaz_amm52x40_x2_ifma256
.rva .LSEH_info_ossl_rsaz_amm52x40_x2_ifma256
.section .xdata
.align 8
.LSEH_info_ossl_rsaz_amm52x40_x1_ifma256:
.byte 9,0,0,0
.rva rsaz_avx_handler
.rva .Lossl_rsaz_amm52x40_x1_ifma256_body,.Lossl_rsaz_amm52x40_x1_ifma256_epilogue
.LSEH_info_ossl_rsaz_amm52x40_x2_ifma256:
.byte 9,0,0,0
.rva rsaz_avx_handler
.rva .Lossl_rsaz_amm52x40_x2_ifma256_body,.Lossl_rsaz_amm52x40_x2_ifma256_epilogue
___
}
}}} else {{{ # fallback for old assembler
$code.=<<___;
.text
.globl ossl_rsaz_amm52x40_x1_ifma256
.globl ossl_rsaz_amm52x40_x2_ifma256
.globl ossl_extract_multiplier_2x40_win5
.type ossl_rsaz_amm52x40_x1_ifma256,\@abi-omnipotent
ossl_rsaz_amm52x40_x1_ifma256:
ossl_rsaz_amm52x40_x2_ifma256:
ossl_extract_multiplier_2x40_win5:
.byte 0x0f,0x0b # ud2
ret
.size ossl_rsaz_amm52x40_x1_ifma256, .-ossl_rsaz_amm52x40_x1_ifma256
___
}}}
$code =~ s/\`([^\`]*)\`/eval $1/gem;
print $code;
close STDOUT or die "error closing STDOUT: $!";

22
crypto/bn/bn_exp.c
View File

@ -1410,12 +1410,20 @@ int BN_mod_exp_mont_consttime_x2(BIGNUM *rr1, const BIGNUM *a1, const BIGNUM *p1
BN_MONT_CTX *mont2 = NULL;
if (ossl_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))) {
(((a1->top == 16) && (p1->top == 16) && (BN_num_bits(m1) == 1024) &&
(a2->top == 16) && (p2->top == 16) && (BN_num_bits(m2) == 1024)) ||
((a1->top == 24) && (p1->top == 24) && (BN_num_bits(m1) == 1536) &&
(a2->top == 24) && (p2->top == 24) && (BN_num_bits(m2) == 1536)) ||
((a1->top == 32) && (p1->top == 32) && (BN_num_bits(m1) == 2048) &&
(a2->top == 32) && (p2->top == 32) && (BN_num_bits(m2) == 2048)))) {
if (bn_wexpand(rr1, 16) == NULL)
int topn = a1->top;
/* Modulus bits of |m1| and |m2| are equal */
int mod_bits = BN_num_bits(m1);
if (bn_wexpand(rr1, topn) == NULL)
goto err;
if (bn_wexpand(rr2, 16) == NULL)
if (bn_wexpand(rr2, topn) == NULL)
goto err;
/* Ensure that montgomery contexts are initialized */
@ -1440,14 +1448,14 @@ int BN_mod_exp_mont_consttime_x2(BIGNUM *rr1, const BIGNUM *a1, const BIGNUM *p1
mont1->RR.d, mont1->n0[0],
rr2->d, a2->d, p2->d, m2->d,
mont2->RR.d, mont2->n0[0],
1024 /* factor bit size */);
mod_bits);
rr1->top = 16;
rr1->top = topn;
rr1->neg = 0;
bn_correct_top(rr1);
bn_check_top(rr1);
rr2->top = 16;
rr2->top = topn;
rr2->neg = 0;
bn_correct_top(rr2);
bn_check_top(rr2);

6
crypto/bn/build.info
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_exp_x2.c rsaz-avx512.s
rsaz-avx2.s rsaz_exp_x2.c rsaz-2k-avx512.s rsaz-3k-avx512.s rsaz-4k-avx512.s
IF[{- $config{target} !~ /^VC/ -}]
$BNASM_x86_64=asm/x86_64-gcc.c $BNASM_x86_64
ELSE
@ -160,7 +160,9 @@ 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[rsaz-2k-avx512.s]=asm/rsaz-2k-avx512.pl
GENERATE[rsaz-3k-avx512.s]=asm/rsaz-3k-avx512.pl
GENERATE[rsaz-4k-avx512.s]=asm/rsaz-4k-avx512.pl
GENERATE[bn-ia64.s]=asm/ia64.S
GENERATE[ia64-mont.s]=asm/ia64-mont.pl

404
crypto/bn/rsaz_exp_x2.c
View File

@ -1,6 +1,6 @@
/*
* Copyright 2020-2021 The OpenSSL Project Authors. All Rights Reserved.
* Copyright (c) 2020, Intel Corporation. All Rights Reserved.
* Copyright (c) 2020-2021, 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
@ -8,7 +8,8 @@
* https://www.openssl.org/source/license.html
*
*
* Originally written by Ilya Albrekht, Sergey Kirillov and Andrey Matyukov
* Originally written by Sergey Kirillov and Andrey Matyukov.
* Special thanks to Ilya Albrekht for his valuable hints.
* Intel Corporation
*
*/
@ -41,8 +42,12 @@ NON_EMPTY_TRANSLATION_UNIT
# 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);
/* Number of registers required to hold |digits_num| amount of qword digits */
# define NUMBER_OF_REGISTERS(digits_num, register_size) \
(((digits_num) * 64 + (register_size) - 1) / (register_size))
static ossl_inline uint64_t get_digit(const uint8_t *in, int in_len);
static ossl_inline void put_digit(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);
@ -54,37 +59,52 @@ 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:
* Naming conventions:
* amm = Almost Montgomery Multiplication
* ams = Almost Montgomery Squaring
* 52x20 - data represented as array of 20 digits in 52-bit radix
* 52xZZ - data represented as array of ZZ digits in 52-bit radix
* _x1_/_x2_ - 1 or 2 independent inputs/outputs
* _256 suffix - uses 256-bit (AVX512VL) registers
* _ifma256 - uses 256-bit wide IFMA ISA (AVX512_IFMA256)
*/
/*AMM = Almost Montgomery Multiplication. */
void ossl_rsaz_amm52x20_x1_256(BN_ULONG *res, const BN_ULONG *base,
const BN_ULONG *exp, const BN_ULONG *m,
BN_ULONG k0);
static 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 ossl_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 ossl_rsaz_amm52x20_x1_ifma256(BN_ULONG *res, const BN_ULONG *a,
const BN_ULONG *b, const BN_ULONG *m,
BN_ULONG k0);
void ossl_rsaz_amm52x20_x2_ifma256(BN_ULONG *out, const BN_ULONG *a,
const BN_ULONG *b, const BN_ULONG *m,
const BN_ULONG k0[2]);
void ossl_extract_multiplier_2x20_win5(BN_ULONG *red_Y,
const BN_ULONG *red_table,
int red_table_idx, int tbl_idx);
int red_table_idx1, int red_table_idx2);
void ossl_rsaz_amm52x30_x1_ifma256(BN_ULONG *res, const BN_ULONG *a,
const BN_ULONG *b, const BN_ULONG *m,
BN_ULONG k0);
void ossl_rsaz_amm52x30_x2_ifma256(BN_ULONG *out, const BN_ULONG *a,
const BN_ULONG *b, const BN_ULONG *m,
const BN_ULONG k0[2]);
void ossl_extract_multiplier_2x30_win5(BN_ULONG *red_Y,
const BN_ULONG *red_table,
int red_table_idx1, int red_table_idx2);
void ossl_rsaz_amm52x40_x1_ifma256(BN_ULONG *res, const BN_ULONG *a,
const BN_ULONG *b, const BN_ULONG *m,
BN_ULONG k0);
void ossl_rsaz_amm52x40_x2_ifma256(BN_ULONG *out, const BN_ULONG *a,
const BN_ULONG *b, const BN_ULONG *m,
const BN_ULONG k0[2]);
void ossl_extract_multiplier_2x40_win5(BN_ULONG *red_Y,
const BN_ULONG *red_table,
int red_table_idx1, int red_table_idx2);
static int RSAZ_mod_exp_x2_ifma256(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],
int modulus_bitsize);
/*
* Dual Montgomery modular exponentiation using prime moduli of the
@ -97,7 +117,10 @@ void ossl_extract_multiplier_2x20_win5(BN_ULONG *red_Y,
*
* Each moduli shall be |factor_size| bit size.
*
* NOTE: currently only 2x1024 case is supported.
* Supported cases:
* - 2x1024
* - 2x1536
* - 2x2048
*
* [out] res|i| - result of modular exponentiation: array of qword values
* in regular (2^64) radix. Size of array shall be enough
@ -126,6 +149,8 @@ int ossl_rsaz_mod_exp_avx512_x2(BN_ULONG *res1,
BN_ULONG k0_2,
int factor_size)
{
typedef void (*AMM)(BN_ULONG *res, const BN_ULONG *a,
const BN_ULONG *b, const BN_ULONG *m, BN_ULONG k0);
int ret = 0;
/*
@ -134,52 +159,60 @@ int ossl_rsaz_mod_exp_avx512_x2(BN_ULONG *res1,
*/
int exp_digits = number_of_digits(factor_size + 2, DIGIT_SIZE);
int coeff_pow = 4 * (DIGIT_SIZE * exp_digits - factor_size);
/* Number of YMM registers required to store exponent's digits */
int ymm_regs_num = NUMBER_OF_REGISTERS(exp_digits, 256 /* ymm bit size */);
/* Capacity of the register set (in qwords) to store exponent */
int regs_capacity = ymm_regs_num * 4;
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;
int storage_len_bytes = 7 * regs_capacity * sizeof(BN_ULONG)
+ 64 /* alignment */;
const BN_ULONG *exp[2] = {0};
BN_ULONG k0[2] = {0};
/* AMM = Almost Montgomery Multiplication */
AMM amm = NULL;
/* Only 1024-bit factor size is supported now */
switch (factor_size) {
case 1024:
amm = ossl_rsaz_amm52x20_x1_256;
exp_x2 = RSAZ_exp52x20_x2_256;
amm = ossl_rsaz_amm52x20_x1_ifma256;
break;
case 1536:
amm = ossl_rsaz_amm52x30_x1_ifma256;
break;
case 2048:
amm = ossl_rsaz_amm52x40_x1_ifma256;
break;
default:
goto err;
}
storage = (BN_ULONG *)OPENSSL_malloc(storage_len_bytes + 64);
storage = (BN_ULONG *)OPENSSL_malloc(storage_len_bytes);
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;
base2_red = storage_aligned + 1 * regs_capacity;
m1_red = storage_aligned + 2 * regs_capacity;
m2_red = storage_aligned + 3 * regs_capacity;
rr1_red = storage_aligned + 4 * regs_capacity;
rr2_red = storage_aligned + 5 * regs_capacity;
coeff_red = storage_aligned + 6 * regs_capacity;
/* 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);
to_words52(base1_red, regs_capacity, base1, factor_size);
to_words52(base2_red, regs_capacity, base2, factor_size);
to_words52(m1_red, regs_capacity, m1, factor_size);
to_words52(m2_red, regs_capacity, m2, factor_size);
to_words52(rr1_red, regs_capacity, rr1, factor_size);
to_words52(rr2_red, regs_capacity, rr2, factor_size);
/*
* Compute target domain Montgomery converters RR' for each modulus
@ -192,10 +225,10 @@ int ossl_rsaz_mod_exp_avx512_x2(BN_ULONG *res1,
* where
* k = 4 * (52 * digits52 - modlen)
* R = 2^(64 * ceil(modlen/64)) mod m
* RR = R^2 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
* EX/ 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 */
@ -213,13 +246,16 @@ int ossl_rsaz_mod_exp_avx512_x2(BN_ULONG *res1,
k0[0] = k0_1;
k0[1] = k0_2;
exp_x2(rr1_red, base1_red, exp, m1_red, rr1_red, k0);
/* Dual (2-exps in parallel) exponentiation */
ret = RSAZ_mod_exp_x2_ifma256(rr1_red, base1_red, exp, m1_red, rr1_red,
k0, factor_size);
if (!ret)
goto err;
/* 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);
@ -229,98 +265,156 @@ err:
}
/*
* Dual 1024-bit w-ary modular exponentiation using prime moduli of the same
* bit size using Almost Montgomery Multiplication, optimized with AVX512_IFMA
* ISA.
* Dual {1024,1536,2048}-bit w-ary modular exponentiation using prime moduli of
* the same bit size using Almost Montgomery Multiplication, optimized with
* AVX512_IFMA256 ISA.
*
* The parameter w (window size) = 5.
*
* [out] res - result of modular exponentiation: 2x20 qword
* [out] res - result of modular exponentiation: 2x{20,30,40} 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.
* [in] base - base (2x{20,30,40} qword values in 2^52 radix)
* [in] exp - array of 2 pointers to {16,24,32} 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] m - moduli (2x{20,30,40} qword values in 2^52 radix)
* [in] rr - Montgomery parameter for 2 moduli:
* RR(1024) = 2^2080 mod m.
* RR(1536) = 2^3120 mod m.
* RR(2048) = 2^4160 mod m.
* (2x{20,30,40} qword values in 2^52 radix)
* [in] k0 - Montgomery parameter for 2 moduli: k0 = -1/m mod 2^64
*
* \return (void).
*/
static 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])
int RSAZ_mod_exp_x2_ifma256(BN_ULONG *out,
const BN_ULONG *base,
const BN_ULONG *exp[2],
const BN_ULONG *m,
const BN_ULONG *rr,
const BN_ULONG k0[2],
int modulus_bitsize)
{
# 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 ossl_rsaz_amm52x20_x2_256
typedef void (*DAMM)(BN_ULONG *res, const BN_ULONG *a,
const BN_ULONG *b, const BN_ULONG *m,
const BN_ULONG k0[2]);
typedef void (*DEXTRACT)(BN_ULONG *res, const BN_ULONG *red_table,
int red_table_idx, int tbl_idx);
int ret = 0;
int idx;
/* Exponent window size */
int exp_win_size = 5;
int exp_win_mask = (1U << exp_win_size) - 1;
/*
* Number of digits (64-bit words) in redundant representation to handle
* modulus bits
*/
int red_digits = 0;
int exp_digits = 0;
BN_ULONG *storage = NULL;
BN_ULONG *storage_aligned = NULL;
int storage_len_bytes = 0;
/* Red(undant) result Y and multiplier X */
BN_ULONG *red_Y = NULL; /* [2][red_digits] */
BN_ULONG *red_X = NULL; /* [2][red_digits] */
/* Pre-computed table of base powers */
BN_ULONG *red_table = NULL; /* [1U << exp_win_size][2][red_digits] */
/* Expanded exponent */
BN_ULONG *expz = NULL; /* [2][exp_digits + 1] */
/* Dual AMM */
DAMM damm = NULL;
/* Extractor from red_table */
DEXTRACT extract = NULL;
/*
* Squaring is done using multiplication now. That can be a subject of
* optimization in future.
*/
# define DAMS(r,a,m,k0) \
ossl_rsaz_amm52x20_x2_256((r),(a),(a),(m),(k0))
# define DAMS(r,a,m,k0) damm((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];
switch (modulus_bitsize) {
case 1024:
red_digits = 20;
exp_digits = 16;
damm = ossl_rsaz_amm52x20_x2_ifma256;
extract = ossl_extract_multiplier_2x20_win5;
break;
case 1536:
/* Extended with 2 digits padding to avoid mask ops in high YMM register */
red_digits = 30 + 2;
exp_digits = 24;
damm = ossl_rsaz_amm52x30_x2_ifma256;
extract = ossl_extract_multiplier_2x30_win5;
break;
case 2048:
red_digits = 40;
exp_digits = 32;
damm = ossl_rsaz_amm52x40_x2_ifma256;
extract = ossl_extract_multiplier_2x40_win5;
break;
default:
goto err;
}
/* Allocate expanded exponent */
ALIGN64 BN_ULONG expz[2][EXP_DIGITS + 1];
storage_len_bytes = (2 * red_digits /* red_Y */
+ 2 * red_digits /* red_X */
+ 2 * red_digits * (1U << exp_win_size) /* red_table */
+ 2 * (exp_digits + 1)) /* expz */
* sizeof(BN_ULONG)
+ 64; /* alignment */
/* Pre-computed table of base powers */
ALIGN64 BN_ULONG red_table[1U << EXP_WIN_SIZE][2][RED_DIGITS];
storage = (BN_ULONG *)OPENSSL_zalloc(storage_len_bytes);
if (storage == NULL)
goto err;
storage_aligned = (BN_ULONG *)ALIGN_OF(storage, 64);
int idx;
memset(red_Y, 0, sizeof(red_Y));
memset(red_table, 0, sizeof(red_table));
memset(red_X, 0, sizeof(red_X));
red_Y = storage_aligned;
red_X = red_Y + 2 * red_digits;
red_table = red_X + 2 * red_digits;
expz = red_table + 2 * red_digits * (1U << exp_win_size);
/*
* 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);
red_X[0 * red_digits] = 1;
red_X[1 * red_digits] = 1;
damm(&red_table[0 * 2 * red_digits], (const BN_ULONG*)red_X, rr, m, k0);
damm(&red_table[1 * 2 * red_digits], 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);
for (idx = 1; idx < (int)((1U << exp_win_size) / 2); idx++) {
DAMS(&red_table[(2 * idx + 0) * 2 * red_digits],
&red_table[(1 * idx) * 2 * red_digits], m, k0);
damm(&red_table[(2 * idx + 1) * 2 * red_digits],
&red_table[(2 * idx) * 2 * red_digits],
&red_table[1 * 2 * red_digits], 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;
memcpy(&expz[0 * (exp_digits + 1)], exp[0], exp_digits * sizeof(BN_ULONG));
expz[1 * (exp_digits + 1) - 1] = 0;
memcpy(&expz[1 * (exp_digits + 1)], exp[1], exp_digits * sizeof(BN_ULONG));
expz[2 * (exp_digits + 1) - 1] = 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 rem = modulus_bitsize % 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_bit_no = modulus_bitsize - 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];
BN_ULONG red_table_idx_0 = expz[exp_chunk_no + 0 * (exp_digits + 1)];
BN_ULONG red_table_idx_1 = expz[exp_chunk_no + 1 * (exp_digits + 1)];
/*
* The function operates with fixed moduli sizes divisible by 64,
* thus table index here is always in supported range [0, EXP_WIN_SIZE).
@ -328,13 +422,10 @@ static void RSAZ_exp52x20_x2_256(BN_ULONG *out, /* [2][20] */
red_table_idx_0 >>= exp_chunk_shift;
red_table_idx_1 >>= exp_chunk_shift;
ossl_extract_multiplier_2x20_win5(red_Y[0], (const BN_ULONG*)red_table,
(int)red_table_idx_0, 0);
ossl_extract_multiplier_2x20_win5(red_Y[1], (const BN_ULONG*)red_table,
(int)red_table_idx_1, 1);
extract(&red_Y[0 * red_digits], (const BN_ULONG*)red_table, (int)red_table_idx_0, (int)red_table_idx_1);
/* Process other exp windows */
for (exp_bit_no -= EXP_WIN_SIZE; exp_bit_no >= 0; exp_bit_no -= EXP_WIN_SIZE) {
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;
@ -342,43 +433,37 @@ static void RSAZ_exp52x20_x2_256(BN_ULONG *out, /* [2][20] */
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 = expz[exp_chunk_no + 0 * (exp_digits + 1)];
T = expz[exp_chunk_no + 1 + 0 * (exp_digits + 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) {
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;
ossl_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 = expz[exp_chunk_no + 1 * (exp_digits + 1)];
T = expz[exp_chunk_no + 1 + 1 * (exp_digits + 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) {
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;
ossl_extract_multiplier_2x20_win5(red_X[1],
(const BN_ULONG*)red_table,
(int)red_table_idx_1, 1);
}
extract(&red_X[0 * red_digits], (const BN_ULONG*)red_table, (int)red_table_idx_0, (int)red_table_idx_1);
}
/* Series of squaring */
@ -388,43 +473,46 @@ static void RSAZ_exp52x20_x2_256(BN_ULONG *out, /* [2][20] */
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);
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.
* may be (modulus_bitsize + 1), 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 [1] for details.
*
* [1] Gueron, S. Efficient software implementations of modular exponentiation.
* DOI: 10.1007/s13389-012-0031-5
*/
/* 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);
memset(red_X, 0, 2 * red_digits * sizeof(BN_ULONG));
red_X[0 * red_digits] = 1;
red_X[1 * red_digits] = 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));
ret = 1;
# undef DAMS
# undef DAMM
# undef EXP_DIGITS
# undef RED_DIGITS
# undef EXP_WIN_MASK
# undef EXP_WIN_SIZE
# undef BITSIZE_MODULUS
err:
if (storage != NULL) {
/* Clear whole storage */
OPENSSL_cleanse(storage, storage_len_bytes);
OPENSSL_free(storage);
}
#undef DAMS
return ret;
}
static ossl_inline uint64_t get_digit52(const uint8_t *in, int in_len)
static ossl_inline uint64_t get_digit(const uint8_t *in, int in_len)
{
uint64_t digit = 0;
assert(in != NULL);
assert(in_len <= 8);
for (; in_len > 0; in_len--) {
digit <<= 8;
@ -458,17 +546,17 @@ static void to_words52(BN_ULONG *out, int out_len,
}
if (in_bitsize > DIGIT_SIZE) {
uint64_t digit = get_digit52(in_str, 7);
uint64_t digit = get_digit(in_str, 7);
out[0] = digit & DIGIT_MASK;
in_str += 6;
in_bitsize -= DIGIT_SIZE;
digit = get_digit52(in_str, BITS2WORD8_SIZE(in_bitsize));
digit = get_digit(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[0] = get_digit(in_str, BITS2WORD8_SIZE(in_bitsize));
out++;
out_len--;
}
@ -480,12 +568,13 @@ static void to_words52(BN_ULONG *out, int out_len,
}
}
static ossl_inline void put_digit52(uint8_t *pStr, int strLen, uint64_t digit)
static ossl_inline void put_digit(uint8_t *out, int out_len, uint64_t digit)
{
assert(pStr != NULL);
assert(out != NULL);
assert(out_len <= 8);
for (; strLen > 0; strLen--) {
*pStr++ = (uint8_t)(digit & 0xFF);
for (; out_len > 0; out_len--) {
*out++ = (uint8_t)(digit & 0xFF);
digit >>= 8;
}
}
@ -508,7 +597,8 @@ static void from_words52(BN_ULONG *out, int out_bitsize, const BN_ULONG *in)
{
uint8_t *out_str = (uint8_t *)out;
for (; out_bitsize >= (2 * DIGIT_SIZE); out_bitsize -= (2 * DIGIT_SIZE), in += 2) {
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;
@ -516,13 +606,13 @@ static void from_words52(BN_ULONG *out, int out_bitsize, const BN_ULONG *in)
}
if (out_bitsize > DIGIT_SIZE) {
put_digit52(out_str, 7, in[0]);
put_digit(out_str, 7, in[0]);
out_str += 6;
out_bitsize -= DIGIT_SIZE;
put_digit52(out_str, BITS2WORD8_SIZE(out_bitsize),
put_digit(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]);
put_digit(out_str, BITS2WORD8_SIZE(out_bitsize), in[0]);
}
}
}

9
test/exptest.c
View File

@ -223,11 +223,12 @@ static int test_mod_exp_x2(int idx)
BIGNUM *m2 = NULL;
int factor_size = 0;
/*
* Currently only 1024-bit factor size is supported.
*/
if (idx <= 100)
factor_size = 1024;
else if (idx <= 200)
factor_size = 1536;
else if (idx <= 300)
factor_size = 2048;
if (!TEST_ptr(ctx = BN_CTX_new()))
goto err;
@ -303,6 +304,6 @@ 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);
ADD_ALL_TESTS(test_mod_exp_x2, 300);
return 1;
}