Dual 1024-bit exponentiation optimization for Intel IceLake CPU

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

Reviewed-by: Paul Dale <pauli@openssl.org>
Reviewed-by: Matt Caswell <matt@openssl.org>
(Merged from https://github.com/openssl/openssl/pull/13750)

7 files changed, 1400 insertions, 10 deletions

diff --git a/crypto/bn/asm/rsaz-avx512.pl b/crypto/bn/asm/rsaz-avx512.pl
new file mode 100644
index 0000000000..063b9d6b5e
--- /dev/null
+++ b/crypto/bn/asm/rsaz-avx512.pl
@@ -0,0 +1,743 @@
+# Copyright 2020 The OpenSSL Project Authors. All Rights Reserved.
+# Copyright (c) 2020, Intel Corporation. All Rights Reserved.
+#
+# Licensed under the Apache License 2.0 (the "License").  You may not use
+# this file except in compliance with the License.  You can obtain a copy
+# in the file LICENSE in the source distribution or at
+# https://www.openssl.org/source/license.html
+#
+#
+# Originally written by Ilya Albrekht, Sergey Kirillov and Andrey Matyukov
+# Intel Corporation
+#
+# December 2020
+#
+# Initial release.
+#
+# Implementation utilizes 256-bit (ymm) registers to avoid frequency scaling issues.
+#
+# IceLake-Client @ 1.3GHz
+# |---------+----------------------+--------------+-------------|
+# |         | OpenSSL 3.0.0-alpha9 | this         | Unit        |
+# |---------+----------------------+--------------+-------------|
+# | rsa2048 | 2 127 659            | 1 015 625    | cycles/sign |
+# |         | 611                  | 1280 / +109% | sign/s      |
+# |---------+----------------------+--------------+-------------|
+#
+
+# $output is the last argument if it looks like a file (it has an extension)
+# $flavour is the first argument if it doesn't look like a file
+$output = $#ARGV >= 0 && $ARGV[$#ARGV] =~ m|\.\w+$| ? pop : undef;
+$flavour = $#ARGV >= 0 && $ARGV[0] !~ m|\.| ? shift : undef;
+
+$win64=0; $win64=1 if ($flavour =~ /[nm]asm|mingw64/ || $output =~ /\.asm$/);
+$avx512ifma=0;
+
+$0 =~ m/(.*[\/\\])[^\/\\]+$/; $dir=$1;
+( $xlate="${dir}x86_64-xlate.pl" and -f $xlate ) or
+( $xlate="${dir}../../perlasm/x86_64-xlate.pl" and -f $xlate) or
+die "can't locate x86_64-xlate.pl";
+
+if (`$ENV{CC} -Wa,-v -c -o /dev/null -x assembler /dev/null 2>&1`
+        =~ /GNU assembler version ([2-9]\.[0-9]+)/) {
+    $avx512ifma = ($1>=2.26);
+}
+
+if (!$avx512 && $win64 && ($flavour =~ /nasm/ || $ENV{ASM} =~ /nasm/) &&
+       `nasm -v 2>&1` =~ /NASM version ([2-9]\.[0-9]+)(?:\.([0-9]+))?/) {
+    $avx512ifma = ($1==2.11 && $2>=8) + ($1>=2.12);
+}
+
+if (!$avx512 && `$ENV{CC} -v 2>&1` =~ /((?:clang|LLVM) version|.*based on LLVM) ([0-9]+\.[0-9]+)/) {
+    $avx512ifma = ($2>=6.0);
+}
+
+open OUT,"| \"$^X\" \"$xlate\" $flavour \"$output\""
+    or die "can't call $xlate: $!";
+*STDOUT=*OUT;
+
+if ($avx512ifma>0) {{{
+@_6_args_universal_ABI = ("%rdi","%rsi","%rdx","%rcx","%r8","%r9");
+
+$code.=<<___;
+.extern OPENSSL_ia32cap_P
+.globl  rsaz_avx512ifma_eligible
+.type   rsaz_avx512ifma_eligible,\@abi-omnipotent
+.align  32
+rsaz_avx512ifma_eligible:
+    mov OPENSSL_ia32cap_P+8(%rip), %ecx
+    xor %eax,%eax
+    and \$`1<<31|1<<21|1<<17|1<<16`, %ecx     # avx512vl + avx512ifma + avx512dq + avx512f
+    cmp \$`1<<31|1<<21|1<<17|1<<16`, %ecx
+    cmove %ecx,%eax
+    ret
+.size   rsaz_avx512ifma_eligible, .-rsaz_avx512ifma_eligible
+___
+
+###############################################################################
+# Almost Montgomery Multiplication (AMM) for 20-digit number in radix 2^52.
+#
+# AMM is defined as presented in the paper
+# "Efficient Software Implementations of Modular Exponentiation" by Shay Gueron.
+#
+# The input and output are presented in 2^52 radix domain, i.e.
+#   |res|, |a|, |b|, |m| are arrays of 20 64-bit qwords with 12 high bits zeroed.
+#   |k0| is a Montgomery coefficient, which is here k0 = -1/m mod 2^64
+#        (note, the implementation counts only 52 bits from it).
+#
+# NB: the AMM implementation does not perform "conditional" subtraction step as
+# specified in the original algorithm as according to the paper "Enhanced Montgomery
+# Multiplication" by Shay Gueron (see Lemma 1), the result will be always < 2*2^1024
+# and can be used as a direct input to the next AMM iteration.
+# This post-condition is true, provided the correct parameter |s| is choosen, i.e.
+# s >= n + 2 * k, which matches our case: 1040 > 1024 + 2 * 1.
+#
+# void RSAZ_amm52x20_x1_256(BN_ULONG *res,
+#                           const BN_ULONG *a,
+#                           const BN_ULONG *b,
+#                           const BN_ULONG *m,
+#                           BN_ULONG k0);
+###############################################################################
+{
+# input parameters ("%rdi","%rsi","%rdx","%rcx","%r8")
+my ($res,$a,$b,$m,$k0) = @_6_args_universal_ABI;
+
+my $mask52     = "%rax";
+my $acc0_0     = "%r9";
+my $acc0_0_low = "%r9d";
+my $acc0_1     = "%r15";
+my $acc0_1_low = "%r15d";
+my $b_ptr      = "%r11";
+
+my $iter = "%ebx";
+
+my $zero = "%ymm0";
+my ($R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0) = ("%ymm1", map("%ymm$_",(16..19)));
+my ($R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1) = ("%ymm2", map("%ymm$_",(20..23)));
+my $Bi = "%ymm3";
+my $Yi = "%ymm4";
+
+# Registers mapping for normalization.
+# We can reuse Bi, Yi registers here.
+my $TMP = $Bi;
+my $mask52x4 = $Yi;
+my ($T0,$T0h,$T1,$T1h,$T2) = map("%ymm$_", (24..28));
+
+sub amm52x20_x1() {
+# _data_offset - offset in the |a| or |m| arrays pointing to the beginning
+#                of data for corresponding AMM operation;
+# _b_offset    - offset in the |b| array pointing to the next qword digit;
+my ($_data_offset,$_b_offset,$_acc,$_R0,$_R0h,$_R1,$_R1h,$_R2,$_k0) = @_;
+my $_R0_xmm = $_R0 =~ s/%y/%x/r;
+$code.=<<___;
+    movq    $_b_offset($b_ptr), %r13             # b[i]
+
+    vpbroadcastq    %r13, $Bi                    # broadcast b[i]
+    movq    $_data_offset($a), %rdx
+    mulx    %r13, %r13, %r12                     # a[0]*b[i] = (t0,t2)
+    addq    %r13, $_acc                          # acc += t0
+    movq    %r12, %r10
+    adcq    \$0, %r10                            # t2 += CF
+
+    movq    $_k0, %r13
+    imulq   $_acc, %r13                          # acc * k0
+    andq    $mask52, %r13                        # yi = (acc * k0) & mask52
+
+    vpbroadcastq    %r13, $Yi                    # broadcast y[i]
+    movq    $_data_offset($m), %rdx
+    mulx    %r13, %r13, %r12                     # yi * m[0] = (t0,t1)
+    addq    %r13, $_acc                          # acc += t0
+    adcq    %r12, %r10                           # t2 += (t1 + CF)
+
+    shrq    \$52, $_acc
+    salq    \$12, %r10
+    or      %r10, $_acc                          # acc = ((acc >> 52) | (t2 << 12))
+
+    vpmadd52luq `$_data_offset+64*0`($a), $Bi, $_R0
+    vpmadd52luq `$_data_offset+64*0+32`($a), $Bi, $_R0h
+    vpmadd52luq `$_data_offset+64*1`($a), $Bi, $_R1
+    vpmadd52luq `$_data_offset+64*1+32`($a), $Bi, $_R1h
+    vpmadd52luq `$_data_offset+64*2`($a), $Bi, $_R2
+
+    vpmadd52luq `$_data_offset+64*0`($m), $Yi, $_R0
+    vpmadd52luq `$_data_offset+64*0+32`($m), $Yi, $_R0h
+    vpmadd52luq `$_data_offset+64*1`($m), $Yi, $_R1
+    vpmadd52luq `$_data_offset+64*1+32`($m), $Yi, $_R1h
+    vpmadd52luq `$_data_offset+64*2`($m), $Yi, $_R2
+
+    # Shift accumulators right by 1 qword, zero extending the highest one
+    valignq     \$1, $_R0, $_R0h, $_R0
+    valignq     \$1, $_R0h, $_R1, $_R0h
+    valignq     \$1, $_R1, $_R1h, $_R1
+    valignq     \$1, $_R1h, $_R2, $_R1h
+    valignq     \$1, $_R2, $zero, $_R2
+
+    vmovq   $_R0_xmm, %r13
+    addq    %r13, $_acc    # acc += R0[0]
+
+    vpmadd52huq `$_data_offset+64*0`($a), $Bi, $_R0
+    vpmadd52huq `$_data_offset+64*0+32`($a), $Bi, $_R0h
+    vpmadd52huq `$_data_offset+64*1`($a), $Bi, $_R1
+    vpmadd52huq `$_data_offset+64*1+32`($a), $Bi, $_R1h
+    vpmadd52huq `$_data_offset+64*2`($a), $Bi, $_R2
+
+    vpmadd52huq `$_data_offset+64*0`($m), $Yi, $_R0
+    vpmadd52huq `$_data_offset+64*0+32`($m), $Yi, $_R0h
+    vpmadd52huq `$_data_offset+64*1`($m), $Yi, $_R1
+    vpmadd52huq `$_data_offset+64*1+32`($m), $Yi, $_R1h
+    vpmadd52huq `$_data_offset+64*2`($m), $Yi, $_R2
+___
+}
+
+# Normalization routine: handles carry bits in R0..R2 QWs and
+# gets R0..R2 back to normalized 2^52 representation.
+#
+# Uses %r8-14,%e[bcd]x
+sub amm52x20_x1_norm {
+my ($_acc,$_R0,$_R0h,$_R1,$_R1h,$_R2) = @_;
+$code.=<<___;
+    # Put accumulator to low qword in R0
+    vpbroadcastq    $_acc, $TMP
+    vpblendd \$3, $TMP, $_R0, $_R0
+
+    # Extract "carries" (12 high bits) from each QW of R0..R2
+    # Save them to LSB of QWs in T0..T2
+    vpsrlq    \$52, $_R0,   $T0
+    vpsrlq    \$52, $_R0h,  $T0h
+    vpsrlq    \$52, $_R1,   $T1
+    vpsrlq    \$52, $_R1h,  $T1h
+    vpsrlq    \$52, $_R2,   $T2
+
+    # "Shift left" T0..T2 by 1 QW
+    valignq \$3, $T1h,  $T2,  $T2
+    valignq \$3, $T1,   $T1h, $T1h
+    valignq \$3, $T0h,  $T1,  $T1
+    valignq \$3, $T0,   $T0h, $T0h
+    valignq \$3, $zero, $T0,  $T0
+
+    # Drop "carries" from R0..R2 QWs
+    vpandq    $mask52x4, $_R0,  $_R0
+    vpandq    $mask52x4, $_R0h, $_R0h
+    vpandq    $mask52x4, $_R1,  $_R1
+    vpandq    $mask52x4, $_R1h, $_R1h
+    vpandq    $mask52x4, $_R2,  $_R2
+
+    # Sum R0..R2 with corresponding adjusted carries
+    vpaddq  $T0,  $_R0,  $_R0
+    vpaddq  $T0h, $_R0h, $_R0h
+    vpaddq  $T1,  $_R1,  $_R1
+    vpaddq  $T1h, $_R1h, $_R1h
+    vpaddq  $T2,  $_R2,  $_R2
+
+    # Now handle carry bits from this addition
+    # Get mask of QWs which 52-bit parts overflow...
+    vpcmpuq   \$1, $_R0,  $mask52x4, %k1 # OP=lt
+    vpcmpuq   \$1, $_R0h, $mask52x4, %k2
+    vpcmpuq   \$1, $_R1,  $mask52x4, %k3
+    vpcmpuq   \$1, $_R1h, $mask52x4, %k4
+    vpcmpuq   \$1, $_R2,  $mask52x4, %k5
+    kmovb   %k1, %r14d                   # k1
+    kmovb   %k2, %r13d                   # k1h
+    kmovb   %k3, %r12d                   # k2
+    kmovb   %k4, %r11d                   # k2h
+    kmovb   %k5, %r10d                   # k3
+
+    # ...or saturated
+    vpcmpuq   \$0, $_R0,  $mask52x4, %k1 # OP=eq
+    vpcmpuq   \$0, $_R0h, $mask52x4, %k2
+    vpcmpuq   \$0, $_R1,  $mask52x4, %k3
+    vpcmpuq   \$0, $_R1h, $mask52x4, %k4
+    vpcmpuq   \$0, $_R2,  $mask52x4, %k5
+    kmovb   %k1, %r9d                    # k4
+    kmovb   %k2, %r8d                    # k4h
+    kmovb   %k3, %ebx                    # k5
+    kmovb   %k4, %ecx                    # k5h
+    kmovb   %k5, %edx                    # k6
+
+    # Get mask of QWs where carries shall be propagated to.
+    # Merge 4-bit masks to 8-bit values to use add with carry.
+    shl   \$4, %r13b
+    or    %r13b, %r14b
+    shl   \$4, %r11b
+    or    %r11b, %r12b
+
+    add   %r14b, %r14b
+    adc   %r12b, %r12b
+    adc   %r10b, %r10b
+
+    shl   \$4, %r8b
+    or    %r8b,%r9b
+    shl   \$4, %cl
+    or    %cl, %bl
+
+    add   %r9b, %r14b
+    adc   %bl, %r12b
+    adc   %dl, %r10b
+
+    xor   %r9b, %r14b
+    xor   %bl, %r12b
+    xor   %dl, %r10b
+
+    kmovb   %r14d, %k1
+    shr     \$4, %r14b
+    kmovb   %r14d, %k2
+    kmovb   %r12d, %k3
+    shr     \$4, %r12b
+    kmovb   %r12d, %k4
+    kmovb   %r10d, %k5
+
+    # Add carries according to the obtained mask
+    vpsubq  $mask52x4, $_R0,  ${_R0}{%k1}
+    vpsubq  $mask52x4, $_R0h, ${_R0h}{%k2}
+    vpsubq  $mask52x4, $_R1,  ${_R1}{%k3}
+    vpsubq  $mask52x4, $_R1h, ${_R1h}{%k4}
+    vpsubq  $mask52x4, $_R2,  ${_R2}{%k5}
+
+    vpandq   $mask52x4, $_R0,  $_R0
+    vpandq   $mask52x4, $_R0h, $_R0h
+    vpandq   $mask52x4, $_R1,  $_R1
+    vpandq   $mask52x4, $_R1h, $_R1h
+    vpandq   $mask52x4, $_R2,  $_R2
+___
+}
+
+$code.=<<___;
+.text
+
+.globl  RSAZ_amm52x20_x1_256
+.type   RSAZ_amm52x20_x1_256,\@function,5
+.align 32
+RSAZ_amm52x20_x1_256:
+.cfi_startproc
+    endbranch
+    push    %rbx
+.cfi_push   %rbx
+    push    %rbp
+.cfi_push   %rbp
+    push    %r12
+.cfi_push   %r12
+    push    %r13
+.cfi_push   %r13
+    push    %r14
+.cfi_push   %r14
+    push    %r15
+.cfi_push   %r15
+.Lrsaz_amm52x20_x1_256_body:
+
+    # Zeroing accumulators
+    vpxord   $zero, $zero, $zero
+    vmovdqa64   $zero, $R0_0
+    vmovdqa64   $zero, $R0_0h
+    vmovdqa64   $zero, $R1_0
+    vmovdqa64   $zero, $R1_0h
+    vmovdqa64   $zero, $R2_0
+
+    xorl    $acc0_0_low, $acc0_0_low
+
+    movq    $b, $b_ptr                       # backup address of b
+    movq    \$0xfffffffffffff, $mask52       # 52-bit mask
+
+    # Loop over 20 digits unrolled by 4
+    mov     \$5, $iter
+
+.align 32
+.Lloop5:
+___
+    foreach my $idx (0..3) {
+        &amm52x20_x1(0,8*$idx,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,$k0);
+    }
+$code.=<<___;
+    lea    `4*8`($b_ptr), $b_ptr
+    dec    $iter
+    jne    .Lloop5
+
+    vmovdqa64   .Lmask52x4(%rip), $mask52x4
+___
+    &amm52x20_x1_norm($acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0);
+$code.=<<___;
+
+    vmovdqu64   $R0_0, ($res)
+    vmovdqu64   $R0_0h, 32($res)
+    vmovdqu64   $R1_0, 64($res)
+    vmovdqu64   $R1_0h, 96($res)
+    vmovdqu64   $R2_0, 128($res)
+
+    vzeroupper
+    mov  0(%rsp),%r15
+.cfi_restore    %r15
+    mov  8(%rsp),%r14
+.cfi_restore    %r14
+    mov  16(%rsp),%r13
+.cfi_restore    %r13
+    mov  24(%rsp),%r12
+.cfi_restore    %r12
+    mov  32(%rsp),%rbp
+.cfi_restore    %rbp
+    mov  40(%rsp),%rbx
+.cfi_restore    %rbx
+    lea  48(%rsp),%rsp
+.cfi_adjust_cfa_offset  -48
+.Lrsaz_amm52x20_x1_256_epilogue:
+    ret
+.cfi_endproc
+.size   RSAZ_amm52x20_x1_256, .-RSAZ_amm52x20_x1_256
+___
+
+$code.=<<___;
+.data
+.align 32
+.Lmask52x4:
+    .quad   0xfffffffffffff
+    .quad   0xfffffffffffff
+    .quad   0xfffffffffffff
+    .quad   0xfffffffffffff
+___
+
+###############################################################################
+# Dual Almost Montgomery Multiplication for 20-digit number in radix 2^52
+#
+# See description of RSAZ_amm52x20_x1_256() above for details about Almost
+# Montgomery Multiplication algorithm and function input parameters description.
+#
+# This function does two AMMs for two independent inputs, hence dual.
+#
+# void RSAZ_amm52x20_x2_256(BN_ULONG out[2][20],
+#                           const BN_ULONG a[2][20],
+#                           const BN_ULONG b[2][20],
+#                           const BN_ULONG m[2][20],
+#                           const BN_ULONG k0[2]);
+###############################################################################
+
+$code.=<<___;
+.text
+
+.globl  RSAZ_amm52x20_x2_256
+.type   RSAZ_amm52x20_x2_256,\@function,5
+.align 32
+RSAZ_amm52x20_x2_256:
+.cfi_startproc
+    endbranch
+    push    %rbx
+.cfi_push   %rbx
+    push    %rbp
+.cfi_push   %rbp
+    push    %r12
+.cfi_push   %r12
+    push    %r13
+.cfi_push   %r13
+    push    %r14
+.cfi_push   %r14
+    push    %r15
+.cfi_push   %r15
+.Lrsaz_amm52x20_x2_256_body:
+
+    # Zeroing accumulators
+    vpxord   $zero, $zero, $zero
+    vmovdqa64   $zero, $R0_0
+    vmovdqa64   $zero, $R0_0h
+    vmovdqa64   $zero, $R1_0
+    vmovdqa64   $zero, $R1_0h
+    vmovdqa64   $zero, $R2_0
+    vmovdqa64   $zero, $R0_1
+    vmovdqa64   $zero, $R0_1h
+    vmovdqa64   $zero, $R1_1
+    vmovdqa64   $zero, $R1_1h
+    vmovdqa64   $zero, $R2_1
+
+    xorl    $acc0_0_low, $acc0_0_low
+    xorl    $acc0_1_low, $acc0_1_low
+
+    movq    $b, $b_ptr                       # backup address of b
+    movq    \$0xfffffffffffff, $mask52       # 52-bit mask
+
+    mov    \$20, $iter
+
+.align 32
+.Lloop20:
+___
+    &amm52x20_x1(   0,   0,$acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0,"($k0)");
+    # 20*8 = offset of the next dimension in two-dimension array
+    &amm52x20_x1(20*8,20*8,$acc0_1,$R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1,"8($k0)");
+$code.=<<___;
+    lea    8($b_ptr), $b_ptr
+    dec    $iter
+    jne    .Lloop20
+
+    vmovdqa64   .Lmask52x4(%rip), $mask52x4
+___
+    &amm52x20_x1_norm($acc0_0,$R0_0,$R0_0h,$R1_0,$R1_0h,$R2_0);
+    &amm52x20_x1_norm($acc0_1,$R0_1,$R0_1h,$R1_1,$R1_1h,$R2_1);
+$code.=<<___;
+
+    vmovdqu64   $R0_0, ($res)
+    vmovdqu64   $R0_0h, 32($res)
+    vmovdqu64   $R1_0, 64($res)
+    vmovdqu64   $R1_0h, 96($res)
+    vmovdqu64   $R2_0, 128($res)
+
+    vmovdqu64   $R0_1, 160($res)
+    vmovdqu64   $R0_1h, 192($res)
+    vmovdqu64   $R1_1, 224($res)
+    vmovdqu64   $R1_1h, 256($res)
+    vmovdqu64   $R2_1, 288($res)
+
+    vzeroupper
+    mov  0(%rsp),%r15
+.cfi_restore    %r15
+    mov  8(%rsp),%r14
+.cfi_restore    %r14
+    mov  16(%rsp),%r13
+.cfi_restore    %r13
+    mov  24(%rsp),%r12
+.cfi_restore    %r12
+    mov  32(%rsp),%rbp
+.cfi_restore    %rbp
+    mov  40(%rsp),%rbx
+.cfi_restore    %rbx
+    lea  48(%rsp),%rsp
+.cfi_adjust_cfa_offset  -48
+.Lrsaz_amm52x20_x2_256_epilogue:
+    ret
+.cfi_endproc
+.size   RSAZ_amm52x20_x2_256, .-RSAZ_amm52x20_x2_256
+___
+}
+
+###############################################################################
+# Constant time extraction from the precomputed table of powers base^i, where
+#    i = 0..2^EXP_WIN_SIZE-1
+#
+# The input |red_table| contains precomputations for two independent base values,
+# so the |tbl_idx| indicates for which base shall we extract the value.
+# |red_table_idx| is a power index.
+#
+# Extracted value (output) is 20 digit number in 2^52 radix.
+#
+# void extract_multiplier_2x20_win5(BN_ULONG *red_Y,
+#                                   const BN_ULONG red_table[1 << EXP_WIN_SIZE][2][20],
+#                                   int red_table_idx,
+#                                   int tbl_idx);           # 0 or 1
+#
+# EXP_WIN_SIZE = 5
+###############################################################################
+{
+# input parameters
+my ($out,$red_tbl,$red_tbl_idx,$tbl_idx) = @_6_args_universal_ABI;
+
+my ($t0,$t1,$t2,$t3,$t4) = map("%ymm$_", (0..4));
+my $t4xmm = $t4 =~ s/%y/%x/r;
+my ($tmp0,$tmp1,$tmp2,$tmp3,$tmp4) = map("%ymm$_", (16..20));
+my ($cur_idx,$idx,$ones) = map("%ymm$_", (21..23));
+
+$code.=<<___;
+.text
+
+.align 32
+.globl  extract_multiplier_2x20_win5
+.type   extract_multiplier_2x20_win5,\@function,4
+extract_multiplier_2x20_win5:
+.cfi_startproc
+    endbranch
+    leaq    ($tbl_idx,$tbl_idx,4), %rax
+    salq    \$5, %rax
+    addq    %rax, $red_tbl
+
+    vmovdqa64   .Lones(%rip), $ones         # broadcast ones
+    vpbroadcastq    $red_tbl_idx, $idx
+    leaq   `(1<<5)*2*20*8`($red_tbl), %rax  # holds end of the tbl
+
+    vpxor   $t4xmm, $t4xmm, $t4xmm
+    vmovdqa64   $t4, $t3                    # zeroing t0..4, cur_idx
+    vmovdqa64   $t4, $t2
+    vmovdqa64   $t4, $t1
+    vmovdqa64   $t4, $t0
+    vmovdqa64   $t4, $cur_idx
+
+.align 32
+.Lloop:
+    vpcmpq  \$0, $cur_idx, $idx, %k1        # mask of (idx == cur_idx)
+    addq    \$320, $red_tbl                 # 320 = 2 * 20 digits * 8 bytes
+    vpaddq  $ones, $cur_idx, $cur_idx       # increment cur_idx
+    vmovdqu64  -320($red_tbl), $tmp0        # load data from red_tbl
+    vmovdqu64  -288($red_tbl), $tmp1
+    vmovdqu64  -256($red_tbl), $tmp2
+    vmovdqu64  -224($red_tbl), $tmp3
+    vmovdqu64  -192($red_tbl), $tmp4
+    vpblendmq  $tmp0, $t0, ${t0}{%k1}       # extract data when mask is not zero
+    vpblendmq  $tmp1, $t1, ${t1}{%k1}
+    vpblendmq  $tmp2, $t2, ${t2}{%k1}
+    vpblendmq  $tmp3, $t3, ${t3}{%k1}
+    vpblendmq  $tmp4, $t4, ${t4}{%k1}
+    cmpq    $red_tbl, %rax
+    jne .Lloop
+
+    vmovdqu64   $t0, ($out)                 # store t0..4
+    vmovdqu64   $t1, 32($out)
+    vmovdqu64   $t2, 64($out)
+    vmovdqu64   $t3, 96($out)
+    vmovdqu64   $t4, 128($out)
+
+    ret
+.cfi_endproc
+.size   extract_multiplier_2x20_win5, .-extract_multiplier_2x20_win5
+___
+$code.=<<___;
+.data
+.align 32
+.Lones:
+    .quad   1,1,1,1
+___
+}
+
+if ($win64) {
+$rec="%rcx";
+$frame="%rdx";
+$context="%r8";
+$disp="%r9";
+
+$code.=<<___
+.extern     __imp_RtlVirtualUnwind
+.type   rsaz_def_handler,\@abi-omnipotent
+.align  16
+rsaz_def_handler:
+    push    %rsi
+    push    %rdi
+    push    %rbx
+    push    %rbp
+    push    %r12
+    push    %r13
+    push    %r14
+    push    %r15
+    pushfq
+    sub     \$64,%rsp
+
+    mov     120($context),%rax # pull context->Rax
+    mov     248($context),%rbx # pull context->Rip
+
+    mov     8($disp),%rsi      # disp->ImageBase
+    mov     56($disp),%r11     # disp->HandlerData
+
+    mov     0(%r11),%r10d      # HandlerData[0]
+    lea     (%rsi,%r10),%r10   # prologue label
+    cmp     %r10,%rbx          # context->Rip<.Lprologue
+    jb  .Lcommon_seh_tail
+
+    mov     152($context),%rax # pull context->Rsp
+
+    mov     4(%r11),%r10d      # HandlerData[1]
+    lea     (%rsi,%r10),%r10   # epilogue label
+    cmp     %r10,%rbx          # context->Rip>=.Lepilogue
+    jae     .Lcommon_seh_tail
+
+    lea     48(%rax),%rax
+
+    mov     -8(%rax),%rbx
+    mov     -16(%rax),%rbp
+    mov     -24(%rax),%r12
+    mov     -32(%rax),%r13
+    mov     -40(%rax),%r14
+    mov     -48(%rax),%r15
+    mov     %rbx,144($context) # restore context->Rbx
+    mov     %rbp,160($context) # restore context->Rbp
+    mov     %r12,216($context) # restore context->R12
+    mov     %r13,224($context) # restore context->R13
+    mov     %r14,232($context) # restore context->R14
+    mov     %r15,240($context) # restore context->R14
+
+.Lcommon_seh_tail:
+    mov     8(%rax),%rdi
+    mov     16(%rax),%rsi
+    mov     %rax,152($context) # restore context->Rsp
+    mov     %rsi,168($context) # restore context->Rsi
+    mov     %rdi,176($context) # restore context->Rdi
+
+    mov     40($disp),%rdi     # disp->ContextRecord
+    mov     $context,%rsi      # context
+    mov     \$154,%ecx         # sizeof(CONTEXT)
+    .long   0xa548f3fc         # cld; rep movsq
+
+    mov     $disp,%rsi
+    xor     %rcx,%rcx          # arg1, UNW_FLAG_NHANDLER
+    mov     8(%rsi),%rdx       # arg2, disp->ImageBase
+    mov     0(%rsi),%r8        # arg3, disp->ControlPc
+    mov     16(%rsi),%r9       # arg4, disp->FunctionEntry
+    mov     40(%rsi),%r10      # disp->ContextRecord
+    lea     56(%rsi),%r11      # &disp->HandlerData
+    lea     24(%rsi),%r12      # &disp->EstablisherFrame
+    mov     %r10,32(%rsp)      # arg5
+    mov     %r11,40(%rsp)      # arg6
+    mov     %r12,48(%rsp)      # arg7
+    mov     %rcx,56(%rsp)      # arg8, (NULL)
+    call    *__imp_RtlVirtualUnwind(%rip)
+
+    mov     \$1,%eax           # ExceptionContinueSearch
+    add     \$64,%rsp
+    popfq
+    pop     %r15
+    pop     %r14
+    pop     %r13
+    pop     %r12
+    pop     %rbp
+    pop     %rbx
+    pop     %rdi
+    pop     %rsi
+    ret
+.size   rsaz_def_handler,.-rsaz_def_handler
+
+.section    .pdata
+.align  4
+    .rva    .LSEH_begin_RSAZ_amm52x20_x1_256
+    .rva    .LSEH_end_RSAZ_amm52x20_x1_256
+    .rva    .LSEH_info_RSAZ_amm52x20_x1_256
+
+    .rva    .LSEH_begin_extract_multiplier_2x20_win5
+    .rva    .LSEH_end_extract_multiplier_2x20_win5
+    .rva    .LSEH_info_extract_multiplier_2x20_win5
+
+    .rva    .LSEH_begin_RSAZ_amm52x20_x2_256
+    .rva    .LSEH_end_RSAZ_amm52x20_x2_256
+    .rva    .LSEH_info_RSAZ_amm52x20_x2_256
+
+.section    .xdata
+.align  8
+.LSEH_info_RSAZ_amm52x20_x1_256:
+    .byte   9,0,0,0
+    .rva    rsaz_def_handler
+    .rva    .Lrsaz_amm52x20_x1_256_body,.Lrsaz_amm52x20_x1_256_epilogue
+.LSEH_info_extract_multiplier_2x20_win5:
+    .byte   9,0,0,0
+    .rva    rsaz_def_handler
+    .rva    .LSEH_begin_extract_multiplier_2x20_win5,.LSEH_begin_extract_multiplier_2x20_win5
+.LSEH_info_RSAZ_amm52x20_x2_256:
+    .byte   9,0,0,0
+    .rva    rsaz_def_handler
+    .rva    .Lrsaz_amm52x20_x2_256_body,.Lrsaz_amm52x20_x2_256_epilogue
+___
+}
+}}} else {{{                # fallback for old assembler
+$code.=<<___;
+.text
+
+.globl  rsaz_avx512ifma_eligible
+.type   rsaz_avx512ifma_eligible,\@abi-omnipotent
+rsaz_avx512ifma_eligible:
+    xor     %eax,%eax
+    ret
+.size   rsaz_avx512ifma_eligible, .-rsaz_avx512ifma_eligible
+
+.globl  RSAZ_amm52x20_x1_256
+.globl  RSAZ_amm52x20_x2_256
+.globl  extract_multiplier_2x20_win5
+.type   RSAZ_amm52x20_x1_256,\@abi-omnipotent
+RSAZ_amm52x20_x1_256:
+RSAZ_amm52x20_x2_256:
+extract_multiplier_2x20_win5:
+    .byte   0x0f,0x0b    # ud2
+    ret
+.size   RSAZ_amm52x20_x1_256, .-RSAZ_amm52x20_x1_256
+___
+}}}
+
+$code =~ s/\`([^\`]*)\`/eval $1/gem;
+print $code;
+close STDOUT or die "error closing STDOUT: $!";
diff --git a/crypto/bn/bn_exp.c b/crypto/bn/bn_exp.c
index 1254415a43..4f6445434b 100644
--- a/crypto/bn/bn_exp.c
+++ b/crypto/bn/bn_exp.c
@@ -1390,3 +1390,85 @@ int BN_mod_exp_simple(BIGNUM *r, const BIGNUM *a, const BIGNUM *p,
     bn_check_top(r);
     return ret;
 }
+
+/*
+ * This is a variant of modular exponentiation optimization that does
+ * parallel 2-primes exponentiation using 256-bit (AVX512VL) AVX512_IFMA ISA
+ * in 52-bit binary redundant representation.
+ * If such instructions are not available, or input data size is not supported,
+ * it falls back to two BN_mod_exp_mont_consttime() calls.
+ */
+int BN_mod_exp_mont_consttime_x2(BIGNUM *rr1, const BIGNUM *a1, const BIGNUM *p1,
+                                 const BIGNUM *m1, BN_MONT_CTX *in_mont1,
+                                 BIGNUM *rr2, const BIGNUM *a2, const BIGNUM *p2,
+                                 const BIGNUM *m2, BN_MONT_CTX *in_mont2,
+                                 BN_CTX *ctx)
+{
+    int ret = 0;
+
+#ifdef RSAZ_ENABLED
+    BN_MONT_CTX *mont1 = NULL;
+    BN_MONT_CTX *mont2 = NULL;
+
+    if (rsaz_avx512ifma_eligible() &&
+        ((a1->top == 16) && (p1->top == 16) && (BN_num_bits(m1) == 1024) &&
+         (a2->top == 16) && (p2->top == 16) && (BN_num_bits(m2) == 1024))) {
+
+        if (bn_wexpand(rr1, 16) == NULL)
+            goto err;
+        if (bn_wexpand(rr2, 16) == NULL)
+            goto err;
+
+        /*  Ensure that montgomery contexts are initialized */
+        if (in_mont1 != NULL) {
+            mont1 = in_mont1;
+        } else {
+            if ((mont1 = BN_MONT_CTX_new()) == NULL)
+                goto err;
+            if (!BN_MONT_CTX_set(mont1, m1, ctx))
+                goto err;
+        }
+        if (in_mont2 != NULL) {
+            mont2 = in_mont2;
+        } else {
+            if ((mont2 = BN_MONT_CTX_new()) == NULL)
+                goto err;
+            if (!BN_MONT_CTX_set(mont2, m2, ctx))
+                goto err;
+        }
+
+        ret = RSAZ_mod_exp_avx512_x2(rr1->d, a1->d, p1->d, m1->d, mont1->RR.d,
+                                     mont1->n0[0],
+                                     rr2->d, a2->d, p2->d, m2->d, mont2->RR.d,
+                                     mont2->n0[0],
+                                     1024 /* factor bit size */);
+
+        rr1->top = 16;
+        rr1->neg = 0;
+        bn_correct_top(rr1);
+        bn_check_top(rr1);
+
+        rr2->top = 16;
+        rr2->neg = 0;
+        bn_correct_top(rr2);
+        bn_check_top(rr2);
+
+        goto err;
+    }
+#endif
+
+    /* rr1 = a1^p1 mod m1 */
+    ret = BN_mod_exp_mont_consttime(rr1, a1, p1, m1, ctx, in_mont1);
+    /* rr2 = a2^p2 mod m2 */
+    ret &= BN_mod_exp_mont_consttime(rr2, a2, p2, m2, ctx, in_mont2);
+
+#ifdef RSAZ_ENABLED
+err:
+    if (in_mont2 == NULL)
+        BN_MONT_CTX_free(mont2);
+    if (in_mont1 == NULL)
+        BN_MONT_CTX_free(mont1);
+#endif
+
+    return ret;
+}
diff --git a/crypto/bn/build.info b/crypto/bn/build.info
index f732be24f8..237d5e90ed 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-avx2.s rsaz_exp_x2.c rsaz-avx512.s
   IF[{- $config{target} !~ /^VC/ -}]
     $BNASM_x86_64=asm/x86_64-gcc.c $BNASM_x86_64
   ELSE
@@ -154,6 +154,7 @@ GENERATE[x86_64-mont5.s]=asm/x86_64-mont5.pl
 GENERATE[x86_64-gf2m.s]=asm/x86_64-gf2m.pl
 GENERATE[rsaz-x86_64.s]=asm/rsaz-x86_64.pl
 GENERATE[rsaz-avx2.s]=asm/rsaz-avx2.pl
+GENERATE[rsaz-avx512.s]=asm/rsaz-avx512.pl
 
 GENERATE[bn-ia64.s]=asm/ia64.S
 GENERATE[ia64-mont.s]=asm/ia64-mont.pl
diff --git a/crypto/bn/rsaz_exp.h b/crypto/bn/rsaz_exp.h
index c05a5d937e..7d3a24b0d8 100644
--- a/crypto/bn/rsaz_exp.h
+++ b/crypto/bn/rsaz_exp.h
@@ -1,6 +1,6 @@
 /*
- * Copyright 2013-2018 The OpenSSL Project Authors. All Rights Reserved.
- * Copyright (c) 2012, Intel Corporation. All Rights Reserved.
+ * Copyright 2013-2020 The OpenSSL Project Authors. All Rights Reserved.
+ * Copyright (c) 2020, Intel Corporation. All Rights Reserved.
  *
  * Licensed under the Apache License 2.0 (the "License").  You may not use
  * this file except in compliance with the License.  You can obtain a copy
@@ -35,6 +35,23 @@ void RSAZ_512_mod_exp(BN_ULONG result[8],
                       const BN_ULONG m_norm[8], BN_ULONG k0,
                       const BN_ULONG RR[8]);
 
+
+int rsaz_avx512ifma_eligible(void);
+
+int RSAZ_mod_exp_avx512_x2(BN_ULONG *res1,
+                           const BN_ULONG *base1,
+                           const BN_ULONG *exponent1,
+                           const BN_ULONG *m1,
+                           const BN_ULONG *RR1,
+                           BN_ULONG k0_1,
+                           BN_ULONG *res2,
+                           const BN_ULONG *base2,
+                           const BN_ULONG *exponent2,
+                           const BN_ULONG *m2,
+                           const BN_ULONG *RR2,
+                           BN_ULONG k0_2,
+                           int factor_size);
+
 # endif
 
 #endif
diff --git a/crypto/bn/rsaz_exp_x2.c b/crypto/bn/rsaz_exp_x2.c
new file mode 100644
index 0000000000..f4c751bef0
--- /dev/null
+++ b/crypto/bn/rsaz_exp_x2.c
@@ -0,0 +1,542 @@
+/*
+ * Copyright 2020 The OpenSSL Project Authors. All Rights Reserved.
+ * Copyright (c) 2020, Intel Corporation. All Rights Reserved.
+ *
+ * Licensed under the Apache License 2.0 (the "License").  You may not use
+ * this file except in compliance with the License.  You can obtain a copy
+ * in the file LICENSE in the source distribution or at
+ * https://www.openssl.org/source/license.html
+ *
+ *
+ * Originally written by Ilya Albrekht, Sergey Kirillov and Andrey Matyukov
+ * Intel Corporation
+ *
+ */
+
+#include <openssl/opensslconf.h>
+#include "rsaz_exp.h"
+
+#ifndef RSAZ_ENABLED
+NON_EMPTY_TRANSLATION_UNIT
+#else
+# include <assert.h>
+# include <string.h>
+
+# if defined(__GNUC__)
+#  define ALIGN64 __attribute__((aligned(64)))
+# elif defined(_MSC_VER)
+#  define ALIGN64 __declspec(align(64))
+# else
+#  define ALIGN64
+# endif
+
+# define ALIGN_OF(ptr, boundary) \
+    ((unsigned char *)(ptr) + (boundary - (((size_t)(ptr)) & (boundary - 1))))
+
+/* Internal radix */
+# define DIGIT_SIZE (52)
+/* 52-bit mask */
+# define DIGIT_MASK ((uint64_t)0xFFFFFFFFFFFFF)
+
+# define BITS2WORD8_SIZE(x)  (((x) + 7) >> 3)
+# define BITS2WORD64_SIZE(x) (((x) + 63) >> 6)
+
+static ossl_inline uint64_t get_digit52(const uint8_t *in, int in_len);
+static ossl_inline void put_digit52(uint8_t *out, int out_len, uint64_t digit);
+static void to_words52(BN_ULONG *out, int out_len, const BN_ULONG *in,
+                       int in_bitsize);
+static void from_words52(BN_ULONG *bn_out, int out_bitsize, const BN_ULONG *in);
+static ossl_inline void set_bit(BN_ULONG *a, int idx);
+
+/* Number of |digit_size|-bit digits in |bitsize|-bit value */
+static ossl_inline int number_of_digits(int bitsize, int digit_size)
+{
+    return (bitsize + digit_size - 1) / digit_size;
+}
+
+typedef void (*AMM52)(BN_ULONG *res, const BN_ULONG *base,
+                      const BN_ULONG *exp, const BN_ULONG *m, BN_ULONG k0);
+typedef void (*EXP52_x2)(BN_ULONG *res, const BN_ULONG *base,
+                         const BN_ULONG *exp[2], const BN_ULONG *m,
+                         const BN_ULONG *rr, const BN_ULONG k0[2]);
+
+/*
+ * For details of the methods declared below please refer to
+ *    crypto/bn/asm/rsaz-avx512.pl
+ *
+ * Naming notes:
+ *  amm = Almost Montgomery Multiplication
+ *  ams = Almost Montgomery Squaring
+ *  52x20 - data represented as array of 20 digits in 52-bit radix
+ *  _x1_/_x2_ - 1 or 2 independent inputs/outputs
+ *  _256 suffix - uses 256-bit (AVX512VL) registers
+ */
+
+/*AMM = Almost Montgomery Multiplication. */
+void RSAZ_amm52x20_x1_256(BN_ULONG *res, const BN_ULONG *base,
+                          const BN_ULONG *exp, const BN_ULONG *m,
+                          BN_ULONG k0);
+void RSAZ_exp52x20_x2_256(BN_ULONG *res, const BN_ULONG *base,
+                      const BN_ULONG *exp[2], const BN_ULONG *m,
+                      const BN_ULONG *rr, const BN_ULONG k0[2]);
+void RSAZ_amm52x20_x2_256(BN_ULONG *out, const BN_ULONG *a,
+                          const BN_ULONG *b, const BN_ULONG *m,
+                          const BN_ULONG k0[2]);
+void extract_multiplier_2x20_win5(BN_ULONG *red_Y,
+                                  const BN_ULONG *red_table,
+                                  int red_table_idx, int tbl_idx);
+
+/*
+ * Dual Montgomery modular exponentiation using prime moduli of the
+ * same bit size, optimized with AVX512 ISA.
+ *
+ * Input and output parameters for each exponentiation are independent and
+ * denoted here by index |i|, i = 1..2.
+ *
+ * Input and output are all in regular 2^64 radix.
+ *
+ * Each moduli shall be |factor_size| bit size.
+ *
+ * NOTE: currently only 2x1024 case is supported.
+ *
+ *  [out] res|i|      - result of modular exponentiation: array of qword values
+ *                      in regular (2^64) radix. Size of array shall be enough
+ *                      to hold |factor_size| bits.
+ *  [in]  base|i|     - base
+ *  [in]  exp|i|      - exponent
+ *  [in]  m|i|        - moduli
+ *  [in]  rr|i|       - Montgomery parameter RR = R^2 mod m|i|
+ *  [in]  k0_|i|      - Montgomery parameter k0 = -1/m|i| mod 2^64
+ *  [in]  factor_size - moduli bit size
+ *
+ * \return 0 in case of failure,
+ *         1 in case of success.
+ */
+int RSAZ_mod_exp_avx512_x2(BN_ULONG *res1,
+                           const BN_ULONG *base1,
+                           const BN_ULONG *exp1,
+                           const BN_ULONG *m1,
+                           const BN_ULONG *rr1,
+                           BN_ULONG k0_1,
+                           BN_ULONG *res2,
+                           const BN_ULONG *base2,
+                           const BN_ULONG *exp2,
+                           const BN_ULONG *m2,
+                           const BN_ULONG *rr2,
+                           BN_ULONG k0_2,
+                           int factor_size)
+{
+    int ret = 0;
+
+    /*
+     * Number of word-size (BN_ULONG) digits to store exponent in redundant
+     * representation.
+     */
+    int exp_digits = number_of_digits(factor_size + 2, DIGIT_SIZE);
+    int coeff_pow = 4 * (DIGIT_SIZE * exp_digits - factor_size);
+    BN_ULONG *base1_red, *m1_red, *rr1_red;
+    BN_ULONG *base2_red, *m2_red, *rr2_red;
+    BN_ULONG *coeff_red;
+    BN_ULONG *storage = NULL;
+    BN_ULONG *storage_aligned = NULL;
+    BN_ULONG storage_len_bytes = 7 * exp_digits * sizeof(BN_ULONG);
+
+    /* AMM = Almost Montgomery Multiplication */
+    AMM52 amm = NULL;
+    /* Dual (2-exps in parallel) exponentiation */
+    EXP52_x2 exp_x2 = NULL;
+
+    const BN_ULONG *exp[2] = {0};
+    BN_ULONG k0[2] = {0};
+
+    /* Only 1024-bit factor size is supported now */
+    switch (factor_size) {
+    case 1024:
+        amm = RSAZ_amm52x20_x1_256;
+        exp_x2 = RSAZ_exp52x20_x2_256;
+        break;
+    default:
+        goto err;
+    }
+
+    storage = (BN_ULONG *)OPENSSL_malloc(storage_len_bytes + 64);
+    if (storage == NULL)
+        goto err;
+    storage_aligned = (BN_ULONG *)ALIGN_OF(storage, 64);
+
+    /* Memory layout for red(undant) representations */
+    base1_red = storage_aligned;
+    base2_red = storage_aligned + 1 * exp_digits;
+    m1_red    = storage_aligned + 2 * exp_digits;
+    m2_red    = storage_aligned + 3 * exp_digits;
+    rr1_red   = storage_aligned + 4 * exp_digits;
+    rr2_red   = storage_aligned + 5 * exp_digits;
+    coeff_red = storage_aligned + 6 * exp_digits;
+
+    /* Convert base_i, m_i, rr_i, from regular to 52-bit radix */
+    to_words52(base1_red, exp_digits, base1, factor_size);
+    to_words52(base2_red, exp_digits, base2, factor_size);
+    to_words52(m1_red, exp_digits, m1, factor_size);
+    to_words52(m2_red, exp_digits, m2, factor_size);
+    to_words52(rr1_red, exp_digits, rr1, factor_size);
+    to_words52(rr2_red, exp_digits, rr2, factor_size);
+
+    /*
+     * Compute target domain Montgomery converters RR' for each modulus
+     * based on precomputed original domain's RR.
+     *
+     * RR -> RR' transformation steps:
+     *  (1) coeff = 2^k
+     *  (2) t = AMM(RR,RR) = RR^2 / R' mod m
+     *  (3) RR' = AMM(t, coeff) = RR^2 * 2^k / R'^2 mod m
+     * where
+     *  k = 4 * (52 * digits52 - modlen)
+     *  R  = 2^(64 * ceil(modlen/64)) mod m
+     *  RR = R^2 mod M
+     *  R' = 2^(52 * ceil(modlen/52)) mod m
+     *
+     *  modlen = 1024: k = 64, RR = 2^2048 mod m, RR' = 2^2080 mod m
+     */
+    memset(coeff_red, 0, exp_digits * sizeof(BN_ULONG));
+    /* (1) in reduced domain representation */
+    set_bit(coeff_red, 64 * (int)(coeff_pow / 52) + coeff_pow % 52);
+
+    amm(rr1_red, rr1_red, rr1_red, m1_red, k0_1);     /* (2) for m1 */
+    amm(rr1_red, rr1_red, coeff_red, m1_red, k0_1);   /* (3) for m1 */
+
+    amm(rr2_red, rr2_red, rr2_red, m2_red, k0_2);     /* (2) for m2 */
+    amm(rr2_red, rr2_red, coeff_red, m2_red, k0_2);   /* (3) for m2 */
+
+    exp[0] = exp1;
+    exp[1] = exp2;
+
+    k0[0] = k0_1;
+    k0[1] = k0_2;
+
+    exp_x2(rr1_red, base1_red, exp, m1_red, rr1_red, k0);
+
+    /* Convert rr_i back to regular radix */
+    from_words52(res1, factor_size, rr1_red);
+    from_words52(res2, factor_size, rr2_red);
+
+    ret = 1;
+err:
+    if (storage != NULL) {
+        OPENSSL_cleanse(storage, storage_len_bytes);
+        OPENSSL_free(storage);
+    }
+    return ret;
+}
+
+/*
+ * Dual 1024-bit w-ary modular exponentiation using prime moduli of the same
+ * bit size using Almost Montgomery Multiplication, optimized with AVX512_IFMA
+ * ISA.
+ *
+ * The parameter w (window size) = 5.
+ *
+ *  [out] res      - result of modular exponentiation: 2x20 qword
+ *                   values in 2^52 radix.
+ *  [in]  base     - base (2x20 qword values in 2^52 radix)
+ *  [in]  exp      - array of 2 pointers to 16 qword values in 2^64 radix.
+ *                   Exponent is not converted to redundant representation.
+ *  [in]  m        - moduli (2x20 qword values in 2^52 radix)
+ *  [in]  rr       - Montgomery parameter for 2 moduli: RR = 2^2080 mod m.
+ *                   (2x20 qword values in 2^52 radix)
+ *  [in]  k0       - Montgomery parameter for 2 moduli: k0 = -1/m mod 2^64
+ *
+ * \return (void).
+ */
+void RSAZ_exp52x20_x2_256(BN_ULONG *out,          /* [2][20] */
+                          const BN_ULONG *base,   /* [2][20] */
+                          const BN_ULONG *exp[2], /* 2x16    */
+                          const BN_ULONG *m,      /* [2][20] */
+                          const BN_ULONG *rr,     /* [2][20] */
+                          const BN_ULONG k0[2])
+{
+# define BITSIZE_MODULUS (1024)
+# define EXP_WIN_SIZE (5)
+# define EXP_WIN_MASK ((1U << EXP_WIN_SIZE) - 1)
+/*
+ * Number of digits (64-bit words) in redundant representation to handle
+ * modulus bits
+ */
+# define RED_DIGITS (20)
+# define EXP_DIGITS (16)
+# define DAMM RSAZ_amm52x20_x2_256
+/*
+ * Squaring is done using multiplication now. That can be a subject of
+ * optimization in future.
+ */
+# define DAMS(r,a,m,k0) \
+              RSAZ_amm52x20_x2_256((r),(a),(a),(m),(k0))
+
+    /* Allocate stack for red(undant) result Y and multiplier X */
+    ALIGN64 BN_ULONG red_Y[2][RED_DIGITS];
+    ALIGN64 BN_ULONG red_X[2][RED_DIGITS];
+
+    /* Allocate expanded exponent */
+    ALIGN64 BN_ULONG expz[2][EXP_DIGITS + 1];
+
+    /* Pre-computed table of base powers */
+    ALIGN64 BN_ULONG red_table[1U << EXP_WIN_SIZE][2][RED_DIGITS];
+
+    int idx;
+
+    memset(red_Y, 0, sizeof(red_Y));
+    memset(red_table, 0, sizeof(red_table));
+    memset(red_X, 0, sizeof(red_X));
+
+    /*
+     * Compute table of powers base^i, i = 0, ..., (2^EXP_WIN_SIZE) - 1
+     *   table[0] = mont(x^0) = mont(1)
+     *   table[1] = mont(x^1) = mont(x)
+     */
+    red_X[0][0] = 1;
+    red_X[1][0] = 1;
+    DAMM(red_table[0][0], (const BN_ULONG*)red_X, rr, m, k0);
+    DAMM(red_table[1][0], base,  rr, m, k0);
+
+    for (idx = 1; idx < (int)((1U << EXP_WIN_SIZE) / 2); idx++) {
+        DAMS(red_table[2 * idx + 0][0], red_table[1 * idx][0], m, k0);
+        DAMM(red_table[2 * idx + 1][0], red_table[2 * idx][0], red_table[1][0], m, k0);
+    }
+
+    /* Copy and expand exponents */
+    memcpy(expz[0], exp[0], EXP_DIGITS * sizeof(BN_ULONG));
+    expz[0][EXP_DIGITS] = 0;
+    memcpy(expz[1], exp[1], EXP_DIGITS * sizeof(BN_ULONG));
+    expz[1][EXP_DIGITS] = 0;
+
+    /* Exponentiation */
+    {
+        int rem = BITSIZE_MODULUS % EXP_WIN_SIZE;
+        int delta = rem ? rem : EXP_WIN_SIZE;
+        BN_ULONG table_idx_mask = EXP_WIN_MASK;
+
+        int exp_bit_no = BITSIZE_MODULUS - delta;
+        int exp_chunk_no = exp_bit_no / 64;
+        int exp_chunk_shift = exp_bit_no % 64;
+
+        /* Process 1-st exp window - just init result */
+        BN_ULONG red_table_idx_0 = expz[0][exp_chunk_no];
+        BN_ULONG red_table_idx_1 = expz[1][exp_chunk_no];
+        /*
+         * The function operates with fixed moduli sizes divisible by 64,
+         * thus table index here is always in supported range [0, EXP_WIN_SIZE).
+         */
+        red_table_idx_0 >>= exp_chunk_shift;
+        red_table_idx_1 >>= exp_chunk_shift;
+
+        extract_multiplier_2x20_win5(red_Y[0], (const BN_ULONG*)red_table, (int)red_table_idx_0, 0);
+        extract_multiplier_2x20_win5(red_Y[1], (const BN_ULONG*)red_table, (int)red_table_idx_1, 1);
+
+        /* Process other exp windows */
+        for (exp_bit_no -= EXP_WIN_SIZE; exp_bit_no >= 0; exp_bit_no -= EXP_WIN_SIZE) {
+            /* Extract pre-computed multiplier from the table */
+            {
+                BN_ULONG T;
+
+                exp_chunk_no = exp_bit_no / 64;
+                exp_chunk_shift = exp_bit_no % 64;
+                {
+                    red_table_idx_0 = expz[0][exp_chunk_no];
+                    T = expz[0][exp_chunk_no + 1];
+
+                    red_table_idx_0 >>= exp_chunk_shift;
+                    /*
+                     * Get additional bits from then next quadword
+                     * when 64-bit boundaries are crossed.
+                     */
+                    if (exp_chunk_shift > 64 - EXP_WIN_SIZE) {
+                        T <<= (64 - exp_chunk_shift);
+                        red_table_idx_0 ^= T;
+                    }
+                    red_table_idx_0 &= table_idx_mask;
+
+                    extract_multiplier_2x20_win5(red_X[0], (const BN_ULONG*)red_table, (int)red_table_idx_0, 0);
+                }
+                {
+                    red_table_idx_1 = expz[1][exp_chunk_no];
+                    T = expz[1][exp_chunk_no + 1];
+
+                    red_table_idx_1 >>= exp_chunk_shift;
+                    /*
+                     * Get additional bits from then next quadword
+                     * when 64-bit boundaries are crossed.
+                     */
+                    if (exp_chunk_shift > 64 - EXP_WIN_SIZE) {
+                        T <<= (64 - exp_chunk_shift);
+                        red_table_idx_1 ^= T;
+                    }
+                    red_table_idx_1 &= table_idx_mask;
+
+                    extract_multiplier_2x20_win5(red_X[1], (const BN_ULONG*)red_table, (int)red_table_idx_1, 1);
+                }
+            }
+
+            /* Series of squaring */
+            DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
+            DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
+            DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
+            DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
+            DAMS((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, m, k0);
+
+            DAMM((BN_ULONG*)red_Y, (const BN_ULONG*)red_Y, (const BN_ULONG*)red_X, m, k0);
+        }
+    }
+
+    /*
+     *
+     * NB: After the last AMM of exponentiation in Montgomery domain, the result
+     * may be 1025-bit, but the conversion out of Montgomery domain performs an
+     * AMM(x,1) which guarantees that the final result is less than |m|, so no
+     * conditional subtraction is needed here. See "Efficient Software
+     * Implementations of Modular Exponentiation" (by Shay Gueron) paper for details.
+     */
+
+    /* Convert result back in regular 2^52 domain */
+    memset(red_X, 0, sizeof(red_X));
+    red_X[0][0] = 1;
+    red_X[1][0] = 1;
+    DAMM(out, (const BN_ULONG*)red_Y, (const BN_ULONG*)red_X, m, k0);
+
+    /* Clear exponents */
+    OPENSSL_cleanse(expz, sizeof(expz));
+    OPENSSL_cleanse(red_Y, sizeof(red_Y));
+
+# undef DAMS
+# undef DAMM
+# undef EXP_DIGITS
+# undef RED_DIGITS
+# undef EXP_WIN_MASK
+# undef EXP_WIN_SIZE
+# undef BITSIZE_MODULUS
+}
+
+static ossl_inline uint64_t get_digit52(const uint8_t *in, int in_len)
+{
+    uint64_t digit = 0;
+
+    assert(in != NULL);
+
+    for (; in_len > 0; in_len--) {
+        digit <<= 8;
+        digit += (uint64_t)(in[in_len - 1]);
+    }
+    return digit;
+}
+
+/*
+ * Convert array of words in regular (base=2^64) representation to array of
+ * words in redundant (base=2^52) one.
+ */
+static void to_words52(BN_ULONG *out, int out_len,
+                       const BN_ULONG *in, int in_bitsize)
+{
+    uint8_t *in_str = NULL;
+
+    assert(out != NULL);
+    assert(in != NULL);
+    /* Check destination buffer capacity */
+    assert(out_len >= number_of_digits(in_bitsize, DIGIT_SIZE));
+
+    in_str = (uint8_t *)in;
+
+    for (; in_bitsize >= (2 * DIGIT_SIZE); in_bitsize -= (2 * DIGIT_SIZE), out += 2) {
+        out[0] = (*(uint64_t *)in_str) & DIGIT_MASK;
+        in_str += 6;
+        out[1] = ((*(uint64_t *)in_str) >> 4) & DIGIT_MASK;
+        in_str += 7;
+        out_len -= 2;
+    }
+
+    if (in_bitsize > DIGIT_SIZE) {
+        uint64_t digit = get_digit52(in_str, 7);
+
+        out[0] = digit & DIGIT_MASK;
+        in_str += 6;
+        in_bitsize -= DIGIT_SIZE;
+        digit = get_digit52(in_str, BITS2WORD8_SIZE(in_bitsize));
+        out[1] = digit >> 4;
+        out += 2;
+        out_len -= 2;
+    } else if (in_bitsize > 0) {
+        out[0] = get_digit52(in_str, BITS2WORD8_SIZE(in_bitsize));
+        out++;
+        out_len--;
+    }
+
+    while (out_len > 0) {
+        *out = 0;
+        out_len--;
+        out++;
+    }
+}
+
+static ossl_inline void put_digit52(uint8_t *pStr, int strLen, uint64_t digit)
+{
+    assert(pStr != NULL);
+
+    for (; strLen > 0; strLen--) {
+        *pStr++ = (uint8_t)(digit & 0xFF);
+        digit >>= 8;
+    }
+}
+
+/*
+ * Convert array of words in redundant (base=2^52) representation to array of
+ * words in regular (base=2^64) one.
+ */
+static void from_words52(BN_ULONG *out, int out_bitsize, const BN_ULONG *in)
+{
+    int i;
+    int out_len = BITS2WORD64_SIZE(out_bitsize);
+
+    assert(out != NULL);
+    assert(in != NULL);
+
+    for (i = 0; i < out_len; i++)
+        out[i] = 0;
+
+    {
+        uint8_t *out_str = (uint8_t *)out;
+
+        for (; out_bitsize >= (2 * DIGIT_SIZE); out_bitsize -= (2 * DIGIT_SIZE), in += 2) {
+            (*(uint64_t *)out_str) = in[0];
+            out_str += 6;
+            (*(uint64_t *)out_str) ^= in[1] << 4;
+            out_str += 7;
+        }
+
+        if (out_bitsize > DIGIT_SIZE) {
+            put_digit52(out_str, 7, in[0]);
+            out_str += 6;
+            out_bitsize -= DIGIT_SIZE;
+            put_digit52(out_str, BITS2WORD8_SIZE(out_bitsize),
+                        (in[1] << 4 | in[0] >> 48));
+        } else if (out_bitsize) {
+            put_digit52(out_str, BITS2WORD8_SIZE(out_bitsize), in[0]);
+        }
+    }
+}
+
+/*
+ * Set bit at index |idx| in the words array |a|.
+ * It does not do any boundaries checks, make sure the index is valid before
+ * calling the function.
+ */
+static ossl_inline void set_bit(BN_ULONG *a, int idx)
+{
+    assert(a != NULL);
+
+    {
+        int i, j;
+
+        i = idx / BN_BITS2;
+        j = idx % BN_BITS2;
+        a[i] |= (((BN_ULONG)1) << j);
+    }
+}
+
+#endif
diff --git a/crypto/rsa/rsa_ossl.c b/crypto/rsa/rsa_ossl.c
index 9f98c037c8..1817392e76 100644
--- a/crypto/rsa/rsa_ossl.c
+++ b/crypto/rsa/rsa_ossl.c
@@ -688,15 +688,20 @@ static int rsa_ossl_mod_exp(BIGNUM *r0, const BIGNUM *I, RSA *rsa, BN_CTX *ctx)
         if (/* m1 = I moq q */
             !bn_from_mont_fixed_top(m1, I, rsa->_method_mod_q, ctx)
             || !bn_to_mont_fixed_top(m1, m1, rsa->_method_mod_q, ctx)
-            /* m1 = m1^dmq1 mod q */
-            || !BN_mod_exp_mont_consttime(m1, m1, rsa->dmq1, rsa->q, ctx,
-                                          rsa->_method_mod_q)
             /* r1 = I mod p */
             || !bn_from_mont_fixed_top(r1, I, rsa->_method_mod_p, ctx)
             || !bn_to_mont_fixed_top(r1, r1, rsa->_method_mod_p, ctx)
-            /* r1 = r1^dmp1 mod p */
-            || !BN_mod_exp_mont_consttime(r1, r1, rsa->dmp1, rsa->p, ctx,
-                                          rsa->_method_mod_p)
+            /*
+             * Use parallel exponentiations optimization if possible,
+             * otherwise fallback to two sequential exponentiations:
+             *    m1 = m1^dmq1 mod q
+             *    r1 = r1^dmp1 mod p
+             */
+            || !BN_mod_exp_mont_consttime_x2(m1, m1, rsa->dmq1, rsa->q,
+                                             rsa->_method_mod_q,
+                                             r1, r1, rsa->dmp1, rsa->p,
+                                             rsa->_method_mod_p,
+                                             ctx)
             /* r1 = (r1 - m1) mod p */
             /*
              * bn_mod_sub_fixed_top is not regular modular subtraction,
diff --git a/crypto/x86_64cpuid.pl b/crypto/x86_64cpuid.pl
index d3e2b9145a..9f7a1b092f 100644
--- a/crypto/x86_64cpuid.pl
+++ b/crypto/x86_64cpuid.pl
@@ -215,7 +215,7 @@ OPENSSL_ia32_cpuid:
 	cmp	\$0xe6,%eax
 	je	.Ldone
 	andl	\$0x3fdeffff,8(%rdi)	# ~(1<<31|1<<30|1<<21|1<<16)
-					# clear AVX512F+BW+VL+FIMA, all of
+					# clear AVX512F+BW+VL+IFMA, all of
 					# them are EVEX-encoded, which requires
 					# ZMM state support even if one uses
 					# only XMM and YMM :-(