Add ec/asm/ecp_nistz256-x86.pl module.

Reviewed-by: Emilia Käsper <emilia@openssl.org>

2 files changed, 1826 insertions, 0 deletions

diff --git a/crypto/ec/Makefile b/crypto/ec/Makefile
index 852183af59..319e0039db 100644
--- a/crypto/ec/Makefile
+++ b/crypto/ec/Makefile
@@ -48,6 +48,9 @@ lib:	$(LIBOBJ)
 	$(RANLIB) $(LIB) || echo Never mind.
 	@touch lib
 
+ecp_nistz256-x86.s:	asm/ecp_nistz256-x86.pl
+	$(PERL) asm/ecp_nistz256-x86.pl $(PERLASM_SCHEME) $(CFLAGS) $(PROCESSOR) > $@
+
 ecp_nistz256-x86_64.s: asm/ecp_nistz256-x86_64.pl
 	$(PERL) asm/ecp_nistz256-x86_64.pl $(PERLASM_SCHEME) > $@
 
diff --git a/crypto/ec/asm/ecp_nistz256-x86.pl b/crypto/ec/asm/ecp_nistz256-x86.pl
new file mode 100755
index 0000000000..0670c690c6
--- /dev/null
+++ b/crypto/ec/asm/ecp_nistz256-x86.pl
@@ -0,0 +1,1823 @@
+#!/usr/bin/env perl
+
+# ====================================================================
+# Written by Andy Polyakov <appro@openssl.org> for the OpenSSL
+# project. The module is, however, dual licensed under OpenSSL and
+# CRYPTOGAMS licenses depending on where you obtain it. For further
+# details see http://www.openssl.org/~appro/cryptogams/.
+# ====================================================================
+#
+# ECP_NISTZ256 module for x86/SSE2.
+#
+# October 2014.
+#
+# Original ECP_NISTZ256 submission targeting x86_64 is detailed in
+# http://eprint.iacr.org/2013/816. In the process of adaptation
+# original .c module was made 32-bit savvy in order to make this
+# implementation possible.
+#
+#		with/without -DECP_NISTZ256_ASM
+# Pentium	+66-163%
+# PIII		+72-172%
+# P4		+65-132%
+# Core2		+90-215%
+# Sandy Bridge	+105-265% (contemporary i[57]-* are all close to this)
+# Atom		+65-155%
+# Opteron	+54-110%
+# Bulldozer	+99-240%
+# VIA Nano	+93-290%
+#
+# Ranges denote minimum and maximum improvement coefficients depending
+# on benchmark. Lower coefficients are for ECDSA sign, server-side
+# operation. Keep in mind that +200% means 3x improvement.
+
+$0 =~ m/(.*[\/\\])[^\/\\]+$/; $dir=$1;
+push(@INC,"${dir}","${dir}../../perlasm");
+require "x86asm.pl";
+
+&asm_init($ARGV[0],"ecp_nistz256-x86.pl",$ARGV[$#ARGV] eq "386");
+
+$sse2=0;
+for (@ARGV) { $sse2=1 if (/-DOPENSSL_IA32_SSE2/); }
+
+&external_label("OPENSSL_ia32cap_P") if ($sse2);
+
+
+########################################################################
+# Convert ecp_nistz256_table.c to layout expected by ecp_nistz_gather_w7
+#
+open TABLE,"<ecp_nistz256_table.c"		or
+open TABLE,"<${dir}../ecp_nistz256_table.c"	or
+die "failed to open ecp_nistz256_table.c:",$!;
+
+use integer;
+
+foreach(<TABLE>) {
+	s/TOBN\(\s*(0x[0-9a-f]+),\s*(0x[0-9a-f]+)\s*\)/push @arr,hex($2),hex($1)/geo;
+}
+close TABLE;
+
+# See ecp_nistz256_table.c for explanation for why it's 64*16*37.
+# 64*16*37-1 is because $#arr returns last valid index or @arr, not
+# amount of elements.
+die "insane number of elements" if ($#arr != 64*16*37-1);
+
+&public_label("ecp_nistz256_precomputed");
+&align(4096);
+&set_label("ecp_nistz256_precomputed");
+
+########################################################################
+# this conversion smashes P256_POINT_AFFINE by individual bytes with
+# 64 byte interval, similar to
+#	1111222233334444
+#	1234123412341234
+for(1..37) {
+	@tbl = splice(@arr,0,64*16);
+	for($i=0;$i<64;$i++) {
+		undef @line;
+		for($j=0;$j<64;$j++) {
+			push @line,(@tbl[$j*16+$i/4]>>(($i%4)*8))&0xff;
+		}
+		&data_byte(join(',',map { sprintf "0x%02x",$_} @line));
+	}
+}
+
+########################################################################
+# Keep in mind that constants are stored least to most significant word
+&static_label("RR");
+&set_label("RR",64);
+&data_word(3,0,-1,-5,-2,-1,-3,4);	# 2^512 mod P-256
+
+&static_label("ONE_mont");
+&set_label("ONE_mont");
+&data_word(1,0,0,-1,-1,-1,-2,0);
+
+&static_label("ONE");
+&set_label("ONE");
+&data_word(1,0,0,0,0,0,0,0);
+&asciz("ECP_NISZ256 for x86/SSE2, CRYPTOGAMS by <appro\@openssl.org>");
+&align(64);
+
+########################################################################
+# void ecp_nistz256_mul_by_2(BN_ULONG edi[8],const BN_ULONG esi[8]);
+&function_begin("ecp_nistz256_mul_by_2");
+	&mov	("esi",&wparam(1));
+	&mov	("edi",&wparam(0));
+	&mov	("ebp","esi");
+########################################################################
+# common pattern for internal functions is that %edi is result pointer,
+# %esi and %ebp are input ones, %ebp being optional. %edi is preserved.
+	&call	("_ecp_nistz256_add");
+&function_end("ecp_nistz256_mul_by_2");
+
+########################################################################
+# void ecp_nistz256_mul_by_3(BN_ULONG edi[8],const BN_ULONG esi[8]);
+&function_begin("ecp_nistz256_mul_by_3");
+	&mov	("esi",&wparam(1));
+					# multiplication by 3 is performed
+					# as 2*n+n, but we can't use output
+					# to store 2*n, because if output
+					# pointer equals to input, then
+					# we'll get 2*n+2*n.
+	&stack_push(8);			# therefore we need to allocate
+					# 256-bit intermediate buffer.
+	&mov	("edi","esp");
+	&mov	("ebp","esi");
+	&call	("_ecp_nistz256_add");
+	&lea	("esi",&DWP(0,"edi"));
+	&mov	("ebp",&wparam(1));
+	&mov	("edi",&wparam(0));
+	&call	("_ecp_nistz256_add");
+	&stack_pop(8);
+&function_end("ecp_nistz256_mul_by_3");
+
+########################################################################
+# void ecp_nistz256_div_by_2(BN_ULONG edi[8],const BN_ULONG esi[8]);
+&function_begin("ecp_nistz256_div_by_2");
+	&mov	("esi",&wparam(1));
+	&mov	("edi",&wparam(0));
+	&call	("_ecp_nistz256_div_by_2");
+&function_end("ecp_nistz256_div_by_2");
+
+&function_begin_B("_ecp_nistz256_div_by_2");
+	# tmp = a is odd ? a+mod : a
+	#
+	# note that because mod has special form, i.e. consists of
+	# 0xffffffff, 1 and 0s, we can conditionally synthesize it by
+	# assigning least significant bit of input to one register,
+	# %ebp, and its negative to another, %edx.
+
+	&mov	("ebp",&DWP(0,"esi"));
+	&xor	("edx","edx");
+	&mov	("ebx",&DWP(4,"esi"));
+	&mov	("eax","ebp");
+	&and	("ebp",1);
+	&mov	("ecx",&DWP(8,"esi"));
+	&sub	("edx","ebp");
+
+	&add	("eax","edx");
+	&adc	("ebx","edx");
+	&mov	(&DWP(0,"edi"),"eax");
+	&adc	("ecx","edx");
+	&mov	(&DWP(4,"edi"),"ebx");
+	&mov	(&DWP(8,"edi"),"ecx");
+
+	&mov	("eax",&DWP(12,"esi"));
+	&mov	("ebx",&DWP(16,"esi"));
+	&adc	("eax",0);
+	&mov	("ecx",&DWP(20,"esi"));
+	&adc	("ebx",0);
+	&mov	(&DWP(12,"edi"),"eax");
+	&adc	("ecx",0);
+	&mov	(&DWP(16,"edi"),"ebx");
+	&mov	(&DWP(20,"edi"),"ecx");
+
+	&mov	("eax",&DWP(24,"esi"));
+	&mov	("ebx",&DWP(28,"esi"));
+	&adc	("eax","ebp");
+	&adc	("ebx","edx");
+	&mov	(&DWP(24,"edi"),"eax");
+	&sbb	("esi","esi");			# broadcast carry bit
+	&mov	(&DWP(28,"edi"),"ebx");
+
+	# ret = tmp >> 1
+
+	&mov	("eax",&DWP(0,"edi"));
+	&mov	("ebx",&DWP(4,"edi"));
+	&mov	("ecx",&DWP(8,"edi"));
+	&mov	("edx",&DWP(12,"edi"));
+
+	&shr	("eax",1);
+	&mov	("ebp","ebx");
+	&shl	("ebx",31);
+	&or	("eax","ebx");
+
+	&shr	("ebp",1);
+	&mov	("ebx","ecx");
+	&shl	("ecx",31);
+	&mov	(&DWP(0,"edi"),"eax");
+	&or	("ebp","ecx");
+	&mov	("eax",&DWP(16,"edi"));
+
+	&shr	("ebx",1);
+	&mov	("ecx","edx");
+	&shl	("edx",31);
+	&mov	(&DWP(4,"edi"),"ebp");
+	&or	("ebx","edx");
+	&mov	("ebp",&DWP(20,"edi"));
+
+	&shr	("ecx",1);
+	&mov	("edx","eax");
+	&shl	("eax",31);
+	&mov	(&DWP(8,"edi"),"ebx");
+	&or	("ecx","eax");
+	&mov	("ebx",&DWP(24,"edi"));
+
+	&shr	("edx",1);
+	&mov	("eax","ebp");
+	&shl	("ebp",31);
+	&mov	(&DWP(12,"edi"),"ecx");
+	&or	("edx","ebp");
+	&mov	("ecx",&DWP(28,"edi"));
+
+	&shr	("eax",1);
+	&mov	("ebp","ebx");
+	&shl	("ebx",31);
+	&mov	(&DWP(16,"edi"),"edx");
+	&or	("eax","ebx");
+
+	&shr	("ebp",1);
+	&mov	("ebx","ecx");
+	&shl	("ecx",31);
+	&mov	(&DWP(20,"edi"),"eax");
+	&or	("ebp","ecx");
+
+	&shr	("ebx",1);
+	&shl	("esi",31);
+	&mov	(&DWP(24,"edi"),"ebp");
+	&or	("ebx","esi");			# handle top-most carry bit
+	&mov	(&DWP(28,"edi"),"ebx");
+
+	&ret	();
+&function_end_B("_ecp_nistz256_div_by_2");
+
+########################################################################
+# void ecp_nistz256_add(BN_ULONG edi[8],const BN_ULONG esi[8],
+#					const BN_ULONG ebp[8]);
+&function_begin("ecp_nistz256_add");
+	&mov	("esi",&wparam(1));
+	&mov	("ebp",&wparam(2));
+	&mov	("edi",&wparam(0));
+	&call	("_ecp_nistz256_add");
+&function_end("ecp_nistz256_add");
+
+&function_begin_B("_ecp_nistz256_add");
+	&mov	("eax",&DWP(0,"esi"));
+	&mov	("ebx",&DWP(4,"esi"));
+	&mov	("ecx",&DWP(8,"esi"));
+	&add	("eax",&DWP(0,"ebp"));
+	&mov	("edx",&DWP(12,"esi"));
+	&adc	("ebx",&DWP(4,"ebp"));
+	&mov	(&DWP(0,"edi"),"eax");
+	&adc	("ecx",&DWP(8,"ebp"));
+	&mov	(&DWP(4,"edi"),"ebx");
+	&adc	("edx",&DWP(12,"ebp"));
+	&mov	(&DWP(8,"edi"),"ecx");
+	&mov	(&DWP(12,"edi"),"edx");
+
+	&mov	("eax",&DWP(16,"esi"));
+	&mov	("ebx",&DWP(20,"esi"));
+	&mov	("ecx",&DWP(24,"esi"));
+	&adc	("eax",&DWP(16,"ebp"));
+	&mov	("edx",&DWP(28,"esi"));
+	&adc	("ebx",&DWP(20,"ebp"));
+	&mov	(&DWP(16,"edi"),"eax");
+	&adc	("ecx",&DWP(24,"ebp"));
+	&mov	(&DWP(20,"edi"),"ebx");
+	&adc	("edx",&DWP(28,"ebp"));
+	&mov	(&DWP(24,"edi"),"ecx");
+	&sbb	("esi","esi");			# broadcast carry bit
+	&mov	(&DWP(28,"edi"),"edx");
+
+	# if a+b carries, subtract modulus.
+	#
+	# Note that because mod has special form, i.e. consists of
+	# 0xffffffff, 1 and 0s, we can conditionally synthesize it by
+	# assigning carry bit to one register, %ebp, and its negative
+	# to another, %esi. But we started by calculating %esi...
+
+	&mov	("eax",&DWP(0,"edi"));
+	&mov	("ebp","esi");
+	&mov	("ebx",&DWP(4,"edi"));
+	&shr	("ebp",31);
+	&mov	("ecx",&DWP(8,"edi"));
+	&sub	("eax","esi");
+	&mov	("edx",&DWP(12,"edi"));
+	&sbb	("ebx","esi");
+	&mov	(&DWP(0,"edi"),"eax");
+	&sbb	("ecx","esi");
+	&mov	(&DWP(4,"edi"),"ebx");
+	&sbb	("edx",0);
+	&mov	(&DWP(8,"edi"),"ecx");
+	&mov	(&DWP(12,"edi"),"edx");
+
+	&mov	("eax",&DWP(16,"edi"));
+	&mov	("ebx",&DWP(20,"edi"));
+	&mov	("ecx",&DWP(24,"edi"));
+	&sbb	("eax",0);
+	&mov	("edx",&DWP(28,"edi"));
+	&sbb	("ebx",0);
+	&mov	(&DWP(16,"edi"),"eax");
+	&sbb	("ecx","ebp");
+	&mov	(&DWP(20,"edi"),"ebx");
+	&sbb	("edx","esi");
+	&mov	(&DWP(24,"edi"),"ecx");
+	&mov	(&DWP(28,"edi"),"edx");
+
+	&ret	();
+&function_end_B("_ecp_nistz256_add");
+
+########################################################################
+# void ecp_nistz256_sub(BN_ULONG edi[8],const BN_ULONG esi[8],
+#					const BN_ULONG ebp[8]);
+&function_begin("ecp_nistz256_sub");
+	&mov	("esi",&wparam(1));
+	&mov	("ebp",&wparam(2));
+	&mov	("edi",&wparam(0));
+	&call	("_ecp_nistz256_sub");
+&function_end("ecp_nistz256_sub");
+
+&function_begin_B("_ecp_nistz256_sub");
+	&mov	("eax",&DWP(0,"esi"));
+	&mov	("ebx",&DWP(4,"esi"));
+	&mov	("ecx",&DWP(8,"esi"));
+	&sub	("eax",&DWP(0,"ebp"));
+	&mov	("edx",&DWP(12,"esi"));
+	&sbb	("ebx",&DWP(4,"ebp"));
+	&mov	(&DWP(0,"edi"),"eax");
+	&sbb	("ecx",&DWP(8,"ebp"));
+	&mov	(&DWP(4,"edi"),"ebx");
+	&sbb	("edx",&DWP(12,"ebp"));
+	&mov	(&DWP(8,"edi"),"ecx");
+	&mov	(&DWP(12,"edi"),"edx");
+
+	&mov	("eax",&DWP(16,"esi"));
+	&mov	("ebx",&DWP(20,"esi"));
+	&mov	("ecx",&DWP(24,"esi"));
+	&sbb	("eax",&DWP(16,"ebp"));
+	&mov	("edx",&DWP(28,"esi"));
+	&sbb	("ebx",&DWP(20,"ebp"));
+	&sbb	("ecx",&DWP(24,"ebp"));
+	&mov	(&DWP(16,"edi"),"eax");
+	&sbb	("edx",&DWP(28,"ebp"));
+	&mov	(&DWP(20,"edi"),"ebx");
+	&sbb	("esi","esi");			# broadcast borrow bit
+	&mov	(&DWP(24,"edi"),"ecx");
+	&mov	(&DWP(28,"edi"),"edx");
+
+	# if a-b borrows, add modulus.
+	#
+	# Note that because mod has special form, i.e. consists of
+	# 0xffffffff, 1 and 0s, we can conditionally synthesize it by
+	# assigning borrow bit to one register, %ebp, and its negative
+	# to another, %esi. But we started by calculating %esi...
+
+	&mov	("eax",&DWP(0,"edi"));
+	&mov	("ebp","esi");
+	&mov	("ebx",&DWP(4,"edi"));
+	&shr	("ebp",31);
+	&mov	("ecx",&DWP(8,"edi"));
+	&add	("eax","esi");
+	&mov	("edx",&DWP(12,"edi"));
+	&adc	("ebx","esi");
+	&mov	(&DWP(0,"edi"),"eax");
+	&adc	("ecx","esi");
+	&mov	(&DWP(4,"edi"),"ebx");
+	&adc	("edx",0);
+	&mov	(&DWP(8,"edi"),"ecx");
+	&mov	(&DWP(12,"edi"),"edx");
+
+	&mov	("eax",&DWP(16,"edi"));
+	&mov	("ebx",&DWP(20,"edi"));
+	&mov	("ecx",&DWP(24,"edi"));
+	&adc	("eax",0);
+	&mov	("edx",&DWP(28,"edi"));
+	&adc	("ebx",0);
+	&mov	(&DWP(16,"edi"),"eax");
+	&adc	("ecx","ebp");
+	&mov	(&DWP(20,"edi"),"ebx");
+	&adc	("edx","esi");
+	&mov	(&DWP(24,"edi"),"ecx");
+	&mov	(&DWP(28,"edi"),"edx");
+
+	&ret	();
+&function_end_B("_ecp_nistz256_sub");
+
+########################################################################
+# void ecp_nistz256_neg(BN_ULONG edi[8],const BN_ULONG esi[8]);
+&function_begin("ecp_nistz256_neg");
+	&mov	("ebp",&wparam(1));
+	&mov	("edi",&wparam(0));
+
+	&xor	("eax","eax");
+	&stack_push(8);
+	&mov	(&DWP(0,"esp"),"eax");
+	&mov	("esi","esp");
+	&mov	(&DWP(4,"esp"),"eax");
+	&mov	(&DWP(8,"esp"),"eax");
+	&mov	(&DWP(12,"esp"),"eax");
+	&mov	(&DWP(16,"esp"),"eax");
+	&mov	(&DWP(20,"esp"),"eax");
+	&mov	(&DWP(24,"esp"),"eax");
+	&mov	(&DWP(28,"esp"),"eax");
+	
+	&call	("_ecp_nistz256_sub");
+
+	&stack_pop(8);
+&function_end("ecp_nistz256_neg");
+
+&function_begin_B("_picup_eax");
+	&mov	("eax",&DWP(0,"esp"));
+	&ret	();
+&function_end_B("_picup_eax");
+
+########################################################################
+# void ecp_nistz256_to_mont(BN_ULONG edi[8],const BN_ULONG esi[8]);
+&function_begin("ecp_nistz256_to_mont");
+	&mov	("esi",&wparam(1));
+	&call	("_picup_eax");
+    &set_label("pic");
+	&lea	("ebp",&DWP(&label("RR")."-".&label("pic"),"eax"));
+						if ($sse2) {
+	&picmeup("eax","OPENSSL_ia32cap_P","eax",&label("pic"));
+	&mov	("eax",&DWP(0,"eax"));		}
+	&mov	("edi",&wparam(0));
+	&call	("_ecp_nistz256_mul_mont");
+&function_end("ecp_nistz256_to_mont");
+
+########################################################################
+# void ecp_nistz256_from_mont(BN_ULONG edi[8],const BN_ULONG esi[8]);
+&function_begin("ecp_nistz256_from_mont");
+	&mov	("esi",&wparam(1));
+	&call	("_picup_eax");
+    &set_label("pic");
+	&lea	("ebp",&DWP(&label("ONE")."-".&label("pic"),"eax"));
+						if ($sse2) {
+	&picmeup("eax","OPENSSL_ia32cap_P","eax",&label("pic"));
+	&mov	("eax",&DWP(0,"eax"));		}
+	&mov	("edi",&wparam(0));
+	&call	("_ecp_nistz256_mul_mont");
+&function_end("ecp_nistz256_from_mont");
+
+########################################################################
+# void ecp_nistz256_mul_mont(BN_ULONG edi[8],const BN_ULONG esi[8],
+#					     const BN_ULONG ebp[8]);
+&function_begin("ecp_nistz256_mul_mont");
+	&mov	("esi",&wparam(1));
+	&mov	("ebp",&wparam(2));
+						if ($sse2) {
+	&call	("_picup_eax");
+    &set_label("pic");
+	&picmeup("eax","OPENSSL_ia32cap_P","eax",&label("pic"));
+	&mov	("eax",&DWP(0,"eax"));		}
+	&mov	("edi",&wparam(0));
+	&call	("_ecp_nistz256_mul_mont");
+&function_end("ecp_nistz256_mul_mont");
+
+########################################################################
+# void ecp_nistz256_sqr_mont(BN_ULONG edi[8],const BN_ULONG esi[8]);
+&function_begin("ecp_nistz256_sqr_mont");
+	&mov	("esi",&wparam(1));
+						if ($sse2) {
+	&call	("_picup_eax");
+    &set_label("pic");
+	&picmeup("eax","OPENSSL_ia32cap_P","eax",&label("pic"));
+	&mov	("eax",&DWP(0,"eax"));		}
+	&mov	("edi",&wparam(0));
+	&mov	("ebp","esi");
+	&call	("_ecp_nistz256_mul_mont");
+&function_end("ecp_nistz256_sqr_mont");
+
+&function_begin_B("_ecp_nistz256_mul_mont");
+						if ($sse2) {
+	&and	("eax",1<<24|1<<26);
+	&cmp	("eax",1<<24|1<<26);		# see if XMM+SSE2 is on
+	&jne	(&label("mul_mont_ialu"));
+
+	########################################
+	# SSE2 code path featuring 32x16-bit
+	# multiplications is ~2x faster than
+	# IALU counterpart (except on Atom)...
+	########################################
+	# stack layout:
+	# +------------------------------------+< %esp
+	# | 7 16-byte temporary XMM words,     |
+	# | "sliding" toward lower address     |
+	# .                                    .
+	# +------------------------------------+
+	# | unused XMM word                    |
+	# +------------------------------------+< +128,%ebx
+	# | 8 16-byte XMM words holding copies |
+	# | of a[i]<<64|a[i]                   |
+	# .                                    .
+	# .                                    .
+	# +------------------------------------+< +256
+	&mov	("edx","esp");
+	&sub	("esp",0x100);
+
+	&movd	("xmm7",&DWP(0,"ebp"));		# b[0] -> 0000.00xy
+	&lea	("ebp",&DWP(4,"ebp"));
+	&pcmpeqd("xmm6","xmm6");
+	&psrlq	("xmm6",48);			# compose 0xffff<<64|0xffff
+
+	&pshuflw("xmm7","xmm7",0b11011100);	# 0000.00xy -> 0000.0x0y
+	&and	("esp",-64);
+	&pshufd	("xmm7","xmm7",0b11011100);	# 0000.0x0y -> 000x.000y
+	&lea	("ebx",&DWP(0x80,"esp"));
+
+	&movd	("xmm0",&DWP(4*0,"esi"));	# a[0] -> 0000.00xy
+	&pshufd	("xmm0","xmm0",0b11001100);	# 0000.00xy -> 00xy.00xy
+	&movd	("xmm1",&DWP(4*1,"esi"));	# a[1] -> ...
+	&movdqa	(&QWP(0x00,"ebx"),"xmm0");	# offload converted a[0]
+	&pmuludq("xmm0","xmm7");		# a[0]*b[0]
+
+	&movd	("xmm2",&DWP(4*2,"esi"));
+	&pshufd	("xmm1","xmm1",0b11001100);
+	&movdqa	(&QWP(0x10,"ebx"),"xmm1");
+	&pmuludq("xmm1","xmm7");		# a[1]*b[0]
+
+	 &movq	("xmm4","xmm0");		# clear upper 64 bits
+	 &pslldq("xmm4",6);
+	 &paddq	("xmm4","xmm0");
+	 &movdqa("xmm5","xmm4");
+	 &psrldq("xmm4",10);			# upper 32 bits of a[0]*b[0]
+	 &pand	("xmm5","xmm6");		# lower 32 bits of a[0]*b[0]
+
+	# Upper half of a[0]*b[i] is carried into next multiplication
+	# iteration, while lower one "participates" in actual reduction.
+	# Normally latter is done by accumulating result of multiplication
+	# of modulus by "magic" digit, but thanks to special form of modulus
+	# and "magic" digit it can be performed only with additions and
+	# subtractions (see note in IALU section below). Note that we are
+	# not bothered with carry bits, they are accumulated in "flatten"
+	# phase after all multiplications and reductions.
+
+	&movd	("xmm3",&DWP(4*3,"esi"));
+	&pshufd	("xmm2","xmm2",0b11001100);
+	&movdqa	(&QWP(0x20,"ebx"),"xmm2");
+	&pmuludq("xmm2","xmm7");		# a[2]*b[0]
+	 &paddq	("xmm1","xmm4");		# a[1]*b[0]+hw(a[0]*b[0]), carry
+	&movdqa	(&QWP(0x00,"esp"),"xmm1");	# t[0]
+
+	&movd	("xmm0",&DWP(4*4,"esi"));
+	&pshufd	("xmm3","xmm3",0b11001100);
+	&movdqa	(&QWP(0x30,"ebx"),"xmm3");
+	&pmuludq("xmm3","xmm7");		# a[3]*b[0]
+	&movdqa	(&QWP(0x10,"esp"),"xmm2");
+
+	&movd	("xmm1",&DWP(4*5,"esi"));
+	&pshufd	("xmm0","xmm0",0b11001100);
+	&movdqa	(&QWP(0x40,"ebx"),"xmm0");
+	&pmuludq("xmm0","xmm7");		# a[4]*b[0]
+	 &paddq	("xmm3","xmm5");		# a[3]*b[0]+lw(a[0]*b[0]), reduction step
+	&movdqa	(&QWP(0x20,"esp"),"xmm3");
+
+	&movd	("xmm2",&DWP(4*6,"esi"));
+	&pshufd	("xmm1","xmm1",0b11001100);
+	&movdqa	(&QWP(0x50,"ebx"),"xmm1");
+	&pmuludq("xmm1","xmm7");		# a[5]*b[0]
+	&movdqa	(&QWP(0x30,"esp"),"xmm0");
+	 &pshufd("xmm4","xmm5",0b10110001);	# xmm4 = xmm5<<32, reduction step
+
+	&movd	("xmm3",&DWP(4*7,"esi"));
+	&pshufd	("xmm2","xmm2",0b11001100);
+	&movdqa	(&QWP(0x60,"ebx"),"xmm2");
+	&pmuludq("xmm2","xmm7");		# a[6]*b[0]
+	&movdqa	(&QWP(0x40,"esp"),"xmm1");
+	 &psubq	("xmm4","xmm5");		# xmm4 = xmm5*0xffffffff, reduction step
+
+	&movd	("xmm0",&DWP(0,"ebp"));		# b[1] -> 0000.00xy
+	&pshufd	("xmm3","xmm3",0b11001100);
+	&movdqa	(&QWP(0x70,"ebx"),"xmm3");
+	&pmuludq("xmm3","xmm7");		# a[7]*b[0]
+
+	&pshuflw("xmm7","xmm0",0b11011100);	# 0000.00xy -> 0000.0x0y
+	&movdqa	("xmm0",&QWP(0x00,"ebx"));	# pre-load converted a[0]
+	&pshufd	("xmm7","xmm7",0b11011100);	# 0000.0x0y -> 000x.000y
+
+	&mov	("ecx",6);
+	&lea	("ebp",&DWP(4,"ebp"));
+	&jmp	(&label("madd_sse2"));
+
+&set_label("madd_sse2",16);
+	 &paddq	("xmm2","xmm5");		# a[6]*b[i-1]+lw(a[0]*b[i-1]), reduction step [modulo-scheduled]
+	 &paddq	("xmm3","xmm4");		# a[7]*b[i-1]+lw(a[0]*b[i-1])*0xffffffff, reduction step [modulo-scheduled]
+	&movdqa	("xmm1",&QWP(0x10,"ebx"));
+	&pmuludq("xmm0","xmm7");		# a[0]*b[i]
+	 &movdqa(&QWP(0x50,"esp"),"xmm2");
+
+	&movdqa	("xmm2",&QWP(0x20,"ebx"));
+	&pmuludq("xmm1","xmm7");		# a[1]*b[i]
+	 &movdqa(&QWP(0x60,"esp"),"xmm3");
+	&paddq	("xmm0",&QWP(0x00,"esp"));
+
+	&movdqa	("xmm3",&QWP(0x30,"ebx"));
+	&pmuludq("xmm2","xmm7");		# a[2]*b[i]
+	 &movq	("xmm4","xmm0");		# clear upper 64 bits
+	 &pslldq("xmm4",6);
+	&paddq	("xmm1",&QWP(0x10,"esp"));
+	 &paddq	("xmm4","xmm0");
+	 &movdqa("xmm5","xmm4");
+	 &psrldq("xmm4",10);			# upper 33 bits of a[0]*b[i]+t[0]
+
+	&movdqa	("xmm0",&QWP(0x40,"ebx"));
+	&pmuludq("xmm3","xmm7");		# a[3]*b[i]
+	 &paddq	("xmm1","xmm4");		# a[1]*b[i]+hw(a[0]*b[i]), carry
+	&paddq	("xmm2",&QWP(0x20,"esp"));
+	&movdqa	(&QWP(0x00,"esp"),"xmm1");
+
+	&movdqa	("xmm1",&QWP(0x50,"ebx"));
+	&pmuludq("xmm0","xmm7");		# a[4]*b[i]
+	&paddq	("xmm3",&QWP(0x30,"esp"));
+	&movdqa	(&QWP(0x10,"esp"),"xmm2");
+	 &pand	("xmm5","xmm6");		# lower 32 bits of a[0]*b[i]
+
+	&movdqa	("xmm2",&DWP(0x60,"ebx"));
+	&pmuludq("xmm1","xmm7");		# a[5]*b[i]
+	 &paddq	("xmm3","xmm5");		# a[3]*b[i]+lw(a[0]*b[i]), reduction step
+	&paddq	("xmm0",&QWP(0x40,"esp"));
+	&movdqa	(&QWP(0x20,"esp"),"xmm3");
+	 &pshufd("xmm4","xmm5",0b10110001);	# xmm4 = xmm5<<32, reduction step
+
+	&movdqa	("xmm3","xmm7");
+	&pmuludq("xmm2","xmm7");		# a[6]*b[i]
+	 &movd	("xmm7",&DWP(0,"ebp"));		# b[i++] -> 0000.00xy
+	 &lea	("ebp",&DWP(4,"ebp"));
+	&paddq	("xmm1",&QWP(0x50,"esp"));
+	 &psubq	("xmm4","xmm5");		# xmm4 = xmm5*0xffffffff, reduction step
+	&movdqa	(&QWP(0x30,"esp"),"xmm0");
+	 &pshuflw("xmm7","xmm7",0b11011100);	# 0000.00xy -> 0000.0x0y
+
+	&pmuludq("xmm3",&QWP(0x70,"ebx"));	# a[7]*b[i]
+	 &pshufd("xmm7","xmm7",0b11011100);	# 0000.0x0y -> 000x.000y
+	 &movdqa("xmm0",&QWP(0x00,"ebx"));	# pre-load converted a[0]
+	&movdqa	(&QWP(0x40,"esp"),"xmm1");
+	&paddq	("xmm2",&QWP(0x60,"esp"));
+
+	&dec	("ecx");
+	&jnz	(&label("madd_sse2"));
+
+	 &paddq	("xmm2","xmm5");		# a[6]*b[6]+lw(a[0]*b[6]), reduction step [modulo-scheduled]
+	 &paddq	("xmm3","xmm4");		# a[7]*b[6]+lw(a[0]*b[6])*0xffffffff, reduction step [modulo-scheduled]
+	&movdqa	("xmm1",&QWP(0x10,"ebx"));
+	&pmuludq("xmm0","xmm7");		# a[0]*b[7]
+	 &movdqa(&QWP(0x50,"esp"),"xmm2");
+
+	&movdqa	("xmm2",&QWP(0x20,"ebx"));
+	&pmuludq("xmm1","xmm7");		# a[1]*b[7]
+	 &movdqa(&QWP(0x60,"esp"),"xmm3");
+	&paddq	("xmm0",&QWP(0x00,"esp"));
+
+	&movdqa	("xmm3",&QWP(0x30,"ebx"));
+	&pmuludq("xmm2","xmm7");		# a[2]*b[7]
+	 &movq	("xmm4","xmm0");		# clear upper 64 bits
+	 &pslldq("xmm4",6);
+	&paddq	("xmm1",&QWP(0x10,"esp"));
+	 &paddq	("xmm4","xmm0");
+	 &movdqa("xmm5","xmm4");
+	 &psrldq("xmm4",10);			# upper 33 bits of a[0]*b[i]+t[0]
+
+	&movdqa	("xmm0",&QWP(0x40,"ebx"));
+	&pmuludq("xmm3","xmm7");		# a[3]*b[7]
+	 &paddq	("xmm1","xmm4");		# a[1]*b[7]+hw(a[0]*b[7]), carry
+	&paddq	("xmm2",&QWP(0x20,"esp"));
+	&movdqa	(&QWP(0x00,"esp"),"xmm1");
+
+	&movdqa	("xmm1",&QWP(0x50,"ebx"));
+	&pmuludq("xmm0","xmm7");		# a[4]*b[7]
+	&paddq	("xmm3",&QWP(0x30,"esp"));
+	&movdqa	(&QWP(0x10,"esp"),"xmm2");
+	 &pand	("xmm5","xmm6");		# lower 32 bits of a[0]*b[i]
+
+	&movdqa	("xmm2",&DWP(0x60,"ebx"));
+	&pmuludq("xmm1","xmm7");		# a[5]*b[7]
+	 &paddq	("xmm3","xmm5");		# reduction step
+	&paddq	("xmm0",&QWP(0x40,"esp"));
+	&movdqa	(&QWP(0x20,"esp"),"xmm3");
+	 &pshufd("xmm4","xmm5",0b10110001);	# xmm4 = xmm5<<32, reduction step
+
+	&movdqa	("xmm3",&QWP(0x70,"ebx"));
+	&pmuludq("xmm2","xmm7");		# a[6]*b[7]
+	&paddq	("xmm1",&QWP(0x50,"esp"));
+	 &psubq	("xmm4","xmm5");		# xmm4 = xmm5*0xffffffff, reduction step
+	&movdqa	(&QWP(0x30,"esp"),"xmm0");
+
+	&pmuludq("xmm3","xmm7");		# a[7]*b[7]
+	&pcmpeqd("xmm7","xmm7");
+	&movdqa	("xmm0",&QWP(0x00,"esp"));
+	&pslldq	("xmm7",8);
+	&movdqa	(&QWP(0x40,"esp"),"xmm1");
+	&paddq	("xmm2",&QWP(0x60,"esp"));
+
+	 &paddq	("xmm2","xmm5");		# a[6]*b[7]+lw(a[0]*b[7]), reduction step
+	 &paddq	("xmm3","xmm4");		# a[6]*b[7]+lw(a[0]*b[7])*0xffffffff, reduction step
+	 &movdqa(&QWP(0x50,"esp"),"xmm2");
+	 &movdqa(&QWP(0x60,"esp"),"xmm3");
+
+	&movdqa	("xmm1",&QWP(0x10,"esp"));
+	&movdqa	("xmm2",&QWP(0x20,"esp"));
+	&movdqa	("xmm3",&QWP(0x30,"esp"));
+
+	&movq	("xmm4","xmm0");		# "flatten"
+	&pand	("xmm0","xmm7");
+	&xor	("ebp","ebp");
+	&pslldq	("xmm4",6);
+	 &movq	("xmm5","xmm1");
+	&paddq	("xmm0","xmm4");
+	 &pand	("xmm1","xmm7");
+	&psrldq	("xmm0",6);
+	&movd	("eax","xmm0");
+	&psrldq	("xmm0",4);
+
+	&paddq	("xmm5","xmm0");
+	&movdqa	("xmm0",&QWP(0x40,"esp"));
+	&sub	("eax",-1);			# start subtracting modulus,
+						# this is used to determine
+						# if result is larger/smaller
+						# than modulus (see below)
+	&pslldq	("xmm5",6);
+	 &movq	("xmm4","xmm2");
+	&paddq	("xmm1","xmm5");
+	 &pand	("xmm2","xmm7");
+	&psrldq	("xmm1",6);
+	&mov	(&DWP(4*0,"edi"),"eax");
+	&movd	("eax","xmm1");
+	&psrldq	("xmm1",4);
+
+	&paddq	("xmm4","xmm1");
+	&movdqa	("xmm1",&QWP(0x50,"esp"));
+	&sbb	("eax",-1);
+	&pslldq	("xmm4",6);
+	 &movq	("xmm5","xmm3");
+	&paddq	("xmm2","xmm4");
+	 &pand	("xmm3","xmm7");
+	&psrldq	("xmm2",6);
+	&mov	(&DWP(4*1,"edi"),"eax");
+	&movd	("eax","xmm2");
+	&psrldq	("xmm2",4);
+
+	&paddq	("xmm5","xmm2");
+	&movdqa	("xmm2",&QWP(0x60,"esp"));
+	&sbb	("eax",-1);
+	&pslldq	("xmm5",6);
+	 &movq	("xmm4","xmm0");
+	&paddq	("xmm3","xmm5");
+	 &pand	("xmm0","xmm7");
+	&psrldq	("xmm3",6);
+	&mov	(&DWP(4*2,"edi"),"eax");
+	&movd	("eax","xmm3");
+	&psrldq	("xmm3",4);
+
+	&paddq	("xmm4","xmm3");
+	&sbb	("eax",0);
+	&pslldq	("xmm4",6);
+	 &movq	("xmm5","xmm1");
+	&paddq	("xmm0","xmm4");
+	 &pand	("xmm1","xmm7");
+	&psrldq	("xmm0",6);
+	&mov	(&DWP(4*3,"edi"),"eax");
+	&movd	("eax","xmm0");
+	&psrldq	("xmm0",4);
+
+	&paddq	("xmm5","xmm0");
+	&sbb	("eax",0);
+	&pslldq	("xmm5",6);
+	 &movq	("xmm4","xmm2");
+	&paddq	("xmm1","xmm5");
+	 &pand	("xmm2","xmm7");
+	&psrldq	("xmm1",6);
+	&movd	("ebx","xmm1");
+	&psrldq	("xmm1",4);
+	&mov	("esp","edx");
+
+	&paddq	("xmm4","xmm1");
+	&pslldq	("xmm4",6);
+	&paddq	("xmm2","xmm4");
+	&psrldq	("xmm2",6);
+	&movd	("ecx","xmm2");
+	&psrldq	("xmm2",4);
+	&sbb	("ebx",0);
+	&movd	("edx","xmm2");
+	&pextrw	("esi","xmm2",2);		# top-most overflow bit
+	&sbb	("ecx",1);
+	&sbb	("edx",-1);
+	&sbb	("esi",0);			# borrow from subtraction
+
+	# Final step is "if result > mod, subtract mod", and at this point
+	# we have result - mod written to output buffer, as well as borrow
+	# bit from this subtraction, and if borrow bit is set, we add
+	# modulus back.
+	#
+	# Note that because mod has special form, i.e. consists of
+	# 0xffffffff, 1 and 0s, we can conditionally synthesize it by
+	# assigning borrow bit to one register, %ebp, and its negative
+	# to another, %esi. But we started by calculating %esi...
+
+	&sub	("ebp","esi");
+	&add	(&DWP(4*0,"edi"),"esi");	# add modulus or zero
+	&adc	(&DWP(4*1,"edi"),"esi");
+	&adc	(&DWP(4*2,"edi"),"esi");
+	&adc	(&DWP(4*3,"edi"),0);
+	&adc	("eax",0);
+	&adc	("ebx",0);
+	&mov	(&DWP(4*4,"edi"),"eax");
+	&adc	("ecx","ebp");
+	&mov	(&DWP(4*5,"edi"),"ebx");
+	&adc	("edx","esi");
+	&mov	(&DWP(4*6,"edi"),"ecx");
+	&mov	(&DWP(4*7,"edi"),"edx");
+
+	&ret	();
+
+&set_label("mul_mont_ialu",16);			}
+
+	########################################
+	# IALU code path suitable for all CPUs.
+	########################################
+	# stack layout:
+	# +------------------------------------+< %esp
+	# | 8 32-bit temporary words, accessed |
+	# | as circular buffer                 |
+	# .                                    .
+	# .                                    .
+	# +------------------------------------+< +32
+	# | offloaded destination pointer      |
+	# +------------------------------------+
+	# | unused                             |
+	# +------------------------------------+< +40
+	&sub	("esp",10*4);
+
+	&mov	("eax",&DWP(0*4,"esi"));		# a[0]
+	&mov	("ebx",&DWP(0*4,"ebp"));		# b[0]
+	&mov	(&DWP(8*4,"esp"),"edi");		# off-load dst ptr
+
+	&mul	("ebx");				# a[0]*b[0]
+	&mov	(&DWP(0*4,"esp"),"eax");		# t[0]
+	&mov	("eax",&DWP(1*4,"esi"));
+	&mov	("ecx","edx")
+
+	&mul	("ebx");				# a[1]*b[0]
+	&add	("ecx","eax");
+	&mov	("eax",&DWP(2*4,"esi"));
+	&adc	("edx",0);
+	&mov	(&DWP(1*4,"esp"),"ecx");		# t[1]
+	&mov	("ecx","edx");
+
+	&mul	("ebx");				# a[2]*b[0]
+	&add	("ecx","eax");
+	&mov	("eax",&DWP(3*4,"esi"));
+	&adc	("edx",0);
+	&mov	(&DWP(2*4,"esp"),"ecx");		# t[2]
+	&mov	("ecx","edx");
+
+	&mul	("ebx");				# a[3]*b[0]
+	&add	("ecx","eax");
+	&mov	("eax",&DWP(4*4,"esi"));
+	&adc	("edx",0);
+	&mov	(&DWP(3*4,"esp"),"ecx");		# t[3]
+	&mov	("ecx","edx");
+
+	&mul	("ebx");				# a[4]*b[0]
+	&add	("ecx","eax");
+	&mov	("eax",&DWP(5*4,"esi"));
+	&adc	("edx",0);
+	&mov	(&DWP(4*4,"esp"),"ecx");		# t[4]
+	&mov	("ecx","edx");
+
+	&mul	("ebx");				# a[5]*b[0]
+	&add	("ecx","eax");
+	&mov	("eax",&DWP(6*4,"esi"));
+	&adc	("edx",0);
+	&mov	(&DWP(5*4,"esp"),"ecx");		# t[5]
+	&mov	("ecx","edx");
+
+	&mul	("ebx");				# a[6]*b[0]
+	&add	("ecx","eax");
+	&mov	("eax",&DWP(7*4,"esi"));
+	&adc	("edx",0);
+	&mov	(&DWP(6*4,"esp"),"ecx");		# t[6]
+	&mov	("ecx","edx");
+
+	&xor	("edi","edi");				# initial top-most carry
+	&mul	("ebx");				# a[7]*b[0]
+	&add	("ecx","eax");				# t[7]
+	&mov	("eax",&DWP(0*4,"esp"));		# t[0]
+	&adc	("edx",0);				# t[8]
+
+for ($i=0;$i<7;$i++) {
+	my $j=$i+1;
+
+	# Reduction iteration is normally performed by accumulating
+	# result of multiplication of modulus by "magic" digit [and
+	# omitting least significant word, which is guaranteed to
+	# be 0], but thanks to special form of modulus and "magic"
+	# digit being equal to least significant word, it can be
+	# performed with additions and subtractions alone. Indeed:
+	#
+	#        ffff.0001.0000.0000.0000.ffff.ffff.ffff
+	# *                                         abcd
+	# + xxxx.xxxx.xxxx.xxxx.xxxx.xxxx.xxxx.xxxx.abcd
+	#
+	# Now observing that ff..ff*x = (2^n-1)*x = 2^n*x-x, we
+	# rewrite above as:
+	#
+	#   xxxx.xxxx.xxxx.xxxx.xxxx.xxxx.xxxx.xxxx.abcd
+	# + abcd.0000.abcd.0000.0000.abcd.0000.0000.0000
+	# -      abcd.0000.0000.0000.0000.0000.0000.abcd
+	#
+	# or marking redundant operations:
+	#
+	#   xxxx.xxxx.xxxx.xxxx.xxxx.xxxx.xxxx.xxxx.----
+	# + abcd.0000.abcd.0000.0000.abcd.----.----.----
+	# -      abcd.----.----.----.----.----.----.----
+
+	&add	(&DWP((($i+3)%8)*4,"esp"),"eax");	# t[3]+=t[0]
+	&adc	(&DWP((($i+4)%8)*4,"esp"),0);		# t[4]+=0
+	&adc	(&DWP((($i+5)%8)*4,"esp"),0);		# t[5]+=0
+	&adc	(&DWP((($i+6)%8)*4,"esp"),"eax");	# t[6]+=t[0]
+	&adc	("ecx",0);				# t[7]+=0
+	&adc	("edx","eax");				# t[8]+=t[0]
+	&adc	("edi",0);				# top-most carry
+	 &mov	("ebx",&DWP($j*4,"ebp"));		# b[i]
+	&sub	("ecx","eax");				# t[7]-=t[0]
+	 &mov	("eax",&DWP(0*4,"esi"));		# a[0]
+	&sbb	("edx",0);				# t[8]-=0
+	&mov	(&DWP((($i+7)%8)*4,"esp"),"ecx");
+	&sbb	("edi",0);				# top-most carry,
+							# keep in mind that
+							# netto result is
+							# *addition* of value
+							# with (abcd<<32)-abcd
+							# on top, so that
+							# underflow is
+							# impossible, because
+							# (abcd<<32)-abcd
+							# doesn't underflow
+	&mov	(&DWP((($i+8)%8)*4,"esp"),"edx");
+
+	&mul	("ebx");				# a[0]*b[i]
+	&add	("eax",&DWP((($j+0)%8)*4,"esp"));
+	&adc	("edx",0);
+	&mov	(&DWP((($j+0)%8)*4,"esp"),"eax");
+	&mov	("eax",&DWP(1*4,"esi"));
+	&mov	("ecx","edx")
+
+	&mul	("ebx");				# a[1]*b[i]
+	&add	("ecx",&DWP((($j+1)%8)*4,"esp"));
+	&adc	("edx",0);
+	&add	("ecx","eax");
+	&adc	("edx",0);
+	&mov	("eax",&DWP(2*4,"esi"));
+	&mov	(&DWP((($j+1)%8)*4,"esp"),"ecx");
+	&mov	("ecx","edx");
+
+	&mul	("ebx");				# a[2]*b[i]
+	&add	("ecx",&DWP((($j+2)%8)*4,"esp"));
+	&adc	("edx",0);
+	&add	("ecx","eax");
+	&adc	("edx",0);
+	&mov	("eax",&DWP(3*4,"esi"));
+	&mov	(&DWP((($j+2)%8)*4,"esp"),"ecx");
+	&mov	("ecx","edx");
+
+	&mul	("ebx");				# a[3]*b[i]</