ml_kem: add AVX2 NTT and inverse NTT

Implement the FIPS 203 transforms in x86_64 perlasm with AVX2 runtime dispatch.

Assisted-by: ChatGPT:gpt-6-Sol
Reviewed-by: Neil Horman <[email protected]>
Reviewed-by: Mounir Idrassi <[email protected]>
Merge-date: Wed Oct 21 16:38:25 2026
Merged-from: https://github.com/openssl/openssl/pull/33023

3 files changed, 497 insertions, 2 deletions

diff --git a/crypto/ml_kem/asm/mlkem_ntt-x86_64.pl b/crypto/ml_kem/asm/mlkem_ntt-x86_64.pl
new file mode 100644
index 00000000000..9c649bb48f3
--- /dev/null
+++ b/crypto/ml_kem/asm/mlkem_ntt-x86_64.pl
@@ -0,0 +1,467 @@
+#! /usr/bin/env perl
+# Copyright 2026 The OpenSSL Project Authors. All Rights Reserved.
+#
+# Licensed under the Apache License 2.0 (the "License").  You may not use
+# this file except in compliance with the License.  You can obtain a copy
+# in the file LICENSE in the source distribution or at
+# https://www.openssl.org/source/license.html
+
+# FIPS 203, Algorithms 9 and 10. Each 16-bit lane contains a canonical
+# coefficient.
+# Shoup multiplication uses floor(zeta * 2^16 / q) to avoid division:
+#   t = a*zeta - floor(a*zeta_shoup/2^16)*q, 0 <= t < 2*q.
+# Both products in the subtraction need only their low 16 bits. The forward
+# transform keeps bounded intermediate values and reduces its final output to
+# [0,q), as required by ml_kem.c.
+
+$output = $#ARGV >= 0 && $ARGV[$#ARGV] =~ m|\.\w+$| ? pop : undef;
+$flavour = $#ARGV >= 0 && $ARGV[0] !~ m|\.| ? shift : undef;
+$win64 = ($flavour =~ /[nm]asm|mingw64/ || $output =~ /\.asm$/);
+$arg = $win64 ? "%rcx" : "%rdi";
+
+$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";
+
+$avx2 = 0;
+if (`$ENV{CC} -Wa,-v -c -o /dev/null -x assembler /dev/null 2>&1`
+    =~ /GNU assembler version ([2-9]\.[0-9]+)/) {
+    $avx2 = ($1 >= 2.22);
+}
+if (!$avx2 && $win64 && ($flavour =~ /nasm/ || $ENV{ASM} =~ /nasm/)
+    && `nasm -v 2>&1` =~ /NASM version ([2-9]\.[0-9]+)(?:\.([0-9]+))?/) {
+    $avx2 = ($1 >= 2.10);
+}
+if (!$avx2 && `$ENV{CC} -v 2>&1`
+    =~ /((?:clang|LLVM) version|.*based on LLVM) ([0-9]+\.[0-9]+)/) {
+    $avx2 = ($2 >= 3.3);
+}
+
+open OUT, "| \"$^X\" \"$xlate\" $flavour \"$output\""
+  or die "can't call $xlate: $!";
+*STDOUT = *OUT;
+
+sub bitreverse7 {
+    my $v = shift;
+    my $r = 0;
+    for (1..7) {
+        $r = ($r << 1) | ($v & 1);
+        $v >>= 1;
+    }
+    return $r;
+}
+
+sub powmod {
+    my ($base, $exponent, $modulus) = @_;
+    my $result = 1;
+    while ($exponent) {
+        $result = ($result * $base) % $modulus if $exponent & 1;
+        $base = ($base * $base) % $modulus;
+        $exponent >>= 1;
+    }
+    return $result;
+}
+
+my @roots = map { powmod(17, bitreverse7($_), 3329) } 0..127;
+my @shoup = map { int($_ * 65536 / 3329) } @roots;
+my @invroots = (0, map { 3329 - $roots[128 - $_] } 1..127);
+my @invshoup = map { int($_ * 65536 / 3329) } @invroots;
+my $code = "";
+# Windows preserves XMM6..XMM15, so clear only the volatile vector registers.
+# Unix has no nonvolatile vector registers and can clear all of them at once.
+my $clear_vectors = $win64
+    ? (join('', map { "    vpxor %ymm$_, %ymm$_, %ymm$_\n" } 0..5)
+       . "    vzeroupper\n")
+    : "    vzeroall\n";
+
+if ($avx2) {
+    $code .= <<___;
+.text
+.extern OPENSSL_ia32cap_P
+.globl mlkem_ntt_avx2_capable
+.type mlkem_ntt_avx2_capable,\@abi-omnipotent
+.align 16
+mlkem_ntt_avx2_capable:
+    mov OPENSSL_ia32cap_P+8(%rip), %eax
+    and \$32, %eax
+    ret
+.Lmlkem_ntt_avx2_capable_end:
+.size mlkem_ntt_avx2_capable, .-mlkem_ntt_avx2_capable
+
+.globl mlkem_ntt_avx2
+.type mlkem_ntt_avx2,\@abi-omnipotent
+.align 32
+mlkem_ntt_avx2:
+    vpbroadcastw .Lmlkem_q(%rip), %ymm4
+___
+
+    # Multiply the peer by its twiddle factor, then do the butterfly. After
+    # stage n, each output is below (2*n+3)*q, hence below 15*q at the end.
+    # Shoup multiplication remains valid for every 16-bit input.
+    # Input: ymm0 = even, ymm1 = peer, ymm2 = zeta, ymm3 = Shoup factor.
+    # Output: ymm5 = new even, ymm0 = new peer. Only ymm0..ymm5 are used,
+    # so no vector register save is needed on either x86-64 ABI.
+    my $butterfly = <<'___';
+    vpmulhuw %ymm3, %ymm1, %ymm5
+    vpmullw %ymm2, %ymm1, %ymm1
+    vpmullw %ymm4, %ymm5, %ymm5
+    vpsubw %ymm5, %ymm1, %ymm1
+    vpaddw %ymm1, %ymm0, %ymm5
+    vpsubw %ymm1, %ymm0, %ymm0
+    vpaddw %ymm4, %ymm0, %ymm0
+    vpaddw %ymm4, %ymm0, %ymm0
+___
+
+    # Stages with distances 128, 64, 32, and 16. One vector holds
+    # 16 independent butterflies; each group has a single twiddle factor.
+    for my $stage (0..3) {
+        my $distance = 128 >> $stage;
+        my $base = 1 << $stage;
+        for my $group (0..($base - 1)) {
+            my $index = $base + $group;
+            for (my $pos = 0; $pos < $distance; $pos += 16) {
+                my $left = 2 * ($group * 2 * $distance + $pos);
+                my $right = $left + 2 * $distance;
+                $code .= "    vmovdqu $left($arg), %ymm0\n";
+                $code .= "    vmovdqu $right($arg), %ymm1\n";
+                $code .= "    vpbroadcastw .Lmlkem_root_$index(%rip), %ymm2\n";
+                $code .= "    vpbroadcastw .Lmlkem_root_$index+2(%rip), %ymm3\n";
+                $code .= $butterfly;
+                $code .= "    vmovdqu %ymm5, $left($arg)\n";
+                $code .= "    vmovdqu %ymm0, $right($arg)\n";
+            }
+        }
+    }
+
+    # For the final three stages, butterflies stay inside a 16-word
+    # vector. Shuffle the partners together, and then blend the results.
+    for my $distance (8, 4, 2) {
+        $code .= <<___;
+    lea .Lmlkem_zetas_$distance(%rip), %rdx
+    xor %eax, %eax
+.Lmlkem_stage_$distance:
+    vmovdqu ($arg,%rax), %ymm0
+___
+        if ($distance == 8) {
+            $code .= <<'___';
+    vperm2i128 $0x11, %ymm0, %ymm0, %ymm1
+    vperm2i128 $0x00, %ymm0, %ymm0, %ymm0
+___
+        } elsif ($distance == 4) {
+            $code .= <<'___';
+    vpshufd $0xee, %ymm0, %ymm1
+    vpshufd $0x44, %ymm0, %ymm0
+___
+        } else {
+            $code .= <<'___';
+    vpshufd $0xf5, %ymm0, %ymm1
+    vpshufd $0xa0, %ymm0, %ymm0
+___
+        }
+        $code .= <<___;
+    vmovdqu (%rdx,%rax), %ymm2
+    vmovdqu 512(%rdx,%rax), %ymm3
+___
+        $code .= $butterfly;
+        if ($distance == 8) {
+            $code .= <<'___';
+    vperm2i128 $0x20, %ymm0, %ymm5, %ymm1
+___
+        } else {
+            my $mask = $distance == 4 ? '0xcc' : '0xaa';
+            $code .= "    vpblendd \$$mask, %ymm0, %ymm5, %ymm1\n";
+        }
+        $code .= <<___;
+    vmovdqu %ymm1, ($arg,%rax)
+    add \$32, %rax
+    cmp \$512, %rax
+    jb .Lmlkem_stage_$distance
+___
+    }
+
+    $code .= <<___;
+    vpbroadcastw .Lmlkem_barrett(%rip), %ymm2
+    xor %eax, %eax
+.Lmlkem_forward_reduce:
+    vmovdqu ($arg,%rax), %ymm0
+    vpmulhuw %ymm2, %ymm0, %ymm1
+    vpmullw %ymm4, %ymm1, %ymm1
+    vpsubw %ymm1, %ymm0, %ymm0
+    vpsubw %ymm4, %ymm0, %ymm1
+    vpminuw %ymm1, %ymm0, %ymm0
+    vmovdqu %ymm0, ($arg,%rax)
+    add \$32, %rax
+    cmp \$512, %rax
+    jb .Lmlkem_forward_reduce
+    xor $arg, $arg
+___
+    $code .= $clear_vectors;
+    $code .= <<'___';
+    ret
+.Lmlkem_ntt_avx2_end:
+.size mlkem_ntt_avx2, .-mlkem_ntt_avx2
+___
+
+    # The inverse uses the negated forward roots in reverse order. This
+    # lets its peer butterfly multiply (even - odd), matching ml_kem.c.
+    # The sum doubles the input bound, while Shoup multiplication keeps the
+    # peer below 2*q. Unix keeps the bound in ymm6. Windows reads it from a
+    # constant table to avoid saving a nonvolatile vector register.
+    my $inverse_butterfly = <<'___';
+    vpaddw %ymm1, %ymm0, %ymm5
+    vpsubw %ymm1, %ymm0, %ymm0
+    vpaddw %ymm6, %ymm0, %ymm0
+    vpmulhuw %ymm3, %ymm0, %ymm1
+    vpmullw %ymm2, %ymm0, %ymm0
+    vpmullw %ymm4, %ymm1, %ymm1
+    vpsubw %ymm1, %ymm0, %ymm0
+___
+
+    $code .= <<___;
+.globl mlkem_inverse_ntt_avx2
+.type mlkem_inverse_ntt_avx2,\@abi-omnipotent
+.align 32
+mlkem_inverse_ntt_avx2:
+    vpbroadcastw .Lmlkem_q(%rip), %ymm4
+___
+
+    for my $distance (2, 4, 8) {
+        my $bound = $distance == 2 ? 'q'
+                  : $distance == 4 ? 'twice_q' : 'four_q';
+        $code .= "    vpbroadcastw .Lmlkem_$bound(%rip), %ymm6\n"
+            unless $win64;
+        $code .= <<___;
+    lea .Lmlkem_invzetas_$distance(%rip), %rdx
+    xor %eax, %eax
+.Lmlkem_invstage_$distance:
+    vmovdqu ($arg,%rax), %ymm0
+___
+        if ($distance == 8) {
+            $code .= <<'___';
+    vperm2i128 $0x11, %ymm0, %ymm0, %ymm1
+    vperm2i128 $0x00, %ymm0, %ymm0, %ymm0
+___
+        } elsif ($distance == 4) {
+            $code .= <<'___';
+    vpshufd $0xee, %ymm0, %ymm1
+    vpshufd $0x44, %ymm0, %ymm0
+___
+        } else {
+            $code .= <<'___';
+    vpshufd $0xf5, %ymm0, %ymm1
+    vpshufd $0xa0, %ymm0, %ymm0
+___
+        }
+        $code .= <<___;
+    vmovdqu (%rdx,%rax), %ymm2
+    vmovdqu 512(%rdx,%rax), %ymm3
+___
+        my $ops = $inverse_butterfly;
+        $ops =~ s/%ymm6/.Lmlkem_vec_$bound(%rip)/g if $win64;
+        $code .= $ops;
+        if ($distance == 8) {
+            $code .= <<'___';
+    vperm2i128 $0x20, %ymm0, %ymm5, %ymm1
+___
+        } else {
+            my $mask = $distance == 4 ? '0xcc' : '0xaa';
+            $code .= "    vpblendd \$$mask, %ymm0, %ymm5, %ymm1\n";
+        }
+        $code .= <<___;
+    vmovdqu %ymm1, ($arg,%rax)
+    add \$32, %rax
+    cmp \$512, %rax
+    jb .Lmlkem_invstage_$distance
+___
+    }
+
+    # After three stages every value is below 8*q. Restore canonical
+    # coefficients before the last four stages so their sums fit 16 bits.
+    $code .= <<___;
+    vpbroadcastw .Lmlkem_barrett(%rip), %ymm2
+    xor %eax, %eax
+.Lmlkem_inverse_reduce:
+    vmovdqu ($arg,%rax), %ymm0
+    vpmulhuw %ymm2, %ymm0, %ymm1
+    vpmullw %ymm4, %ymm1, %ymm1
+    vpsubw %ymm1, %ymm0, %ymm0
+    vpsubw %ymm4, %ymm0, %ymm1
+    vpminuw %ymm1, %ymm0, %ymm0
+    vmovdqu %ymm0, ($arg,%rax)
+    add \$32, %rax
+    cmp \$512, %rax
+    jb .Lmlkem_inverse_reduce
+___
+
+    for my $stage (0..3) {
+        my $distance = 16 << $stage;
+        my $base = 129 - 256 / $distance;
+        my $bound = ('q', 'twice_q', 'four_q', 'eight_q')[$stage];
+        $code .= "    vpbroadcastw .Lmlkem_$bound(%rip), %ymm6\n"
+            unless $win64;
+        for my $group (0..(128 / $distance - 1)) {
+            my $index = $base + $group;
+            for (my $pos = 0; $pos < $distance; $pos += 16) {
+                my $left = 2 * ($group * 2 * $distance + $pos);
+                my $right = $left + 2 * $distance;
+                $code .= "    vmovdqu $left($arg), %ymm0\n";
+                $code .= "    vmovdqu $right($arg), %ymm1\n";
+                $code .= "    vpbroadcastw .Lmlkem_invroot_$index(%rip), %ymm2\n";
+                $code .= "    vpbroadcastw .Lmlkem_invroot_$index+2(%rip), %ymm3\n";
+                my $ops = $inverse_butterfly;
+                $ops =~ s/%ymm6/.Lmlkem_vec_$bound(%rip)/g if $win64;
+                $code .= $ops;
+                $code .= "    vmovdqu %ymm5, $left($arg)\n";
+                $code .= "    vmovdqu %ymm0, $right($arg)\n";
+            }
+        }
+    }
+
+    # Multiply every output coefficient by 128^-1 mod q.
+    $code .= <<___;
+    vpbroadcastw .Lmlkem_inverse_degree(%rip), %ymm2
+    vpbroadcastw .Lmlkem_inverse_degree+2(%rip), %ymm3
+    xor %eax, %eax
+.Lmlkem_inverse_scale:
+    vmovdqu ($arg,%rax), %ymm0
+    vpmulhuw %ymm3, %ymm0, %ymm1
+    vpmullw %ymm2, %ymm0, %ymm0
+    vpmullw %ymm4, %ymm1, %ymm1
+    vpsubw %ymm1, %ymm0, %ymm0
+    vpsubw %ymm4, %ymm0, %ymm1
+    vpminuw %ymm1, %ymm0, %ymm0
+    vmovdqu %ymm0, ($arg,%rax)
+    add \$32, %rax
+    cmp \$512, %rax
+    jb .Lmlkem_inverse_scale
+    xor $arg, $arg
+___
+    $code .= $clear_vectors;
+    $code .= <<'___';
+    ret
+.Lmlkem_inverse_ntt_avx2_end:
+.size mlkem_inverse_ntt_avx2, .-mlkem_inverse_ntt_avx2
+
+.section .rodata
+.align 32
+.Lmlkem_q:
+    .long 3329
+.Lmlkem_twice_q:
+    .long 6658
+.Lmlkem_four_q:
+    .long 13316
+.Lmlkem_eight_q:
+    .long 26632
+.Lmlkem_barrett:
+    .long 19
+___
+    if ($win64) {
+        for my $bound (['q', 3329], ['twice_q', 6658],
+                       ['four_q', 13316], ['eight_q', 26632]) {
+            my $pair = $bound->[1] + ($bound->[1] << 16);
+            $code .= ".align 32\n.Lmlkem_vec_" . $bound->[0] . ":\n";
+            $code .= "    .long " . join(", ", ($pair) x 8) . "\n";
+        }
+    }
+    $code .= ".Lmlkem_inverse_degree:\n";
+    $code .= "    .long " . (3303 + (int(3303 * 65536 / 3329) << 16)) . "\n";
+    for my $index (1..15) {
+        $code .= ".Lmlkem_root_$index:\n";
+        $code .= "    .long " . ($roots[$index] + ($shoup[$index] << 16)) . "\n";
+    }
+    for my $index (113..127) {
+        $code .= ".Lmlkem_invroot_$index:\n";
+        $code .= "    .long " . ($invroots[$index]
+            + ($invshoup[$index] << 16)) . "\n";
+    }
+    for my $distance (8, 4, 2) {
+        my $base = 128 / $distance;
+        $code .= ".align 32\n.Lmlkem_zetas_$distance:\n";
+        for my $factor (0, 1) {
+            for my $chunk (0..15) {
+                my @words;
+                for my $lane (0..15) {
+                    my $index = $base + $chunk * (16 / (2 * $distance))
+                        + int($lane / (2 * $distance));
+                    push @words, $factor ? $shoup[$index] : $roots[$index];
+                }
+                my @packed = map { $words[2 * $_]
+                    + ($words[2 * $_ + 1] << 16) } 0..7;
+                $code .= "    .long " . join(", ", @packed) . "\n";
+            }
+        }
+    }
+    for my $distance (2, 4, 8) {
+        my $base = $distance == 2 ? 1 : $distance == 4 ? 65 : 97;
+        $code .= ".align 32\n.Lmlkem_invzetas_$distance:\n";
+        for my $factor (0, 1) {
+            for my $chunk (0..15) {
+                my @words;
+                for my $lane (0..15) {
+                    my $index = $base + $chunk * (16 / (2 * $distance))
+                        + int($lane / (2 * $distance));
+                    push @words, $factor ? $invshoup[$index]
+                        : $invroots[$index];
+                }
+                my @packed = map { $words[2 * $_]
+                    + ($words[2 * $_ + 1] << 16) } 0..7;
+                $code .= "    .long " . join(", ", @packed) . "\n";
+            }
+        }
+    }
+} else {
+    $code .= <<'___';
+.text
+.globl mlkem_ntt_avx2_capable
+.type mlkem_ntt_avx2_capable,@abi-omnipotent
+mlkem_ntt_avx2_capable:
+    xor %eax, %eax
+    ret
+.Lmlkem_ntt_avx2_capable_end:
+.size mlkem_ntt_avx2_capable, .-mlkem_ntt_avx2_capable
+.globl mlkem_ntt_avx2
+.type mlkem_ntt_avx2,@abi-omnipotent
+mlkem_ntt_avx2:
+    .byte 0x0f,0x0b
+.Lmlkem_ntt_avx2_end:
+.size mlkem_ntt_avx2, .-mlkem_ntt_avx2
+.globl mlkem_inverse_ntt_avx2
+.type mlkem_inverse_ntt_avx2,@abi-omnipotent
+mlkem_inverse_ntt_avx2:
+    .byte 0x0f,0x0b
+.Lmlkem_inverse_ntt_avx2_end:
+.size mlkem_inverse_ntt_avx2, .-mlkem_inverse_ntt_avx2
+___
+}
+
+if ($win64) {
+    # All three routines are leaf functions. They do not change RSP or any
+    # nonvolatile register, so a zero-code UNWIND_INFO describes each one.
+    $code .= <<'___';
+.section .pdata
+.align 4
+    .rva mlkem_ntt_avx2_capable
+    .rva .Lmlkem_ntt_avx2_capable_end
+    .rva .Lmlkem_ntt_avx2_capable_unwind
+    .rva mlkem_ntt_avx2
+    .rva .Lmlkem_ntt_avx2_end
+    .rva .Lmlkem_ntt_avx2_unwind
+    .rva mlkem_inverse_ntt_avx2
+    .rva .Lmlkem_inverse_ntt_avx2_end
+    .rva .Lmlkem_inverse_ntt_avx2_unwind
+
+.section .xdata
+.align 4
+.Lmlkem_ntt_avx2_capable_unwind:
+    .byte 1,0,0,0
+.Lmlkem_ntt_avx2_unwind:
+    .byte 1,0,0,0
+.Lmlkem_inverse_ntt_avx2_unwind:
+    .byte 1,0,0,0
+___
+}
+
+print $code;
+close STDOUT or die "error closing STDOUT: $!";
diff --git a/crypto/ml_kem/build.info b/crypto/ml_kem/build.info
index a1bd38568ca..0b6ab49e839 100644
--- a/crypto/ml_kem/build.info
+++ b/crypto/ml_kem/build.info
@@ -2,6 +2,8 @@ LIBS = ../../libcrypto
 
 $MLKEMASM=
 IF[{- !$disabled{asm} -}]
