-
Notifications
You must be signed in to change notification settings - Fork 6
/
divfixedv64.asm
496 lines (443 loc) · 17.9 KB
/
divfixedv64.asm
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
;************************* divfixedv64.asm *********************************
; Author: Agner Fog
; Date created: 2011-07-25
; Last modified: 2012-03-10
;
; Function prototypes:
; void setdivisorV8i16(__m128i buf[2], int16_t d);
; void setdivisorV8u16(__m128i buf[2], uint16_t d);
; void setdivisorV4i32(__m128i buf[2], int32_t d);
; void setdivisorV4u32(__m128i buf[2], uint32_t d);
;
; __m128i dividefixedV8i16(const __m128i buf[2], __m128i x);
; __m128i dividefixedV8u16(const __m128i buf[2], __m128i x);;
; __m128i dividefixedV4i32(const __m128i buf[2], __m128i x);
; __m128i dividefixedV4u32(const __m128i buf[2], __m128i x);
;
; Alternative versions for VectorClass.h:
; (These versions pack all parameters into a single register)
; __m128i setdivisor8s(int16_t d);
; __m128i setdivisor8us(uint16_t d);
; __m128i setdivisor4i(int32_t d);
; __m128i setdivisor4ui(uint32_t d);
;
; Description:
; Functions for integer vector division by the same divisor, signed
; and unsigned 16-bit and 32-bit integer versions.
;
; The setdivisor functions calculate the reciprocal divisor and shift counts,
; the dividefixed functions do the division by multiplication and shift of the
; vector elements of packed 16-bit or 32-bit signed or unsigned integers.
;
; The divisor must be positive. A zero divisor generated a divide by zero error.
; A negative divisor generates a division overflow error. To divide by a negative
; divisor, change the sign of the divisor and the result.
;
; The methods used are described in this article:
; T. Granlund and P. L. Montgomery: Division by Invariant Integers Using Multiplication,
; Proceedings of the SIGPLAN 1994 Conference on Programming Language Design and Implementation.
; http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.1.2556
;
; Mathematical formula, unsigned division:
; x = dividend
; d = divisor
; n = integer size, bits
; L = ceil(log2(d))
; m = 1 + 2^n * (2^L-d) / d [2^L should overflow to 0 if L = n]
; sh1 = min(L,1)
; sh2 = max(L-1,0)
; t = m*x >> n [high part of unsigned multiplication]
; x/d = (((x-t) >> sh1) + t) >> sh2
;
; Mathematical formula, signed division:
; x = dividend
; d = abs(divisor)
; n = integer size, bits
; L = ceil(log2(d))
; L = max(L,1)
; m = 1 + 2^(n+L-1)/d - 2^n [division should overflow to 0 if d = 1]
; sh1 = L-1
; q = x + (m*x >> n) [high part of signed multiplication]
; q = (q >> sh1) - (x<0 ? -1 : 0)
; if (divisor < 0) q = -q [negative divisor not supported in present implementation]
; x/d = q
;
; Copyright (c) 2011 - 2012 GNU General Public License www.gnu.org/licenses
;******************************************************************************
default rel
%IFDEF WINDOWS
%define par1 rcx ; function parameter 1
%define par1d ecx
%define par1w cx
%define par2 rdx ; function parameter 2
%define par2d edx
%define par2w dx
%define buf r8 ; pointer to buffer
%ENDIF
%IFDEF UNIX
%define par1 rdi ; function parameter 1
%define par1d edi
%define par1w di
%define par2 rsi ; function parameter 2
%define par2d esi
%define par2w si
%define buf rdi ; pointer to buffer
%ENDIF
; Imported from instrset64.asm:
extern InstructionSet ; Instruction set for CPU dispatcher
section .text align = 16
;******************************************************************************
; 16 bit signed integers
;******************************************************************************
; extern "C" __m128i setdivisor8s(int16_t d);
; vector of 8 x 16 bit signed integers
global setdivisor8s: function
setdivisor8s:
push rbx
movsx ebx, par1w ; d
dec ebx
mov ecx, -1 ; value for bsr if ebx = 0
bsr ecx, ebx ; floor(log2(d-1))
inc ebx
js H120 ; Generate error if d < 0. (error for d=0 will come in the div instruction)
inc ecx ; L = ceil(log2(d))
sub ecx, 1 ; shift count = L - 1
adc ecx, 0 ; avoid negative shift count
xor eax, eax
mov edx, 1
cmp ebx, edx
je H110 ; avoid division overflow when d = 1
shl edx, cl
div bx ; 2^(16+L-1)/d
H110: inc eax
movd xmm0, eax ; multiplier
pshuflw xmm0, xmm0, 0 ; broadcast into lower 4 words
movd xmm1, ecx ; shift count
punpcklqdq xmm0, xmm1 ; insert shift count into upper half
pop rbx
ret
H120: ; d < 0 not supported. Generate error
mov edx, 1
div edx
ud2
; setdivisor8s end
; extern "C" void setdivisorV8i16(__m128i buf[2], int16_t d);
; vector of 8 x 16 bit signed integers
global setdivisorV8i16: function
setdivisorV8i16:
push par1 ; buf
mov par1d, par2d ; d
call setdivisor8s
pop rax ; buf
punpcklqdq xmm0, xmm0 ; copy multiplier into upper 4 words
movdqa [rax], xmm0 ; multiplier
movdqa [rax+16], xmm1 ; shift count is still in xmm1
ret
; setdivisorV8i16 end
; extern "C" int dividefixedV8i16(const __m128i buf[2], __m128i x);
global dividefixedV8i16: function
dividefixedV8i16:
; buf = par1
; x = xmm0 (UNIX) or [par2] (Windows)
%IFDEF WINDOWS
movdqa xmm0, [par2] ; x
%ENDIF
movdqa xmm1, xmm0 ; x
pmulhw xmm0, [par1] ; multiply high signed words
paddw xmm0, xmm1
movd xmm2, [par1+16] ; shift count
psraw xmm0, xmm2 ; shift right arithmetic
psraw xmm1, 15 ; sign of x
psubw xmm0, xmm1
ret
;dividefixedV8i16 end
;******************************************************************************
; 16 bit unsigned integers
;******************************************************************************
; extern "C" __m128i setdivisor8us(uint16_t d);
; vector of 8 x 16 bit unsigned integers
align 16
global setdivisor8us: function
setdivisor8us:
push rbx
movzx ebx, par1w ; d
dec ebx
mov ecx, -1 ; value for bsr if ebx = 0
bsr ecx, ebx ; floor(log2(d-1))
inc ebx
inc ecx ; L = ceil(log2(d))
mov edx, 1
shl edx, cl ; 2^L [32-bit shift to allow overflow]
sub edx, ebx
xor eax, eax
div bx
inc eax
movd xmm0, eax
pshuflw xmm0, xmm0, 0 ; broadcast into lower 4 words
sub ecx, 1
setae dl
movzx edx, dl ; shift 1
seta al
neg al
and al,cl
movzx eax, al ; shift 2
movd xmm1, edx ; shift 1
movd xmm2, eax ; shift 2
punpckldq xmm1, xmm2 ; combine into two dwords
punpcklqdq xmm0, xmm1 ; multipliers, shift1, shift2
pop rbx
ret
; setdivisor8us end
;extern "C" void setdivisorV8u16(__m128i buf[2], uint16_t d);
; 8 x 16 bit unsigned
global setdivisorV8u16: function
setdivisorV8u16:
push par1 ; buf
mov par1d, par2d ; d
call setdivisor8us
pop rax ; buf
punpcklqdq xmm0, xmm0 ; copy multiplier into upper 4 words
movdqa [rax], xmm0 ; multiplier
movdqa [rax+16], xmm1 ; shift counts are still in xmm1
ret
; setdivisorV8u16 end
;extern "C" __m128i dividefixedV8u16(const __m128i buf[2], __m128i x);
global dividefixedV8u16: function
align 16
dividefixedV8u16:
; buf = par1
; x = xmm0 (UNIX) or [par2] (Windows)
%IFDEF WINDOWS
movdqa xmm0, [par2] ; x
%ENDIF
movdqa xmm1, xmm0 ; x
pmulhuw xmm0, [par1] ; multiply high unsigned words
psubw xmm1, xmm0
movd xmm2, [par1+16] ; shift1
psrlw xmm1, xmm2
paddw xmm0, xmm1
movd xmm2, [par1+20] ; shift2
psrlw xmm0, xmm2
ret
;dividefixedV8u16 end
;******************************************************************************
; 32 bit signed integers
;******************************************************************************
; extern "C" __m128i setdivisor4i(int32_t d);
; vector of 4 x 32 bit signed integers
align 16
global setdivisor4i: function
setdivisor4i:
push rbx
mov ebx, par1d ; d
dec ebx
mov ecx, -1 ; value for bsr if ebx = 0
bsr ecx, ebx ; floor(log2(d-1))
inc ebx
js K120 ; Generate error if d < 0. (error for d=0 will come in the div instruction)
inc ecx ; L = ceil(log2(d))
sub ecx, 1 ; shift count = L - 1
adc ecx, 0 ; avoid negative shift count
xor eax, eax
mov edx, 1
cmp ebx, edx
je K110 ; avoid division overflow when d = 1
shl edx, cl
div ebx ; 2^(16+L-1)/d
K110: inc eax
movd xmm0, eax ; multiplier
pshufd xmm0, xmm0, 0 ; broadcast into 4 dwords
movd xmm1, ecx ; shift count
punpcklqdq xmm0, xmm1 ; insert shift count into upper half
pop rbx
ret
K120: ; d < 0 not supported. Generate error
mov edx, 1
div edx
ud2
; setdivisor4i end
; extern "C" void setdivisorV4i32(__m128i buf[2], int32_t d);
; vector of 4 x 32 bit signed integers
global setdivisorV4i32: function
setdivisorV4i32:
push par1 ; buf
mov par1d, par2d ; d
call setdivisor4i
pop rax ; buf
punpcklqdq xmm0, xmm0 ; copy multiplier into upper 4 words
movdqa [rax], xmm0 ; multiplier
movdqa [rax+16], xmm1 ; shift count is still in xmm1
ret
; setdivisorV4i32 end
; extern "C" int dividefixedV4i32(const __m128i buf[2], __m128i x);
global dividefixedV4i32: function
; Direct entries to CPU-specific versions
global dividefixedV4i32SSE2: function
global dividefixedV4i32SSE41: function
align 8
dividefixedV4i32: ; function dispatching
jmp near [dividefixedV4i32Dispatch] ; Go to appropriate version, depending on instruction set
align 16
dividefixedV4i32SSE41:
; buf = par1
; x = xmm0 (UNIX) or [par2] (Windows)
%IFDEF WINDOWS
movdqa xmm0,[par2] ; x
%ENDIF
movdqa xmm1, xmm0 ; x
movdqa xmm2, xmm0 ; x
movdqa xmm3, [par1] ; multiplier
pmuldq xmm0, xmm3 ; 32 x 32 -> 64 bit signed multiplication of x[0] and x[2]
psrlq xmm0, 32 ; high dword of result 0 and 2
psrlq xmm1, 32 ; get x[1] and x[3] into position for multiplication
pmuldq xmm1, xmm3 ; 32 x 32 -> 64 bit signed multiplication of x[1] and x[3]
pcmpeqd xmm3, xmm3
psllq xmm3, 32 ; generate mask of dword 1 and 3
pand xmm1, xmm3 ; high dword of result 1 and 3
por xmm0, xmm1 ; combine all four results into one vector
paddd xmm0, xmm2
movd xmm3, [par1+16] ; shift count
psrad xmm0, xmm3 ; shift right arithmetic
psrad xmm2, 31 ; sign of x
psubd xmm0, xmm2
ret
;dividefixedV4i32SSE41 end
dividefixedV4i32SSE2:
; I have tried to change sign and use pmuludq, but get rounding error (gives 9/10 = 1).
; This solution, with 4 separate multiplications, is probably faster anyway despite store forwarding stall
push rbp
mov rbp, rsp
%IFDEF WINDOWS
movdqa xmm0,[par2] ; x
mov buf, par1
%ENDIF
sub rsp, 16 ; allocate stack space
and rsp, -16 ; stack should be aligned already. align anyway to be safe
movdqa [rsp], xmm0 ; store x
movdqa xmm2, xmm0 ; x
mov ecx, [buf] ; multiplier
; do four signed high multiplications
mov eax, [rsp]
imul ecx
mov [rsp], edx
mov eax, [rsp+4]
imul ecx
mov [rsp+4], edx
mov eax, [rsp+8]
imul ecx
mov [rsp+8], edx
mov eax, [rsp+12]
imul ecx
mov [rsp+12], edx
movdqa xmm0, [rsp] ; x*m vector
paddd xmm0, xmm2
movd xmm3, [buf+16] ; shift count
psrad xmm0, xmm3 ; shift right arithmetic
psrad xmm2, 31 ; sign of x
psubd xmm0, xmm2
mov rsp, rbp
pop rbp
ret
;dividefixedV4i32SSE2 end
; ********************************************************************************
; CPU dispatching for dividefixedV4i32. This is executed only once
; ********************************************************************************
dividefixedV4i32CPUDispatch:
; get supported instruction set
push par1
push par2
call InstructionSet
pop par2
pop par1
; Point to generic version
lea r8, [dividefixedV4i32SSE2]
cmp eax, 8 ; check if PMULDQ supported
jb Q100
; SSE4.1 supported
; Point to SSE4.1 version of strstr
lea r8, [dividefixedV4i32SSE41]
Q100: mov [dividefixedV4i32Dispatch], r8
; Continue in appropriate version
jmp r8
SECTION .data
; Pointer to appropriate versions. Initially point to dispatcher
dividefixedV4i32Dispatch Dq dividefixedV4i32CPUDispatch
section .text
;******************************************************************************
; 32 bit unsigned integers
;******************************************************************************
; extern "C" __m128i setdivisor4ui(uint32_t d);
; vector of 4 x 32 bit unsigned integers
align 16
global setdivisor4ui: function
setdivisor4ui:
push rbx
mov ebx, par1d ; d
dec ebx
mov ecx, -1 ; value for bsr if ebx = 0
bsr ecx, ebx ; floor(log2(d-1))
inc ebx
inc ecx ; L = ceil(log2(d))
mov edx, 1
shl rdx, cl ; 2^L [64 bit shift to allow overflow]
sub edx, ebx
xor eax, eax
div ebx
inc eax
movd xmm0, eax
pshufd xmm0, xmm0, 0 ; broadcast into 4 dwords
sub ecx, 1
setae dl
movzx edx, dl ; shift1
seta al
neg al
and al,cl
movzx eax, al
movd xmm1, edx ; shift 1
movd xmm2, eax ; shift 2
punpckldq xmm1, xmm2 ; combine into two dwords
punpcklqdq xmm0, xmm1 ; multipliers, shift1, shift2
pop rbx
ret
; setdivisor4ui end
;extern "C" void setdivisorV4u32(__m128i buf[2], uint32_t d);
; 4 x 32 bit unsigned
global setdivisorV4u32: function
setdivisorV4u32:
push par1 ; buf
mov par1d, par2d ; d
call setdivisor4ui
pop rax ; buf
punpcklqdq xmm0, xmm0 ; copy multiplier into upper 4 words
movdqa [rax], xmm0 ; multiplier
movdqa [rax+16], xmm1 ; shift counts are still in xmm1
ret
; setdivisorV4u32 end
;extern "C" __m128i dividefixedV4u32(const __m128i buf[2], __m128i x);
global dividefixedV4u32: function
align 16
dividefixedV4u32:
; buf = par1
; x = xmm0 (UNIX) or [par2] (Windows)
%IFDEF WINDOWS
movdqa xmm0,[par2] ; x
%ENDIF
movdqa xmm1, xmm0 ; x
movdqa xmm2, xmm0 ; x
movdqa xmm3, [par1] ; multiplier
pmuludq xmm0, xmm3 ; 32 x 32 -> 64 bit unsigned multiplication of x[0] and x[2]
psrlq xmm0, 32 ; high dword of result 0 and 2
psrlq xmm1, 32 ; get x[1] and x[3] into position for multiplication
pmuludq xmm1, xmm3 ; 32 x 32 -> 64 bit unsigned multiplication of x[1] and x[3]
pcmpeqd xmm3, xmm3
psllq xmm3, 32 ; generate mask of dword 1 and 3
pand xmm1, xmm3 ; high dword of result 1 and 3
por xmm0, xmm1 ; combine all four results into one vector
psubd xmm2, xmm0
movd xmm3, [par1+16] ; shift1
psrld xmm2, xmm3
paddd xmm0, xmm2
movd xmm3, [par1+20] ; shift2
psrld xmm0, xmm3
ret
;dividefixedV4u32 end