--- /dev/null
+dnl Intel Pentium-4 mpn_submul_1 -- Multiply a limb vector with a limb and
+dnl subtract the result from a second limb vector.
+
+dnl Copyright 2001, 2002 Free Software Foundation, Inc.
+dnl
+dnl This file is part of the GNU MP Library.
+dnl
+dnl The GNU MP Library is free software; you can redistribute it and/or
+dnl modify it under the terms of the GNU Lesser General Public License as
+dnl published by the Free Software Foundation; either version 3 of the
+dnl License, or (at your option) any later version.
+dnl
+dnl The GNU MP Library is distributed in the hope that it will be useful,
+dnl but WITHOUT ANY WARRANTY; without even the implied warranty of
+dnl MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+dnl Lesser General Public License for more details.
+dnl
+dnl You should have received a copy of the GNU Lesser General Public License
+dnl along with the GNU MP Library. If not, see http://www.gnu.org/licenses/.
+
+include(`../config.m4')
+
+
+C P4: 7 cycles/limb, unstable timing, at least on early Pentium4 silicon
+C (stepping 10).
+
+
+C mp_limb_t mpn_submul_1 (mp_ptr dst, mp_srcptr src, mp_size_t size,
+C mp_limb_t mult);
+C mp_limb_t mpn_submul_1c (mp_ptr dst, mp_srcptr src, mp_size_t size,
+C mp_limb_t mult, mp_limb_t carry);
+C
+C This code is not particularly good at 7 c/l. The dependent chain is only
+C 4 c/l and there's only 4 MMX unit instructions, so it's not clear why that
+C speed isn't achieved.
+C
+C The arrangements made here to get a two instruction dependent chain are
+C slightly subtle. In the loop the carry (or borrow rather) is a negative
+C so that a paddq can be used to give a low limb ready to store, and a high
+C limb ready to become the new carry after a psrlq.
+C
+C If the carry was a simple twos complement negative then the psrlq shift
+C would need to bring in 0 bits or 1 bits according to whether the high was
+C zero or non-zero, since a non-zero value would represent a negative
+C needing sign extension. That wouldn't be particularly easy to arrange and
+C certainly would add an instruction to the dependent chain, so instead an
+C offset is applied so that the high limb will be 0xFFFFFFFF+c. With c in
+C the range -0xFFFFFFFF to 0, the value 0xFFFFFFFF+c is in the range 0 to
+C 0xFFFFFFFF and is therefore always positive and can always have 0 bits
+C shifted in, which is what psrlq does.
+C
+C The extra 0xFFFFFFFF must be subtracted before c is used, but that can be
+C done off the dependent chain. The total adjustment then is to add
+C 0xFFFFFFFF00000000 to offset the new carry, and subtract
+C 0x00000000FFFFFFFF to remove the offset from the current carry, for a net
+C add of 0xFFFFFFFE00000001. In the code this is applied to the destination
+C limb when fetched.
+C
+C It's also possible to view the 0xFFFFFFFF adjustment as a ones-complement
+C negative, which is how it's undone for the return value, but that doesn't
+C seem as clear.
+
+defframe(PARAM_CARRY, 20)
+defframe(PARAM_MULTIPLIER,16)
+defframe(PARAM_SIZE, 12)
+defframe(PARAM_SRC, 8)
+defframe(PARAM_DST, 4)
+
+ TEXT
+ ALIGN(16)
+
+PROLOGUE(mpn_submul_1c)
+deflit(`FRAME',0)
+ movd PARAM_CARRY, %mm1
+ jmp L(start_1c)
+EPILOGUE()
+
+PROLOGUE(mpn_submul_1)
+deflit(`FRAME',0)
+ pxor %mm1, %mm1 C initial borrow
+
+L(start_1c):
+ movl PARAM_SRC, %eax
+ pcmpeqd %mm0, %mm0
+
+ movd PARAM_MULTIPLIER, %mm7
+ pcmpeqd %mm6, %mm6
+
+ movl PARAM_DST, %edx
+ psrlq $32, %mm0 C 0x00000000FFFFFFFF
+
+ movl PARAM_SIZE, %ecx
+ psllq $32, %mm6 C 0xFFFFFFFF00000000
+
+ psubq %mm0, %mm6 C 0xFFFFFFFE00000001
+
+ psubq %mm1, %mm0 C 0xFFFFFFFF - borrow
+
+
+ C eax src, incrementing
+ C ebx
+ C ecx loop counter, decrementing
+ C edx dst, incrementing
+ C
+ C mm0 0xFFFFFFFF - borrow
+ C mm6 0xFFFFFFFE00000001
+ C mm7 multiplier
+
+L(loop):
+ movd (%eax), %mm1 C src
+ leal 4(%eax), %eax
+ movd (%edx), %mm2 C dst
+ paddq %mm6, %mm2 C add 0xFFFFFFFE00000001
+ pmuludq %mm7, %mm1
+ psubq %mm1, %mm2 C prod
+ paddq %mm2, %mm0 C borrow
+ subl $1, %ecx
+ movd %mm0, (%edx) C result
+ psrlq $32, %mm0
+ leal 4(%edx), %edx
+ jnz L(loop)
+
+ movd %mm0, %eax
+ notl %eax
+ emms
+ ret
+
+EPILOGUE()