+  $MLKEMASM_x86_64=mlkem_ntt-x86_64.s
+  $MLKEMDEF_x86_64=MLKEM_NTT_X86_64_ASM
   $MLKEMDEF_ppc64=MLKEM_NTT_PPC_ASM
 
   IF[{- $target{sys_id} ne "AIX" && $target{sys_id} ne "MACOSX" -}]
@@ -36,3 +38,5 @@ ENDIF
 # need to be applied to all affected libraries and modules.
 DEFINE[../../libcrypto]=$MLKEMDEF
 DEFINE[../../providers/libfips.a]=$MLKEMDEF
+
+GENERATE[mlkem_ntt-x86_64.s]=asm/mlkem_ntt-x86_64.pl
diff --git a/crypto/ml_kem/ml_kem.c b/crypto/ml_kem/ml_kem.c
index aff9798ecd4..ba10e11615f 100644
--- a/crypto/ml_kem/ml_kem.c
+++ b/crypto/ml_kem/ml_kem.c
@@ -448,7 +448,7 @@ static CRYPTO_ONCE ml_kem_ntt_once = CRYPTO_ONCE_STATIC_INIT;
  *  scalar_inverse_ntt_demontgomerize
  *      Inverse NTT whose input may still carry an inverse-Montgomery factor
  *      R^-1 from a preceding inner_product_montgomery call.  On the generic
