Ticket #1: submul_1.yasm

;; *********************************************************************
;;  Intel64 mpn_submul_1 -- Multiply a limb vector with a limb and
;;  subtract the result from a second limb vector.
;;
;;  Copyright (C) 2006  Jason Worth Martin <jason.worth.martin@gmail.com>
;;
;;  This program is free software; you can redistribute it and/or modify
;;  it under the terms of the GNU General Public License as published by
;;  the Free Software Foundation; either version 2 of the License, or
;;  (at your option) any later version.
;;
;;  This program is distributed in the hope that it will be useful,
;;  but WITHOUT ANY WARRANTY; without even the implied warranty of
;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;  GNU General Public License for more details.
;;
;;  You should have received a copy of the GNU General Public License along
;;  with this program; if not, write to the Free Software Foundation, Inc.,
;;  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
;;
;; **************************************************************************
;;
;;
;; CREDITS
;;
;; This code is based largely on Pierrick Gaudry's excellent assembly
;; support for the AMD64 architecture.  (Note that Intel64 and AMD64,
;; while using the same instruction set, have very different
;; microarchitectures.  So, this code performs very poorly on AMD64
;; machines even though it is near-optimal on Intel64.)
;;
;; Roger Golliver works for Intel and provided insightful improvements
;; particularly in using the "lea" instruction to perform additions
;; and register-to-register moves.
;;
;; Eric Bainville has a brilliant exposition of optimizing arithmetic for
;; AMD64 (http://www.bealto.it).  I adapted many of the ideas he
;; describes to Intel64.
;;
;; Agner Fog is a demigod in the x86 world.  If you are reading assembly
;; code files and you haven't heard of Agner Fog, then take a minute to
;; look over his software optimization manuals (http://www.agner.org/).
;; They are superb.
;;
;; *********************************************************************


;; With a 4-way unroll the code has
;;
;;              cycles/limb
;; Hammer:           4.6
;; Woodcrest:        4.6
;;
;; With increased unrolling, it appears to converge to 4 cycles/limb
;; on Intel Core 2 machines.  I believe that this is optimal, however
;; it requires such absurd unrolling that it becomes unusable for all
;; but the largest inputs.  A 4-way unroll seems like a good balance
;; to me because then commonly used input sizes (e.g. 1024bit Public
;; keys) still benifit from the speed up.

;;
;; This is just a check to see if we are in my code testing sandbox
;; or if we are actually in the GMP source tree
;;
%ifdef  __JWM_Test_Code__
%include 'yasm_mac.inc'
%else
%include '../yasm_mac.inc'
%endif



;; *********************************************************************
;; *********************************************************************
;;
;; Here are the important macro parameters for the code
;;
;;      BpL is Bytes per Limb (8 since this is 64bit code)
;;
;;      UNROLL_SIZE is a power of 2 for which we unroll the code.
;;                  possible values are 2,4,8,15,..., 256.  A reasonable
;;                  value is probably 4.  If really large inputs
;;                  are expected, then 16 is probably good.  Larger
;;                  values are really only useful for flashy
;;                  benchmarks and testing asymptotic behavior.
;;
;;      THRESHOLD is the minimum number of limbs needed before we bother
;;                using the complicated loop.  A reasonable value is
;;                2*UNROLL_SIZE + 6
;;
;; *********************************************************************
;; *********************************************************************
%define BpL     8
%define UNROLL_SIZE     4
%define UNROLL_MASK     UNROLL_SIZE-1
%define THRESHOLD       2*UNROLL_SIZE+6

;; Here is a convenient Macro for addressing
;; memory.  Entries of the form
;;
;;      ADDR(ptr,index,displacement)
;;
;; get converted to
;;
;;      [displacement*BpL + ptr + index*BpL]
;;
%define ADDR(a,b,c)     [c*BpL+a+b*BpL]


;; Register     Usage
;; --------     -----
;; rax          low word from mul
;; rbx*
;; rcx          s2limb
;; rdx          high word from mul
;; rsi          s1p
;; rdi          rp
;; rbp*         Base Pointer
;; rsp*         Stack Pointer
;; r8           A_x
;; r9           A_y
;; r10          A_z
;; r11          B_x
;; r12*         B_y
;; r13*         B_z
;; r14*         temp
;; r15*         index
;;
;; * indicates that the register must be
;; preserved for the caller.
%define s2limb  rcx
%define s1p     rsi
%define rp      rdi
%define A_x     r8
%define A_y     r9
%define A_z     r10
%define B_x     r11
%define B_y     r12
%define B_z     r13
%define temp    r14
%define index   r15

        
;; INPUT PARAMETERS
;; rp           rdi
;; s1p          rsi
;; n            rdx
;; s2limb       rcx
        BITS    64
G_EXPORT mpn_submul_1
                                        ;; Compare the limb count
                                        ;; with the threshold value.
                                        ;; If the limb count is small
                                        ;; we just use the small loop,
                                        ;; otherwise we jump to the
                                        ;; more complicated loop.
        cmp     rdx,THRESHOLD
        jge     L_mpn_submul_1_main_loop_prep
        mov     r11,rdx
        lea     rsi,[rsi+rdx*8]
        lea     rdi,[rdi+rdx*8]
        neg     r11
        xor     r8, r8
        xor     r10, r10
        jmp     L_mpn_submul_1_small_loop
        
        align   16
L_mpn_submul_1_small_loop:
        mov     rax,[rsi+r11*8]
        mul     rcx
        add     rax,r8
        adc     rdx,r10
        sub     [rdi+r11*8],rax
        mov     r8,r10
        adc     r8,rdx
        inc     r11
        jne     L_mpn_submul_1_small_loop

        mov     rax,r8
        ret

L_mpn_submul_1_main_loop_prep:
        push    r15
        push    r14
        push    r13
        push    r12
                                ;; If n is even, we need to do three
                                ;; pre-multiplies, if n is odd we only
                                ;; need to do two.
        mov     temp,rdx
        mov     index,0
        mov     A_x,0
        mov     A_y,0
        and     rdx,1
        jnz     L_mpn_submul_1_odd_n

                                        ;; Case n is even
        mov     rax,ADDR(s1p,index,0)
        mul     s2limb
        sub     ADDR(rp,index,0),rax
        adc     A_x,rdx
        add     index,1
                                        ;; At this point
                                        ;;  temp = n (even)
                                        ;; index = 1

L_mpn_submul_1_odd_n:
                                        ;; Now
                                        ;; temp = n
                                        ;; index = 1 if n even
                                        ;;       = 0 if n odd
                                        ;;
        mov     rax,ADDR(s1p,index,0)
        mul     s2limb
        mov     A_z,ADDR(rp,index,0)
        add     A_x,rax
        adc     A_y,rdx

        mov     rax,ADDR(s1p,index,1)
        mul     s2limb
        mov     B_z,ADDR(rp,index,1)
        mov     B_x,rax
        mov     B_y,rdx
        mov     rax,ADDR(s1p,index,2)

        add     index,3
        lea     s1p,ADDR(s1p,temp,-1)
        lea     rp,ADDR(rp,temp,-1)
        neg     temp
        add     index,temp
                                ;; At this point:
                                ;; s1p   = address of last s1limb
                                ;; rp    = address of last rplimb
                                ;; temp  = -n
                                ;; index = 4 - n if n even
                                ;;       = 3 - n if n odd
                                ;;
                                ;; So, index is a (negative) even
                                ;; number.
                                ;;
                                ;; *****************************************
                                ;; ATTENTION:
                                ;;
                                ;; From here on, I will use array
                                ;; indexing notation in the comments
                                ;; because it is convenient.  So, I
                                ;; will pretend that index is positive
                                ;; because then a comment like
                                ;;      B_z = rp[index-1]
                                ;; is easier to read.
                                ;; However, keep in mind that index is
                                ;; actually a negative number indexing
                                ;; back from the end of the array.
                                ;; This is a common trick to remove one
                                ;; compare operation from the main loop.
                                ;; *****************************************

                                ;;
                                ;; Now we enter a spin-up loop the
                                ;; will make sure that the index is
                                ;; a multiple of UNROLL_SIZE before
                                ;; going to our main unrolled loop.
        mov     temp,index
        neg     temp
        and     temp,UNROLL_MASK
        jz      L_mpn_submul_1_main_loop
        shr     temp,1
L_mpn_submul_1_main_loop_spin_up:
                                ;; At this point we should have:
                                ;;
                                ;; A_x = low_mul[index-2] + carry_in
                                ;; A_y = high_mul[index-2] + CF
                                ;; A_z = rp[index-2]
                                ;;
                                ;; B_x = low_mul[index-1]
                                ;; B_y = high_mul[index-1]
                                ;; B_z = rp[index-1]
                                ;;
                                ;; rax = s1p[index]
        mul     s2limb
        sub     A_z,A_x
        lea     A_x,[rax]
        mov     rax,ADDR(s1p,index,1)
        adc     B_x,A_y
        mov     ADDR(rp,index,-2),A_z
        mov     A_z,ADDR(rp,index,0)
        adc     B_y,0
        lea     A_y,[rdx]
                                ;; At this point we should have:
                                ;;
                                ;; B_x = low_mul[index-1] + carry_in
                                ;; B_y = high_mul[index-1] + CF
                                ;; B_z = rp[index-1]
                                ;;
                                ;; A_x = low_mul[index]
                                ;; A_y = high_mul[index]
                                ;; A_z = rp[index]
                                ;;
                                ;; rax = s1p[index+1]
        mul     s2limb
        sub     B_z,B_x
        lea     B_x,[rax]
        mov     rax,ADDR(s1p,index,2)
        adc     A_x,B_y
        mov     ADDR(rp,index,-1),B_z
        mov     B_z,ADDR(rp,index,1)
        adc     A_y,0
        lea     B_y,[rdx]

        add     index,2
        sub     temp,1
        jnz     L_mpn_submul_1_main_loop_spin_up
        
        jmp     L_mpn_submul_1_main_loop
        
        align   16
L_mpn_submul_1_main_loop:
                                ;; The code here is really the same
                                ;; logic as the spin-up loop.  It's
                                ;; just been unrolled.
%assign unroll_index 0
%rep    UNROLL_SIZE/2
        mul     s2limb
        sub     A_z,A_x
        lea     A_x,[rax]
        mov     rax,ADDR(s1p,index,(2*unroll_index+1))
        adc     B_x,A_y
        mov     ADDR(rp,index,(2*unroll_index-2)),A_z
        mov     A_z,ADDR(rp,index,(2*unroll_index))
        adc     B_y,0
        lea     A_y,[rdx]

        mul     s2limb
        sub     B_z,B_x
        lea     B_x,[rax]
        mov     rax,ADDR(s1p,index,(2*unroll_index+2))
        adc     A_x,B_y
        mov     ADDR(rp,index,(2*unroll_index-1)),B_z
        mov     B_z,ADDR(rp,index,(2*unroll_index+1))
        adc     A_y,0
        lea     B_y,[rdx]
%assign unroll_index    unroll_index+1
%endrep


        add     index,UNROLL_SIZE
        jnz     L_mpn_submul_1_main_loop

L_mpn_submul_1_finish:  
                                ;; At this point we should have:
                                ;;
                                ;; index = n-1
                                ;;
                                ;; A_x = low_mul[index-2] + carry_in
                                ;; A_y = high_mul[index-2] + CF
                                ;; A_z = rp[index-2]
                                ;;
                                ;; B_x = low_mul[index-1]
                                ;; B_y = high_mul[index-1]
                                ;; B_z = rp[index-1]
                                ;;
                                ;; rax = s1p[index]
        mul     s2limb
        sub     A_z,A_x
        mov     A_x,rax
        mov     ADDR(rp,index,-2),A_z
        mov     A_z,ADDR(rp,index,0)
        adc     B_x,A_y
        adc     B_y,0
        mov     A_y,rdx
                                ;; At this point we should have:
                                ;;
                                ;; index = n-1
                                ;;
                                ;; B_x = low_mul[index-1] + carry_in
                                ;; B_y = high_mul[index-1] + CF
                                ;; B_z = rp[index-1]
                                ;;
                                ;; A_x = low_mul[index]
                                ;; A_y = high_mul[index]
                                ;; A_z = rp[index]
        sub     B_z,B_x
        mov     ADDR(rp,index,-1),B_z
        adc     A_x,B_y
        adc     A_y,0
                                ;; At this point we should have:
                                ;;
                                ;; index = n-1
                                ;;
                                ;; A_x = low_mul[index] + carry_in
                                ;; A_y = high_mul[index] + CF
                                ;; A_z = rp[index]
                                ;;
        sub     A_z,A_x
        mov     ADDR(rp,index,0),A_z
        adc     A_y,0

        mov     rax,A_y
        pop     r12
        pop     r13
        pop     r14
        pop     r15
        ret