mpn/x86_64/mode1o.asm

   1 dnl  AMD64 mpn_modexact_1_odd -- exact division style remainder.
   2
   3 dnl  Copyright 2000, 2001, 2002, 2003, 2004, 2005, 2006 Free Software
   4 dnl  Foundation, Inc.
   5 dnl
   6 dnl  This file is part of the GNU MP Library.
   7 dnl
   8 dnl  The GNU MP Library is free software; you can redistribute it and/or
   9 dnl  modify it under the terms of the GNU Lesser General Public License as
  10 dnl  published by the Free Software Foundation; either version 3 of the
  11 dnl  License, or (at your option) any later version.
  12 dnl
  13 dnl  The GNU MP Library is distributed in the hope that it will be useful,
  14 dnl  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 dnl  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16 dnl  Lesser General Public License for more details.
  17 dnl
  18 dnl  You should have received a copy of the GNU Lesser General Public License
  19 dnl  along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.
  20
  21 include(`../config.m4')
  22
  23
  24 C            cycles/limb
  25 C K8,K9:        10
  26 C K10:          10
  27 C P4:           33
  28 C P6 core2:     13
  29 C P6 corei7:    14.5
  30 C P6 Atom:      35
  31
  32
  33 C mp_limb_t mpn_modexact_1_odd (mp_srcptr src, mp_size_t size,
  34 C                               mp_limb_t divisor);
  35 C mp_limb_t mpn_modexact_1c_odd (mp_srcptr src, mp_size_t size,
  36 C                                mp_limb_t divisor, mp_limb_t carry);
  37 C
  38 C
  39 C The dependent chain in the main loop is
  40 C
  41 C                            cycles
  42 C       subq    %rdx, %rax      1
  43 C       imulq   %r9, %rax       4
  44 C       mulq    %r8             5
  45 C                             ----
  46 C       total                  10
  47 C
  48 C The movq load from src seems to need to be scheduled back before the jz to
  49 C achieve this speed, out-of-order execution apparently can't completely
  50 C hide the latency otherwise.
  51 C
  52 C The l=src[i]-cbit step is rotated back too, since that allows us to avoid
  53 C it for the first iteration (where there's no cbit).
  54 C
  55 C The code alignment used (32-byte) for the loop also seems necessary.
  56 C Without that the non-PIC case has adcq crossing the 0x60 offset,
  57 C apparently making it run at 11 cycles instead of 10.
  58 C
  59 C Not done:
  60 C
  61 C divq for size==1 was measured at about 79 cycles, compared to the inverse
  62 C at about 25 cycles (both including function call overheads), so that's not
  63 C used.
  64 C
  65 C Enhancements:
  66 C
  67 C For PIC, we shouldn't really need the GOT fetch for binvert_limb_table,
  68 C it'll be in rodata or text in libgmp.so and can be accessed directly %rip
  69 C relative.  This would be for small model only (something we don't
  70 C presently detect, but which is all that gcc 3.3.3 supports), since 8-byte
  71 C PC-relative relocations are apparently not available.  Some rough
  72 C experiments with binutils 2.13 looked worrylingly like it might come out
  73 C with an unwanted text segment relocation though, even with ".protected".
  74
  75
  76 ASM_START()
  77         TEXT
  78         ALIGN(32)
  79 PROLOGUE(mpn_modexact_1_odd)
  80
  81         movl    $0, %ecx
  82
  83 PROLOGUE(mpn_modexact_1c_odd)
  84
  85         C rdi   src
  86         C rsi   size
  87         C rdx   divisor
  88         C rcx   carry
  89
  90         movq    %rdx, %r8               C d
  91         shrl    %edx                    C d/2
  92 ifdef(`PIC',`
  93         movq    binvert_limb_table@GOTPCREL(%rip), %r9
  94 ',`
  95         movabsq $binvert_limb_table, %r9
  96 ')
  97
  98         andl    $127, %edx
  99         movq    %rcx, %r10              C initial carry
 100
 101         movzbl  (%r9,%rdx), %edx        C inv 8 bits
 102
 103         movq    (%rdi), %rax            C src[0]
 104         leaq    (%rdi,%rsi,8), %r11     C src end
 105         movq    %r8, %rdi               C d, made available to imull
 106
 107         leal    (%rdx,%rdx), %ecx       C 2*inv
 108         imull   %edx, %edx              C inv*inv
 109
 110         negq    %rsi                    C -size
 111
 112         imull   %edi, %edx              C inv*inv*d
 113
 114         subl    %edx, %ecx              C inv = 2*inv - inv*inv*d, 16 bits
 115
 116         leal    (%rcx,%rcx), %edx       C 2*inv
 117         imull   %ecx, %ecx              C inv*inv
 118
 119         imull   %edi, %ecx              C inv*inv*d
 120
 121         subl    %ecx, %edx              C inv = 2*inv - inv*inv*d, 32 bits
 122         xorl    %ecx, %ecx              C initial cbit
 123
 124         leaq    (%rdx,%rdx), %r9        C 2*inv
 125         imulq   %rdx, %rdx              C inv*inv
 126
 127         imulq   %r8, %rdx               C inv*inv*d
 128
 129         subq    %rdx, %r9               C inv = 2*inv - inv*inv*d, 64 bits
 130         movq    %r10, %rdx              C initial climb
 131
 132         ASSERT(e,`      C d*inv == 1 mod 2^64
 133         movq    %r8, %r10
 134         imulq   %r9, %r10
 135         cmpq    $1, %r10')
 136
 137         incq    %rsi
 138         jz      L(one)
 139
 140
 141         ALIGN(16)
 142 L(top):
 143         C rax   l = src[i]-cbit
 144         C rcx   new cbit, 0 or 1
 145         C rdx   climb, high of last product
 146         C rsi   counter, limbs, negative
 147         C rdi
 148         C r8    divisor
 149         C r9    inverse
 150         C r11   src end ptr
 151
 152         subq    %rdx, %rax              C l = src[i]-cbit - climb
 153
 154         adcq    $0, %rcx                C more cbit
 155         imulq   %r9, %rax               C q = l * inverse
 156
 157         mulq    %r8                     C climb = high (q * d)
 158
 159         movq    (%r11,%rsi,8), %rax     C src[i+1]
 160         subq    %rcx, %rax              C next l = src[i+1] - cbit
 161         setc    %cl                     C new cbit
 162
 163         incq    %rsi
 164         jnz     L(top)
 165
 166
 167 L(one):
 168         subq    %rdx, %rax              C l = src[i]-cbit - climb
 169
 170         adcq    $0, %rcx                C more cbit
 171         imulq   %r9, %rax               C q = l * inverse
 172
 173         mulq    %r8                     C climb = high (q * d)
 174
 175         leaq    (%rcx,%rdx), %rax       C climb+cbit
 176         ret
 177
 178 EPILOGUE(mpn_modexact_1c_odd)
 179 EPILOGUE(mpn_modexact_1_odd)