- *      and PPC paths the two pointers are identical (both point at the same
+ *      PPC, and AVX2 paths the two pointers are identical (both point at the same
  *      fully-reduced implementation), because the generic inner_product
  *      already produces fully-reduced outputs via Barrett reduction.  On the
  *      s390x/vec128 path the two are distinct:
@@ -482,12 +482,28 @@ static ml_kem_scalar_ntt_fn scalar_ntt = ossl_ml_kem_scalar_ntt_generic;
 static ml_kem_scalar_inverse_ntt_fn scalar_inverse_ntt = ossl_ml_kem_scalar_inverse_ntt_generic;
 /*
  * scalar_inverse_ntt_demontgomerize: used after inner_product_montgomery.  On
- * generic/PPC this is the same function as scalar_inverse_ntt.  On s390x it is
+ * generic/PPC/AVX2 this is the same function as scalar_inverse_ntt.  On s390x it is
  * a specialised variant that also removes the inverse-Montgomery factor R^-1
  * left by inner_product_montgomery_vec128.
  */
 static ml_kem_scalar_inverse_ntt_demontgomerize_fn scalar_inverse_ntt_demontgomerize = ossl_ml_kem_scalar_inverse_ntt_generic;
 
+#if defined(MLKEM_NTT_X86_64_ASM)
+int mlkem_ntt_avx2_capable(void);
+void mlkem_ntt_avx2(uint16_t *c);
+void mlkem_inverse_ntt_avx2(uint16_t *c);
+
+static void scalar_ntt_avx2(scalar *s)
+{
+    mlkem_ntt_avx2(s->c);
+}
+
+static void scalar_inverse_ntt_avx2(scalar *s)
+{
+    mlkem_inverse_ntt_avx2(s->c);
+}
+#endif
+
 #if defined(MLKEM_NTT_PPC_ASM) && defined(_ARCH_PPC64)
 /*
  * PPC64LE Platform supports.
@@ -569,6 +585,14 @@ static ml_kem_matrix_mult_intt_fn matrix_mult_intt = matrix_mult_intt_generic;
 
 static void ml_kem_ntt_init(void)
 {
+#if defined(MLKEM_NTT_X86_64_ASM)
+    if (mlkem_ntt_avx2_capable()) {
+        scalar_ntt = scalar_ntt_avx2;
+        scalar_inverse_ntt = scalar_inverse_ntt_avx2;
+        scalar_inverse_ntt_demontgomerize = scalar_inverse_ntt_avx2;
+    }
+#endif
+
 /*
  * Initialize NTT function pointers to PPC64le implementations if available.
  * Scalar implementations are used by default.