mpn/x86/pentium4/sse2/bdiv_q_1.asm

   1 dnl  Intel Pentium-4 mpn_divexact_1 -- mpn by limb exact division.
   2
   3 dnl  Copyright 2001, 2002, 2007, 2011 Free Software Foundation, Inc.
   4 dnl
   5 dnl  This file is part of the GNU MP Library.
   6 dnl
   7 dnl  Rearranged from mpn/x86/pentium4/sse2/dive_1.asm by Marco Bodrato.
   8 dnl
   9 dnl  The GNU MP Library is free software; you can redistribute it and/or
  10 dnl  modify it under the terms of the GNU Lesser General Public License as
  11 dnl  published by the Free Software Foundation; either version 3 of the
  12 dnl  License, or (at your option) any later version.
  13 dnl
  14 dnl  The GNU MP Library is distributed in the hope that it will be useful,
  15 dnl  but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 dnl  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17 dnl  Lesser General Public License for more details.
  18 dnl
  19 dnl  You should have received a copy of the GNU Lesser General Public License
  20 dnl  along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.
  21
  22 include(`../config.m4')
  23
  24
  25 C P4: 19.0 cycles/limb
  26
  27 C Pairs of movd's are used to avoid unaligned loads.  Despite the loads not
  28 C being on the dependent chain and there being plenty of cycles available,
  29 C using an unaligned movq on every second iteration measured about 23 c/l.
  30 C
  31
  32 defframe(PARAM_SHIFT,  24)
  33 defframe(PARAM_INVERSE,20)
  34 defframe(PARAM_DIVISOR,16)
  35 defframe(PARAM_SIZE,   12)
  36 defframe(PARAM_SRC,    8)
  37 defframe(PARAM_DST,    4)
  38
  39         TEXT
  40
  41 C mp_limb_t
  42 C mpn_pi1_bdiv_q_1 (mp_ptr dst, mp_srcptr src, mp_size_t size, mp_limb_t divisor,
  43 C                   mp_limb_t inverse, int shift)
  44         ALIGN(32)
  45 PROLOGUE(mpn_pi1_bdiv_q_1)
  46 deflit(`FRAME',0)
  47
  48         movl    PARAM_SIZE, %edx
  49
  50         movl    PARAM_SRC, %eax
  51
  52         movl    PARAM_DIVISOR, %ecx
  53
  54         movd    %ecx, %mm6
  55         movl    PARAM_SHIFT, %ecx
  56
  57         movd    %ecx, %mm7              C shift
  58
  59         C
  60
  61         movl    PARAM_INVERSE, %ecx
  62         movd    %ecx, %mm5              C inv
  63
  64         movl    PARAM_DST, %ecx
  65         pxor    %mm1, %mm1              C initial carry limb
  66         pxor    %mm0, %mm0              C initial carry bit
  67
  68         subl    $1, %edx
  69         jz      L(done)
  70
  71         pcmpeqd %mm4, %mm4
  72         psrlq   $32, %mm4               C 0x00000000FFFFFFFF
  73
  74 C The dependent chain here is as follows.
  75 C
  76 C                                       latency
  77 C       psubq    s = (src-cbit) - climb    2
  78 C       pmuludq  q = s*inverse             8
  79 C       pmuludq  prod = q*divisor          8
  80 C       psrlq    climb = high(prod)        2
  81 C                                         --
  82 C                                         20
  83 C
  84 C Yet the loop measures 19.0 c/l, so obviously there's something gained
  85 C there over a straight reading of the chip documentation.
  86
  87 L(top):
  88         C eax   src, incrementing
  89         C ebx
  90         C ecx   dst, incrementing
  91         C edx   counter, size-1 iterations
  92         C
  93         C mm0   carry bit
  94         C mm1   carry limb
  95         C mm4   0x00000000FFFFFFFF
  96         C mm5   inverse
  97         C mm6   divisor
  98         C mm7   shift
  99
 100         movd    (%eax), %mm2
 101         movd    4(%eax), %mm3
 102         addl    $4, %eax
 103         punpckldq %mm3, %mm2
 104
 105         psrlq   %mm7, %mm2
 106         pand    %mm4, %mm2              C src
 107         psubq   %mm0, %mm2              C src - cbit
 108
 109         psubq   %mm1, %mm2              C src - cbit - climb
 110         movq    %mm2, %mm0
 111         psrlq   $63, %mm0               C new cbit
 112
 113         pmuludq %mm5, %mm2              C s*inverse
 114         movd    %mm2, (%ecx)            C q
 115         addl    $4, %ecx
 116
 117         movq    %mm6, %mm1
 118         pmuludq %mm2, %mm1              C q*divisor
 119         psrlq   $32, %mm1               C new climb
 120
 121 L(entry):
 122         subl    $1, %edx
 123         jnz     L(top)
 124
 125 L(done):
 126         movd    (%eax), %mm2
 127         psrlq   %mm7, %mm2              C src
 128         psubq   %mm0, %mm2              C src - cbit
 129
 130         psubq   %mm1, %mm2              C src - cbit - climb
 131
 132         pmuludq %mm5, %mm2              C s*inverse
 133         movd    %mm2, (%ecx)            C q
 134
 135         emms
 136         ret
 137
 138 EPILOGUE()
 139
 140         ALIGN(16)
 141 C mp_limb_t mpn_bdiv_q_1 (mp_ptr dst, mp_srcptr src, mp_size_t size,
 142 C                           mp_limb_t divisor);
 143 C
 144 PROLOGUE(mpn_bdiv_q_1)
 145 deflit(`FRAME',0)
 146
 147         movl    PARAM_SIZE, %edx
 148
 149         movl    PARAM_DIVISOR, %ecx
 150
 151         C eax   src
 152         C ebx
 153         C ecx   divisor
 154         C edx   size-1
 155
 156         movl    %ecx, %eax
 157         bsfl    %ecx, %ecx              C trailing twos
 158
 159         shrl    %cl, %eax               C d = divisor without twos
 160         movd    %eax, %mm6
 161         movd    %ecx, %mm7              C shift
 162
 163         shrl    %eax                    C d/2
 164
 165         andl    $127, %eax              C d/2, 7 bits
 166
 167 ifdef(`PIC',`
 168         LEA(    binvert_limb_table, %ecx)
 169         movzbl  (%eax,%ecx), %eax               C inv 8 bits
 170 ',`
 171         movzbl  binvert_limb_table(%eax), %eax  C inv 8 bits
 172 ')
 173
 174         C
 175
 176         movd    %eax, %mm5              C inv
 177
 178         movd    %eax, %mm0              C inv
 179
 180         pmuludq %mm5, %mm5              C inv*inv
 181
 182         C
 183
 184         pmuludq %mm6, %mm5              C inv*inv*d
 185         paddd   %mm0, %mm0              C 2*inv
 186
 187         C
 188
 189         psubd   %mm5, %mm0              C inv = 2*inv - inv*inv*d
 190         pxor    %mm5, %mm5
 191
 192         paddd   %mm0, %mm5
 193         pmuludq %mm0, %mm0              C inv*inv
 194
 195         pcmpeqd %mm4, %mm4
 196         psrlq   $32, %mm4               C 0x00000000FFFFFFFF
 197
 198         C
 199
 200         pmuludq %mm6, %mm0              C inv*inv*d
 201         paddd   %mm5, %mm5              C 2*inv
 202
 203         movl    PARAM_SRC, %eax
 204         movl    PARAM_DST, %ecx
 205         pxor    %mm1, %mm1              C initial carry limb
 206
 207         C
 208
 209         psubd   %mm0, %mm5              C inv = 2*inv - inv*inv*d
 210
 211         ASSERT(e,`      C expect d*inv == 1 mod 2^GMP_LIMB_BITS
 212         pushl   %eax    FRAME_pushl()
 213         movq    %mm6, %mm0
 214         pmuludq %mm5, %mm0
 215         movd    %mm0, %eax
 216         cmpl    $1, %eax
 217         popl    %eax    FRAME_popl()')
 218
 219         pxor    %mm0, %mm0              C initial carry bit
 220         jmp     L(entry)
 221
 222 EPILOGUE()