diff --git a/CHANGES.md b/CHANGES.md index 940d450fdf..7724011022 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -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 ----------- diff --git a/crypto/bn/asm/rsaz-avx512.pl b/crypto/bn/asm/rsaz-2k-avx512.pl similarity index 70% rename from crypto/bn/asm/rsaz-avx512.pl rename to crypto/bn/asm/rsaz-2k-avx512.pl index d031caa88e..aec5470387 100644 --- a/crypto/bn/asm/rsaz-avx512.pl +++ b/crypto/bn/asm/rsaz-2k-avx512.pl @@ -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 ___ }}} diff --git a/crypto/bn/asm/rsaz-3k-avx512.pl b/crypto/bn/asm/rsaz-3k-avx512.pl new file mode 100644 index 0000000000..e294afd294 --- /dev/null +++ b/crypto/bn/asm/rsaz-3k-avx512.pl @@ -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: $!"; diff --git a/crypto/bn/asm/rsaz-4k-avx512.pl b/crypto/bn/asm/rsaz-4k-avx512.pl new file mode 100644 index 0000000000..fb5bf10198 --- /dev/null +++ b/crypto/bn/asm/rsaz-4k-avx512.pl @@ -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: $!"; diff --git a/crypto/bn/bn_exp.c b/crypto/bn/bn_exp.c index 5329cd12a9..bb20e1683e 100644 --- a/crypto/bn/bn_exp.c +++ b/crypto/bn/bn_exp.c @@ -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); diff --git a/crypto/bn/build.info b/crypto/bn/build.info index 9330274aef..b0fcb7ab97 100644 --- a/crypto/bn/build.info +++ b/crypto/bn/build.info @@ -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 diff --git a/crypto/bn/rsaz_exp_x2.c b/crypto/bn/rsaz_exp_x2.c index b7d11180f8..93f104f68e 100644 --- a/crypto/bn/rsaz_exp_x2.c +++ b/crypto/bn/rsaz_exp_x2.c @@ -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]); } } } diff --git a/test/exptest.c b/test/exptest.c index 84d972afe3..675984c8cb 100644 --- a/test/exptest.c +++ b/test/exptest.c @@ -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; }