xref: /qemu/target/arm/tcg/vec_helper.c (revision 7782a9afec81d1efe23572135c1ed777691ccde5)
1 /*
2  * ARM AdvSIMD / SVE Vector Operations
3  *
4  * Copyright (c) 2018 Linaro
5  *
6  * This library is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Lesser General Public
8  * License as published by the Free Software Foundation; either
9  * version 2 of the License, or (at your option) any later version.
10  *
11  * This library is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this library; if not, see <http://www.gnu.org/licenses/>.
18  */
19 
20 #include "qemu/osdep.h"
21 #include "cpu.h"
22 #include "exec/helper-proto.h"
23 #include "tcg/tcg-gvec-desc.h"
24 #include "fpu/softfloat.h"
25 #include "vec_internal.h"
26 
27 /* Note that vector data is stored in host-endian 64-bit chunks,
28    so addressing units smaller than that needs a host-endian fixup.  */
29 #ifdef HOST_WORDS_BIGENDIAN
30 #define H1(x)  ((x) ^ 7)
31 #define H2(x)  ((x) ^ 3)
32 #define H4(x)  ((x) ^ 1)
33 #else
34 #define H1(x)  (x)
35 #define H2(x)  (x)
36 #define H4(x)  (x)
37 #endif
38 
39 /* Signed saturating rounding doubling multiply-accumulate high half, 16-bit */
40 static int16_t do_sqrdmlah_h(int16_t src1, int16_t src2, int16_t src3,
41                              bool neg, bool round, uint32_t *sat)
42 {
43     /*
44      * Simplify:
45      * = ((a3 << 16) + ((e1 * e2) << 1) + (1 << 15)) >> 16
46      * = ((a3 << 15) + (e1 * e2) + (1 << 14)) >> 15
47      */
48     int32_t ret = (int32_t)src1 * src2;
49     if (neg) {
50         ret = -ret;
51     }
52     ret += ((int32_t)src3 << 15) + (round << 14);
53     ret >>= 15;
54 
55     if (ret != (int16_t)ret) {
56         *sat = 1;
57         ret = (ret < 0 ? INT16_MIN : INT16_MAX);
58     }
59     return ret;
60 }
61 
62 uint32_t HELPER(neon_qrdmlah_s16)(CPUARMState *env, uint32_t src1,
63                                   uint32_t src2, uint32_t src3)
64 {
65     uint32_t *sat = &env->vfp.qc[0];
66     uint16_t e1 = do_sqrdmlah_h(src1, src2, src3, false, true, sat);
67     uint16_t e2 = do_sqrdmlah_h(src1 >> 16, src2 >> 16, src3 >> 16,
68                                 false, true, sat);
69     return deposit32(e1, 16, 16, e2);
70 }
71 
72 void HELPER(gvec_qrdmlah_s16)(void *vd, void *vn, void *vm,
73                               void *vq, uint32_t desc)
74 {
75     uintptr_t opr_sz = simd_oprsz(desc);
76     int16_t *d = vd;
77     int16_t *n = vn;
78     int16_t *m = vm;
79     uintptr_t i;
80 
81     for (i = 0; i < opr_sz / 2; ++i) {
82         d[i] = do_sqrdmlah_h(n[i], m[i], d[i], false, true, vq);
83     }
84     clear_tail(d, opr_sz, simd_maxsz(desc));
85 }
86 
87 uint32_t HELPER(neon_qrdmlsh_s16)(CPUARMState *env, uint32_t src1,
88                                   uint32_t src2, uint32_t src3)
89 {
90     uint32_t *sat = &env->vfp.qc[0];
91     uint16_t e1 = do_sqrdmlah_h(src1, src2, src3, true, true, sat);
92     uint16_t e2 = do_sqrdmlah_h(src1 >> 16, src2 >> 16, src3 >> 16,
93                                 true, true, sat);
94     return deposit32(e1, 16, 16, e2);
95 }
96 
97 void HELPER(gvec_qrdmlsh_s16)(void *vd, void *vn, void *vm,
98                               void *vq, uint32_t desc)
99 {
100     uintptr_t opr_sz = simd_oprsz(desc);
101     int16_t *d = vd;
102     int16_t *n = vn;
103     int16_t *m = vm;
104     uintptr_t i;
105 
106     for (i = 0; i < opr_sz / 2; ++i) {
107         d[i] = do_sqrdmlah_h(n[i], m[i], d[i], true, true, vq);
108     }
109     clear_tail(d, opr_sz, simd_maxsz(desc));
110 }
111 
112 void HELPER(neon_sqdmulh_h)(void *vd, void *vn, void *vm,
113                             void *vq, uint32_t desc)
114 {
115     intptr_t i, opr_sz = simd_oprsz(desc);
116     int16_t *d = vd, *n = vn, *m = vm;
117 
118     for (i = 0; i < opr_sz / 2; ++i) {
119         d[i] = do_sqrdmlah_h(n[i], m[i], 0, false, false, vq);
120     }
121     clear_tail(d, opr_sz, simd_maxsz(desc));
122 }
123 
124 void HELPER(neon_sqrdmulh_h)(void *vd, void *vn, void *vm,
125                              void *vq, uint32_t desc)
126 {
127     intptr_t i, opr_sz = simd_oprsz(desc);
128     int16_t *d = vd, *n = vn, *m = vm;
129 
130     for (i = 0; i < opr_sz / 2; ++i) {
131         d[i] = do_sqrdmlah_h(n[i], m[i], 0, false, true, vq);
132     }
133     clear_tail(d, opr_sz, simd_maxsz(desc));
134 }
135 
136 /* Signed saturating rounding doubling multiply-accumulate high half, 32-bit */
137 static int32_t do_sqrdmlah_s(int32_t src1, int32_t src2, int32_t src3,
138                              bool neg, bool round, uint32_t *sat)
139 {
140     /* Simplify similarly to int_qrdmlah_s16 above.  */
141     int64_t ret = (int64_t)src1 * src2;
142     if (neg) {
143         ret = -ret;
144     }
145     ret += ((int64_t)src3 << 31) + (round << 30);
146     ret >>= 31;
147 
148     if (ret != (int32_t)ret) {
149         *sat = 1;
150         ret = (ret < 0 ? INT32_MIN : INT32_MAX);
151     }
152     return ret;
153 }
154 
155 uint32_t HELPER(neon_qrdmlah_s32)(CPUARMState *env, int32_t src1,
156                                   int32_t src2, int32_t src3)
157 {
158     uint32_t *sat = &env->vfp.qc[0];
159     return do_sqrdmlah_s(src1, src2, src3, false, true, sat);
160 }
161 
162 void HELPER(gvec_qrdmlah_s32)(void *vd, void *vn, void *vm,
163                               void *vq, uint32_t desc)
164 {
165     uintptr_t opr_sz = simd_oprsz(desc);
166     int32_t *d = vd;
167     int32_t *n = vn;
168     int32_t *m = vm;
169     uintptr_t i;
170 
171     for (i = 0; i < opr_sz / 4; ++i) {
172         d[i] = do_sqrdmlah_s(n[i], m[i], d[i], false, true, vq);
173     }
174     clear_tail(d, opr_sz, simd_maxsz(desc));
175 }
176 
177 uint32_t HELPER(neon_qrdmlsh_s32)(CPUARMState *env, int32_t src1,
178                                   int32_t src2, int32_t src3)
179 {
180     uint32_t *sat = &env->vfp.qc[0];
181     return do_sqrdmlah_s(src1, src2, src3, true, true, sat);
182 }
183 
184 void HELPER(gvec_qrdmlsh_s32)(void *vd, void *vn, void *vm,
185                               void *vq, uint32_t desc)
186 {
187     uintptr_t opr_sz = simd_oprsz(desc);
188     int32_t *d = vd;
189     int32_t *n = vn;
190     int32_t *m = vm;
191     uintptr_t i;
192 
193     for (i = 0; i < opr_sz / 4; ++i) {
194         d[i] = do_sqrdmlah_s(n[i], m[i], d[i], true, true, vq);
195     }
196     clear_tail(d, opr_sz, simd_maxsz(desc));
197 }
198 
199 void HELPER(neon_sqdmulh_s)(void *vd, void *vn, void *vm,
200                             void *vq, uint32_t desc)
201 {
202     intptr_t i, opr_sz = simd_oprsz(desc);
203     int32_t *d = vd, *n = vn, *m = vm;
204 
205     for (i = 0; i < opr_sz / 4; ++i) {
206         d[i] = do_sqrdmlah_s(n[i], m[i], 0, false, false, vq);
207     }
208     clear_tail(d, opr_sz, simd_maxsz(desc));
209 }
210 
211 void HELPER(neon_sqrdmulh_s)(void *vd, void *vn, void *vm,
212                              void *vq, uint32_t desc)
213 {
214     intptr_t i, opr_sz = simd_oprsz(desc);
215     int32_t *d = vd, *n = vn, *m = vm;
216 
217     for (i = 0; i < opr_sz / 4; ++i) {
218         d[i] = do_sqrdmlah_s(n[i], m[i], 0, false, true, vq);
219     }
220     clear_tail(d, opr_sz, simd_maxsz(desc));
221 }
222 
223 /* Integer 8 and 16-bit dot-product.
224  *
225  * Note that for the loops herein, host endianness does not matter
226  * with respect to the ordering of data within the 64-bit lanes.
227  * All elements are treated equally, no matter where they are.
228  */
229 
230 void HELPER(gvec_sdot_b)(void *vd, void *vn, void *vm, uint32_t desc)
231 {
232     intptr_t i, opr_sz = simd_oprsz(desc);
233     uint32_t *d = vd;
234     int8_t *n = vn, *m = vm;
235 
236     for (i = 0; i < opr_sz / 4; ++i) {
237         d[i] += n[i * 4 + 0] * m[i * 4 + 0]
238               + n[i * 4 + 1] * m[i * 4 + 1]
239               + n[i * 4 + 2] * m[i * 4 + 2]
240               + n[i * 4 + 3] * m[i * 4 + 3];
241     }
242     clear_tail(d, opr_sz, simd_maxsz(desc));
243 }
244 
245 void HELPER(gvec_udot_b)(void *vd, void *vn, void *vm, uint32_t desc)
246 {
247     intptr_t i, opr_sz = simd_oprsz(desc);
248     uint32_t *d = vd;
249     uint8_t *n = vn, *m = vm;
250 
251     for (i = 0; i < opr_sz / 4; ++i) {
252         d[i] += n[i * 4 + 0] * m[i * 4 + 0]
253               + n[i * 4 + 1] * m[i * 4 + 1]
254               + n[i * 4 + 2] * m[i * 4 + 2]
255               + n[i * 4 + 3] * m[i * 4 + 3];
256     }
257     clear_tail(d, opr_sz, simd_maxsz(desc));
258 }
259 
260 void HELPER(gvec_sdot_h)(void *vd, void *vn, void *vm, uint32_t desc)
261 {
262     intptr_t i, opr_sz = simd_oprsz(desc);
263     uint64_t *d = vd;
264     int16_t *n = vn, *m = vm;
265 
266     for (i = 0; i < opr_sz / 8; ++i) {
267         d[i] += (int64_t)n[i * 4 + 0] * m[i * 4 + 0]
268               + (int64_t)n[i * 4 + 1] * m[i * 4 + 1]
269               + (int64_t)n[i * 4 + 2] * m[i * 4 + 2]
270               + (int64_t)n[i * 4 + 3] * m[i * 4 + 3];
271     }
272     clear_tail(d, opr_sz, simd_maxsz(desc));
273 }
274 
275 void HELPER(gvec_udot_h)(void *vd, void *vn, void *vm, uint32_t desc)
276 {
277     intptr_t i, opr_sz = simd_oprsz(desc);
278     uint64_t *d = vd;
279     uint16_t *n = vn, *m = vm;
280 
281     for (i = 0; i < opr_sz / 8; ++i) {
282         d[i] += (uint64_t)n[i * 4 + 0] * m[i * 4 + 0]
283               + (uint64_t)n[i * 4 + 1] * m[i * 4 + 1]
284               + (uint64_t)n[i * 4 + 2] * m[i * 4 + 2]
285               + (uint64_t)n[i * 4 + 3] * m[i * 4 + 3];
286     }
287     clear_tail(d, opr_sz, simd_maxsz(desc));
288 }
289 
290 void HELPER(gvec_sdot_idx_b)(void *vd, void *vn, void *vm, uint32_t desc)
291 {
292     intptr_t i, segend, opr_sz = simd_oprsz(desc), opr_sz_4 = opr_sz / 4;
293     intptr_t index = simd_data(desc);
294     uint32_t *d = vd;
295     int8_t *n = vn;
296     int8_t *m_indexed = (int8_t *)vm + index * 4;
297 
298     /* Notice the special case of opr_sz == 8, from aa64/aa32 advsimd.
299      * Otherwise opr_sz is a multiple of 16.
300      */
301     segend = MIN(4, opr_sz_4);
302     i = 0;
303     do {
304         int8_t m0 = m_indexed[i * 4 + 0];
305         int8_t m1 = m_indexed[i * 4 + 1];
306         int8_t m2 = m_indexed[i * 4 + 2];
307         int8_t m3 = m_indexed[i * 4 + 3];
308 
309         do {
310             d[i] += n[i * 4 + 0] * m0
311                   + n[i * 4 + 1] * m1
312                   + n[i * 4 + 2] * m2
313                   + n[i * 4 + 3] * m3;
314         } while (++i < segend);
315         segend = i + 4;
316     } while (i < opr_sz_4);
317 
318     clear_tail(d, opr_sz, simd_maxsz(desc));
319 }
320 
321 void HELPER(gvec_udot_idx_b)(void *vd, void *vn, void *vm, uint32_t desc)
322 {
323     intptr_t i, segend, opr_sz = simd_oprsz(desc), opr_sz_4 = opr_sz / 4;
324     intptr_t index = simd_data(desc);
325     uint32_t *d = vd;
326     uint8_t *n = vn;
327     uint8_t *m_indexed = (uint8_t *)vm + index * 4;
328 
329     /* Notice the special case of opr_sz == 8, from aa64/aa32 advsimd.
330      * Otherwise opr_sz is a multiple of 16.
331      */
332     segend = MIN(4, opr_sz_4);
333     i = 0;
334     do {
335         uint8_t m0 = m_indexed[i * 4 + 0];
336         uint8_t m1 = m_indexed[i * 4 + 1];
337         uint8_t m2 = m_indexed[i * 4 + 2];
338         uint8_t m3 = m_indexed[i * 4 + 3];
339 
340         do {
341             d[i] += n[i * 4 + 0] * m0
342                   + n[i * 4 + 1] * m1
343                   + n[i * 4 + 2] * m2
344                   + n[i * 4 + 3] * m3;
345         } while (++i < segend);
346         segend = i + 4;
347     } while (i < opr_sz_4);
348 
349     clear_tail(d, opr_sz, simd_maxsz(desc));
350 }
351 
352 void HELPER(gvec_sdot_idx_h)(void *vd, void *vn, void *vm, uint32_t desc)
353 {
354     intptr_t i, opr_sz = simd_oprsz(desc), opr_sz_8 = opr_sz / 8;
355     intptr_t index = simd_data(desc);
356     uint64_t *d = vd;
357     int16_t *n = vn;
358     int16_t *m_indexed = (int16_t *)vm + index * 4;
359 
360     /* This is supported by SVE only, so opr_sz is always a multiple of 16.
361      * Process the entire segment all at once, writing back the results
362      * only after we've consumed all of the inputs.
363      */
364     for (i = 0; i < opr_sz_8 ; i += 2) {
365         uint64_t d0, d1;
366 
367         d0  = n[i * 4 + 0] * (int64_t)m_indexed[i * 4 + 0];
368         d0 += n[i * 4 + 1] * (int64_t)m_indexed[i * 4 + 1];
369         d0 += n[i * 4 + 2] * (int64_t)m_indexed[i * 4 + 2];
370         d0 += n[i * 4 + 3] * (int64_t)m_indexed[i * 4 + 3];
371         d1  = n[i * 4 + 4] * (int64_t)m_indexed[i * 4 + 0];
372         d1 += n[i * 4 + 5] * (int64_t)m_indexed[i * 4 + 1];
373         d1 += n[i * 4 + 6] * (int64_t)m_indexed[i * 4 + 2];
374         d1 += n[i * 4 + 7] * (int64_t)m_indexed[i * 4 + 3];
375 
376         d[i + 0] += d0;
377         d[i + 1] += d1;
378     }
379 
380     clear_tail(d, opr_sz, simd_maxsz(desc));
381 }
382 
383 void HELPER(gvec_udot_idx_h)(void *vd, void *vn, void *vm, uint32_t desc)
384 {
385     intptr_t i, opr_sz = simd_oprsz(desc), opr_sz_8 = opr_sz / 8;
386     intptr_t index = simd_data(desc);
387     uint64_t *d = vd;
388     uint16_t *n = vn;
389     uint16_t *m_indexed = (uint16_t *)vm + index * 4;
390 
391     /* This is supported by SVE only, so opr_sz is always a multiple of 16.
392      * Process the entire segment all at once, writing back the results
393      * only after we've consumed all of the inputs.
394      */
395     for (i = 0; i < opr_sz_8 ; i += 2) {
396         uint64_t d0, d1;
397 
398         d0  = n[i * 4 + 0] * (uint64_t)m_indexed[i * 4 + 0];
399         d0 += n[i * 4 + 1] * (uint64_t)m_indexed[i * 4 + 1];
400         d0 += n[i * 4 + 2] * (uint64_t)m_indexed[i * 4 + 2];
401         d0 += n[i * 4 + 3] * (uint64_t)m_indexed[i * 4 + 3];
402         d1  = n[i * 4 + 4] * (uint64_t)m_indexed[i * 4 + 0];
403         d1 += n[i * 4 + 5] * (uint64_t)m_indexed[i * 4 + 1];
404         d1 += n[i * 4 + 6] * (uint64_t)m_indexed[i * 4 + 2];
405         d1 += n[i * 4 + 7] * (uint64_t)m_indexed[i * 4 + 3];
406 
407         d[i + 0] += d0;
408         d[i + 1] += d1;
409     }
410 
411     clear_tail(d, opr_sz, simd_maxsz(desc));
412 }
413 
414 void HELPER(gvec_fcaddh)(void *vd, void *vn, void *vm,
415                          void *vfpst, uint32_t desc)
416 {
417     uintptr_t opr_sz = simd_oprsz(desc);
418     float16 *d = vd;
419     float16 *n = vn;
420     float16 *m = vm;
421     float_status *fpst = vfpst;
422     uint32_t neg_real = extract32(desc, SIMD_DATA_SHIFT, 1);
423     uint32_t neg_imag = neg_real ^ 1;
424     uintptr_t i;
425 
426     /* Shift boolean to the sign bit so we can xor to negate.  */
427     neg_real <<= 15;
428     neg_imag <<= 15;
429 
430     for (i = 0; i < opr_sz / 2; i += 2) {
431         float16 e0 = n[H2(i)];
432         float16 e1 = m[H2(i + 1)] ^ neg_imag;
433         float16 e2 = n[H2(i + 1)];
434         float16 e3 = m[H2(i)] ^ neg_real;
435 
436         d[H2(i)] = float16_add(e0, e1, fpst);
437         d[H2(i + 1)] = float16_add(e2, e3, fpst);
438     }
439     clear_tail(d, opr_sz, simd_maxsz(desc));
440 }
441 
442 void HELPER(gvec_fcadds)(void *vd, void *vn, void *vm,
443                          void *vfpst, uint32_t desc)
444 {
445     uintptr_t opr_sz = simd_oprsz(desc);
446     float32 *d = vd;
447     float32 *n = vn;
448     float32 *m = vm;
449     float_status *fpst = vfpst;
450     uint32_t neg_real = extract32(desc, SIMD_DATA_SHIFT, 1);
451     uint32_t neg_imag = neg_real ^ 1;
452     uintptr_t i;
453 
454     /* Shift boolean to the sign bit so we can xor to negate.  */
455     neg_real <<= 31;
456     neg_imag <<= 31;
457 
458     for (i = 0; i < opr_sz / 4; i += 2) {
459         float32 e0 = n[H4(i)];
460         float32 e1 = m[H4(i + 1)] ^ neg_imag;
461         float32 e2 = n[H4(i + 1)];
462         float32 e3 = m[H4(i)] ^ neg_real;
463 
464         d[H4(i)] = float32_add(e0, e1, fpst);
465         d[H4(i + 1)] = float32_add(e2, e3, fpst);
466     }
467     clear_tail(d, opr_sz, simd_maxsz(desc));
468 }
469 
470 void HELPER(gvec_fcaddd)(void *vd, void *vn, void *vm,
471                          void *vfpst, uint32_t desc)
472 {
473     uintptr_t opr_sz = simd_oprsz(desc);
474     float64 *d = vd;
475     float64 *n = vn;
476     float64 *m = vm;
477     float_status *fpst = vfpst;
478     uint64_t neg_real = extract64(desc, SIMD_DATA_SHIFT, 1);
479     uint64_t neg_imag = neg_real ^ 1;
480     uintptr_t i;
481 
482     /* Shift boolean to the sign bit so we can xor to negate.  */
483     neg_real <<= 63;
484     neg_imag <<= 63;
485 
486     for (i = 0; i < opr_sz / 8; i += 2) {
487         float64 e0 = n[i];
488         float64 e1 = m[i + 1] ^ neg_imag;
489         float64 e2 = n[i + 1];
490         float64 e3 = m[i] ^ neg_real;
491 
492         d[i] = float64_add(e0, e1, fpst);
493         d[i + 1] = float64_add(e2, e3, fpst);
494     }
495     clear_tail(d, opr_sz, simd_maxsz(desc));
496 }
497 
498 void HELPER(gvec_fcmlah)(void *vd, void *vn, void *vm,
499                          void *vfpst, uint32_t desc)
500 {
501     uintptr_t opr_sz = simd_oprsz(desc);
502     float16 *d = vd;
503     float16 *n = vn;
504     float16 *m = vm;
505     float_status *fpst = vfpst;
506     intptr_t flip = extract32(desc, SIMD_DATA_SHIFT, 1);
507     uint32_t neg_imag = extract32(desc, SIMD_DATA_SHIFT + 1, 1);
508     uint32_t neg_real = flip ^ neg_imag;
509     uintptr_t i;
510 
511     /* Shift boolean to the sign bit so we can xor to negate.  */
512     neg_real <<= 15;
513     neg_imag <<= 15;
514 
515     for (i = 0; i < opr_sz / 2; i += 2) {
516         float16 e2 = n[H2(i + flip)];
517         float16 e1 = m[H2(i + flip)] ^ neg_real;
518         float16 e4 = e2;
519         float16 e3 = m[H2(i + 1 - flip)] ^ neg_imag;
520 
521         d[H2(i)] = float16_muladd(e2, e1, d[H2(i)], 0, fpst);
522         d[H2(i + 1)] = float16_muladd(e4, e3, d[H2(i + 1)], 0, fpst);
523     }
524     clear_tail(d, opr_sz, simd_maxsz(desc));
525 }
526 
527 void HELPER(gvec_fcmlah_idx)(void *vd, void *vn, void *vm,
528                              void *vfpst, uint32_t desc)
529 {
530     uintptr_t opr_sz = simd_oprsz(desc);
531     float16 *d = vd;
532     float16 *n = vn;
533     float16 *m = vm;
534     float_status *fpst = vfpst;
535     intptr_t flip = extract32(desc, SIMD_DATA_SHIFT, 1);
536     uint32_t neg_imag = extract32(desc, SIMD_DATA_SHIFT + 1, 1);
537     intptr_t index = extract32(desc, SIMD_DATA_SHIFT + 2, 2);
538     uint32_t neg_real = flip ^ neg_imag;
539     intptr_t elements = opr_sz / sizeof(float16);
540     intptr_t eltspersegment = 16 / sizeof(float16);
541     intptr_t i, j;
542 
543     /* Shift boolean to the sign bit so we can xor to negate.  */
544     neg_real <<= 15;
545     neg_imag <<= 15;
546 
547     for (i = 0; i < elements; i += eltspersegment) {
548         float16 mr = m[H2(i + 2 * index + 0)];
549         float16 mi = m[H2(i + 2 * index + 1)];
550         float16 e1 = neg_real ^ (flip ? mi : mr);
551         float16 e3 = neg_imag ^ (flip ? mr : mi);
552 
553         for (j = i; j < i + eltspersegment; j += 2) {
554             float16 e2 = n[H2(j + flip)];
555             float16 e4 = e2;
556 
557             d[H2(j)] = float16_muladd(e2, e1, d[H2(j)], 0, fpst);
558             d[H2(j + 1)] = float16_muladd(e4, e3, d[H2(j + 1)], 0, fpst);
559         }
560     }
561     clear_tail(d, opr_sz, simd_maxsz(desc));
562 }
563 
564 void HELPER(gvec_fcmlas)(void *vd, void *vn, void *vm,
565                          void *vfpst, uint32_t desc)
566 {
567     uintptr_t opr_sz = simd_oprsz(desc);
568     float32 *d = vd;
569     float32 *n = vn;
570     float32 *m = vm;
571     float_status *fpst = vfpst;
572     intptr_t flip = extract32(desc, SIMD_DATA_SHIFT, 1);
573     uint32_t neg_imag = extract32(desc, SIMD_DATA_SHIFT + 1, 1);
574     uint32_t neg_real = flip ^ neg_imag;
575     uintptr_t i;
576 
577     /* Shift boolean to the sign bit so we can xor to negate.  */
578     neg_real <<= 31;
579     neg_imag <<= 31;
580 
581     for (i = 0; i < opr_sz / 4; i += 2) {
582         float32 e2 = n[H4(i + flip)];
583         float32 e1 = m[H4(i + flip)] ^ neg_real;
584         float32 e4 = e2;
585         float32 e3 = m[H4(i + 1 - flip)] ^ neg_imag;
586 
587         d[H4(i)] = float32_muladd(e2, e1, d[H4(i)], 0, fpst);
588         d[H4(i + 1)] = float32_muladd(e4, e3, d[H4(i + 1)], 0, fpst);
589     }
590     clear_tail(d, opr_sz, simd_maxsz(desc));
591 }
592 
593 void HELPER(gvec_fcmlas_idx)(void *vd, void *vn, void *vm,
594                              void *vfpst, uint32_t desc)
595 {
596     uintptr_t opr_sz = simd_oprsz(desc);
597     float32 *d = vd;
598     float32 *n = vn;
599     float32 *m = vm;
600     float_status *fpst = vfpst;
601     intptr_t flip = extract32(desc, SIMD_DATA_SHIFT, 1);
602     uint32_t neg_imag = extract32(desc, SIMD_DATA_SHIFT + 1, 1);
603     intptr_t index = extract32(desc, SIMD_DATA_SHIFT + 2, 2);
604     uint32_t neg_real = flip ^ neg_imag;
605     intptr_t elements = opr_sz / sizeof(float32);
606     intptr_t eltspersegment = 16 / sizeof(float32);
607     intptr_t i, j;
608 
609     /* Shift boolean to the sign bit so we can xor to negate.  */
610     neg_real <<= 31;
611     neg_imag <<= 31;
612 
613     for (i = 0; i < elements; i += eltspersegment) {
614         float32 mr = m[H4(i + 2 * index + 0)];
615         float32 mi = m[H4(i + 2 * index + 1)];
616         float32 e1 = neg_real ^ (flip ? mi : mr);
617         float32 e3 = neg_imag ^ (flip ? mr : mi);
618 
619         for (j = i; j < i + eltspersegment; j += 2) {
620             float32 e2 = n[H4(j + flip)];
621             float32 e4 = e2;
622 
623             d[H4(j)] = float32_muladd(e2, e1, d[H4(j)], 0, fpst);
624             d[H4(j + 1)] = float32_muladd(e4, e3, d[H4(j + 1)], 0, fpst);
625         }
626     }
627     clear_tail(d, opr_sz, simd_maxsz(desc));
628 }
629 
630 void HELPER(gvec_fcmlad)(void *vd, void *vn, void *vm,
631                          void *vfpst, uint32_t desc)
632 {
633     uintptr_t opr_sz = simd_oprsz(desc);
634     float64 *d = vd;
635     float64 *n = vn;
636     float64 *m = vm;
637     float_status *fpst = vfpst;
638     intptr_t flip = extract32(desc, SIMD_DATA_SHIFT, 1);
639     uint64_t neg_imag = extract32(desc, SIMD_DATA_SHIFT + 1, 1);
640     uint64_t neg_real = flip ^ neg_imag;
641     uintptr_t i;
642 
643     /* Shift boolean to the sign bit so we can xor to negate.  */
644     neg_real <<= 63;
645     neg_imag <<= 63;
646 
647     for (i = 0; i < opr_sz / 8; i += 2) {
648         float64 e2 = n[i + flip];
649         float64 e1 = m[i + flip] ^ neg_real;
650         float64 e4 = e2;
651         float64 e3 = m[i + 1 - flip] ^ neg_imag;
652 
653         d[i] = float64_muladd(e2, e1, d[i], 0, fpst);
654         d[i + 1] = float64_muladd(e4, e3, d[i + 1], 0, fpst);
655     }
656     clear_tail(d, opr_sz, simd_maxsz(desc));
657 }
658 
659 /*
660  * Floating point comparisons producing an integer result (all 1s or all 0s).
661  * Note that EQ doesn't signal InvalidOp for QNaNs but GE and GT do.
662  * Softfloat routines return 0/1, which we convert to the 0/-1 Neon requires.
663  */
664 static uint16_t float16_ceq(float16 op1, float16 op2, float_status *stat)
665 {
666     return -float16_eq_quiet(op1, op2, stat);
667 }
668 
669 static uint32_t float32_ceq(float32 op1, float32 op2, float_status *stat)
670 {
671     return -float32_eq_quiet(op1, op2, stat);
672 }
673 
674 static uint16_t float16_cge(float16 op1, float16 op2, float_status *stat)
675 {
676     return -float16_le(op2, op1, stat);
677 }
678 
679 static uint32_t float32_cge(float32 op1, float32 op2, float_status *stat)
680 {
681     return -float32_le(op2, op1, stat);
682 }
683 
684 static uint16_t float16_cgt(float16 op1, float16 op2, float_status *stat)
685 {
686     return -float16_lt(op2, op1, stat);
687 }
688 
689 static uint32_t float32_cgt(float32 op1, float32 op2, float_status *stat)
690 {
691     return -float32_lt(op2, op1, stat);
692 }
693 
694 static uint16_t float16_acge(float16 op1, float16 op2, float_status *stat)
695 {
696     return -float16_le(float16_abs(op2), float16_abs(op1), stat);
697 }
698 
699 static uint32_t float32_acge(float32 op1, float32 op2, float_status *stat)
700 {
701     return -float32_le(float32_abs(op2), float32_abs(op1), stat);
702 }
703 
704 static uint16_t float16_acgt(float16 op1, float16 op2, float_status *stat)
705 {
706     return -float16_lt(float16_abs(op2), float16_abs(op1), stat);
707 }
708 
709 static uint32_t float32_acgt(float32 op1, float32 op2, float_status *stat)
710 {
711     return -float32_lt(float32_abs(op2), float32_abs(op1), stat);
712 }
713 
714 static int16_t vfp_tosszh(float16 x, void *fpstp)
715 {
716     float_status *fpst = fpstp;
717     if (float16_is_any_nan(x)) {
718         float_raise(float_flag_invalid, fpst);
719         return 0;
720     }
721     return float16_to_int16_round_to_zero(x, fpst);
722 }
723 
724 static uint16_t vfp_touszh(float16 x, void *fpstp)
725 {
726     float_status *fpst = fpstp;
727     if (float16_is_any_nan(x)) {
728         float_raise(float_flag_invalid, fpst);
729         return 0;
730     }
731     return float16_to_uint16_round_to_zero(x, fpst);
732 }
733 
734 #define DO_2OP(NAME, FUNC, TYPE) \
735 void HELPER(NAME)(void *vd, void *vn, void *stat, uint32_t desc)  \
736 {                                                                 \
737     intptr_t i, oprsz = simd_oprsz(desc);                         \
738     TYPE *d = vd, *n = vn;                                        \
739     for (i = 0; i < oprsz / sizeof(TYPE); i++) {                  \
740         d[i] = FUNC(n[i], stat);                                  \
741     }                                                             \
742     clear_tail(d, oprsz, simd_maxsz(desc));                       \
743 }
744 
745 DO_2OP(gvec_frecpe_h, helper_recpe_f16, float16)
746 DO_2OP(gvec_frecpe_s, helper_recpe_f32, float32)
747 DO_2OP(gvec_frecpe_d, helper_recpe_f64, float64)
748 
749 DO_2OP(gvec_frsqrte_h, helper_rsqrte_f16, float16)
750 DO_2OP(gvec_frsqrte_s, helper_rsqrte_f32, float32)
751 DO_2OP(gvec_frsqrte_d, helper_rsqrte_f64, float64)
752 
753 DO_2OP(gvec_sitos, helper_vfp_sitos, int32_t)
754 DO_2OP(gvec_uitos, helper_vfp_uitos, uint32_t)
755 DO_2OP(gvec_tosizs, helper_vfp_tosizs, float32)
756 DO_2OP(gvec_touizs, helper_vfp_touizs, float32)
757 DO_2OP(gvec_sstoh, int16_to_float16, int16_t)
758 DO_2OP(gvec_ustoh, uint16_to_float16, uint16_t)
759 DO_2OP(gvec_tosszh, vfp_tosszh, float16)
760 DO_2OP(gvec_touszh, vfp_touszh, float16)
761 
762 #define WRAP_CMP0_FWD(FN, CMPOP, TYPE)                          \
763     static TYPE TYPE##_##FN##0(TYPE op, float_status *stat)     \
764     {                                                           \
765         return TYPE##_##CMPOP(op, TYPE##_zero, stat);           \
766     }
767 
768 #define WRAP_CMP0_REV(FN, CMPOP, TYPE)                          \
769     static TYPE TYPE##_##FN##0(TYPE op, float_status *stat)    \
770     {                                                           \
771         return TYPE##_##CMPOP(TYPE##_zero, op, stat);           \
772     }
773 
774 #define DO_2OP_CMP0(FN, CMPOP, DIRN)                    \
775     WRAP_CMP0_##DIRN(FN, CMPOP, float16)                \
776     WRAP_CMP0_##DIRN(FN, CMPOP, float32)                \
777     DO_2OP(gvec_f##FN##0_h, float16_##FN##0, float16)   \
778     DO_2OP(gvec_f##FN##0_s, float32_##FN##0, float32)
779 
780 DO_2OP_CMP0(cgt, cgt, FWD)
781 DO_2OP_CMP0(cge, cge, FWD)
782 DO_2OP_CMP0(ceq, ceq, FWD)
783 DO_2OP_CMP0(clt, cgt, REV)
784 DO_2OP_CMP0(cle, cge, REV)
785 
786 #undef DO_2OP
787 #undef DO_2OP_CMP0
788 
789 /* Floating-point trigonometric starting value.
790  * See the ARM ARM pseudocode function FPTrigSMul.
791  */
792 static float16 float16_ftsmul(float16 op1, uint16_t op2, float_status *stat)
793 {
794     float16 result = float16_mul(op1, op1, stat);
795     if (!float16_is_any_nan(result)) {
796         result = float16_set_sign(result, op2 & 1);
797     }
798     return result;
799 }
800 
801 static float32 float32_ftsmul(float32 op1, uint32_t op2, float_status *stat)
802 {
803     float32 result = float32_mul(op1, op1, stat);
804     if (!float32_is_any_nan(result)) {
805         result = float32_set_sign(result, op2 & 1);
806     }
807     return result;
808 }
809 
810 static float64 float64_ftsmul(float64 op1, uint64_t op2, float_status *stat)
811 {
812     float64 result = float64_mul(op1, op1, stat);
813     if (!float64_is_any_nan(result)) {
814         result = float64_set_sign(result, op2 & 1);
815     }
816     return result;
817 }
818 
819 static float16 float16_abd(float16 op1, float16 op2, float_status *stat)
820 {
821     return float16_abs(float16_sub(op1, op2, stat));
822 }
823 
824 static float32 float32_abd(float32 op1, float32 op2, float_status *stat)
825 {
826     return float32_abs(float32_sub(op1, op2, stat));
827 }
828 
829 /*
830  * Reciprocal step. These are the AArch32 version which uses a
831  * non-fused multiply-and-subtract.
832  */
833 static float16 float16_recps_nf(float16 op1, float16 op2, float_status *stat)
834 {
835     op1 = float16_squash_input_denormal(op1, stat);
836     op2 = float16_squash_input_denormal(op2, stat);
837 
838     if ((float16_is_infinity(op1) && float16_is_zero(op2)) ||
839         (float16_is_infinity(op2) && float16_is_zero(op1))) {
840         return float16_two;
841     }
842     return float16_sub(float16_two, float16_mul(op1, op2, stat), stat);
843 }
844 
845 static float32 float32_recps_nf(float32 op1, float32 op2, float_status *stat)
846 {
847     op1 = float32_squash_input_denormal(op1, stat);
848     op2 = float32_squash_input_denormal(op2, stat);
849 
850     if ((float32_is_infinity(op1) && float32_is_zero(op2)) ||
851         (float32_is_infinity(op2) && float32_is_zero(op1))) {
852         return float32_two;
853     }
854     return float32_sub(float32_two, float32_mul(op1, op2, stat), stat);
855 }
856 
857 /* Reciprocal square-root step. AArch32 non-fused semantics. */
858 static float16 float16_rsqrts_nf(float16 op1, float16 op2, float_status *stat)
859 {
860     op1 = float16_squash_input_denormal(op1, stat);
861     op2 = float16_squash_input_denormal(op2, stat);
862 
863     if ((float16_is_infinity(op1) && float16_is_zero(op2)) ||
864         (float16_is_infinity(op2) && float16_is_zero(op1))) {
865         return float16_one_point_five;
866     }
867     op1 = float16_sub(float16_three, float16_mul(op1, op2, stat), stat);
868     return float16_div(op1, float16_two, stat);
869 }
870 
871 static float32 float32_rsqrts_nf(float32 op1, float32 op2, float_status *stat)
872 {
873     op1 = float32_squash_input_denormal(op1, stat);
874     op2 = float32_squash_input_denormal(op2, stat);
875 
876     if ((float32_is_infinity(op1) && float32_is_zero(op2)) ||
877         (float32_is_infinity(op2) && float32_is_zero(op1))) {
878         return float32_one_point_five;
879     }
880     op1 = float32_sub(float32_three, float32_mul(op1, op2, stat), stat);
881     return float32_div(op1, float32_two, stat);
882 }
883 
884 #define DO_3OP(NAME, FUNC, TYPE) \
885 void HELPER(NAME)(void *vd, void *vn, void *vm, void *stat, uint32_t desc) \
886 {                                                                          \
887     intptr_t i, oprsz = simd_oprsz(desc);                                  \
888     TYPE *d = vd, *n = vn, *m = vm;                                        \
889     for (i = 0; i < oprsz / sizeof(TYPE); i++) {                           \
890         d[i] = FUNC(n[i], m[i], stat);                                     \
891     }                                                                      \
892     clear_tail(d, oprsz, simd_maxsz(desc));                                \
893 }
894 
895 DO_3OP(gvec_fadd_h, float16_add, float16)
896 DO_3OP(gvec_fadd_s, float32_add, float32)
897 DO_3OP(gvec_fadd_d, float64_add, float64)
898 
899 DO_3OP(gvec_fsub_h, float16_sub, float16)
900 DO_3OP(gvec_fsub_s, float32_sub, float32)
901 DO_3OP(gvec_fsub_d, float64_sub, float64)
902 
903 DO_3OP(gvec_fmul_h, float16_mul, float16)
904 DO_3OP(gvec_fmul_s, float32_mul, float32)
905 DO_3OP(gvec_fmul_d, float64_mul, float64)
906 
907 DO_3OP(gvec_ftsmul_h, float16_ftsmul, float16)
908 DO_3OP(gvec_ftsmul_s, float32_ftsmul, float32)
909 DO_3OP(gvec_ftsmul_d, float64_ftsmul, float64)
910 
911 DO_3OP(gvec_fabd_h, float16_abd, float16)
912 DO_3OP(gvec_fabd_s, float32_abd, float32)
913 
914 DO_3OP(gvec_fceq_h, float16_ceq, float16)
915 DO_3OP(gvec_fceq_s, float32_ceq, float32)
916 
917 DO_3OP(gvec_fcge_h, float16_cge, float16)
918 DO_3OP(gvec_fcge_s, float32_cge, float32)
919 
920 DO_3OP(gvec_fcgt_h, float16_cgt, float16)
921 DO_3OP(gvec_fcgt_s, float32_cgt, float32)
922 
923 DO_3OP(gvec_facge_h, float16_acge, float16)
924 DO_3OP(gvec_facge_s, float32_acge, float32)
925 
926 DO_3OP(gvec_facgt_h, float16_acgt, float16)
927 DO_3OP(gvec_facgt_s, float32_acgt, float32)
928 
929 DO_3OP(gvec_fmax_h, float16_max, float16)
930 DO_3OP(gvec_fmax_s, float32_max, float32)
931 
932 DO_3OP(gvec_fmin_h, float16_min, float16)
933 DO_3OP(gvec_fmin_s, float32_min, float32)
934 
935 DO_3OP(gvec_fmaxnum_h, float16_maxnum, float16)
936 DO_3OP(gvec_fmaxnum_s, float32_maxnum, float32)
937 
938 DO_3OP(gvec_fminnum_h, float16_minnum, float16)
939 DO_3OP(gvec_fminnum_s, float32_minnum, float32)
940 
941 DO_3OP(gvec_recps_nf_h, float16_recps_nf, float16)
942 DO_3OP(gvec_recps_nf_s, float32_recps_nf, float32)
943 
944 DO_3OP(gvec_rsqrts_nf_h, float16_rsqrts_nf, float16)
945 DO_3OP(gvec_rsqrts_nf_s, float32_rsqrts_nf, float32)
946 
947 #ifdef TARGET_AARCH64
948 
949 DO_3OP(gvec_recps_h, helper_recpsf_f16, float16)
950 DO_3OP(gvec_recps_s, helper_recpsf_f32, float32)
951 DO_3OP(gvec_recps_d, helper_recpsf_f64, float64)
952 
953 DO_3OP(gvec_rsqrts_h, helper_rsqrtsf_f16, float16)
954 DO_3OP(gvec_rsqrts_s, helper_rsqrtsf_f32, float32)
955 DO_3OP(gvec_rsqrts_d, helper_rsqrtsf_f64, float64)
956 
957 #endif
958 #undef DO_3OP
959 
960 /* Non-fused multiply-add (unlike float16_muladd etc, which are fused) */
961 static float16 float16_muladd_nf(float16 dest, float16 op1, float16 op2,
962                                  float_status *stat)
963 {
964     return float16_add(dest, float16_mul(op1, op2, stat), stat);
965 }
966 
967 static float32 float32_muladd_nf(float32 dest, float32 op1, float32 op2,
968                                  float_status *stat)
969 {
970     return float32_add(dest, float32_mul(op1, op2, stat), stat);
971 }
972 
973 static float16 float16_mulsub_nf(float16 dest, float16 op1, float16 op2,
974                                  float_status *stat)
975 {
976     return float16_sub(dest, float16_mul(op1, op2, stat), stat);
977 }
978 
979 static float32 float32_mulsub_nf(float32 dest, float32 op1, float32 op2,
980                                  float_status *stat)
981 {
982     return float32_sub(dest, float32_mul(op1, op2, stat), stat);
983 }
984 
985 /* Fused versions; these have the semantics Neon VFMA/VFMS want */
986 static float16 float16_muladd_f(float16 dest, float16 op1, float16 op2,
987                                 float_status *stat)
988 {
989     return float16_muladd(op1, op2, dest, 0, stat);
990 }
991 
992 static float32 float32_muladd_f(float32 dest, float32 op1, float32 op2,
993                                  float_status *stat)
994 {
995     return float32_muladd(op1, op2, dest, 0, stat);
996 }
997 
998 static float16 float16_mulsub_f(float16 dest, float16 op1, float16 op2,
999                                  float_status *stat)
1000 {
1001     return float16_muladd(float16_chs(op1), op2, dest, 0, stat);
1002 }
1003 
1004 static float32 float32_mulsub_f(float32 dest, float32 op1, float32 op2,
1005                                  float_status *stat)
1006 {
1007     return float32_muladd(float32_chs(op1), op2, dest, 0, stat);
1008 }
1009 
1010 #define DO_MULADD(NAME, FUNC, TYPE)                                     \
1011 void HELPER(NAME)(void *vd, void *vn, void *vm, void *stat, uint32_t desc) \
1012 {                                                                          \
1013     intptr_t i, oprsz = simd_oprsz(desc);                                  \
1014     TYPE *d = vd, *n = vn, *m = vm;                                        \
1015     for (i = 0; i < oprsz / sizeof(TYPE); i++) {                           \
1016         d[i] = FUNC(d[i], n[i], m[i], stat);                               \
1017     }                                                                      \
1018     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1019 }
1020 
1021 DO_MULADD(gvec_fmla_h, float16_muladd_nf, float16)
1022 DO_MULADD(gvec_fmla_s, float32_muladd_nf, float32)
1023 
1024 DO_MULADD(gvec_fmls_h, float16_mulsub_nf, float16)
1025 DO_MULADD(gvec_fmls_s, float32_mulsub_nf, float32)
1026 
1027 DO_MULADD(gvec_vfma_h, float16_muladd_f, float16)
1028 DO_MULADD(gvec_vfma_s, float32_muladd_f, float32)
1029 
1030 DO_MULADD(gvec_vfms_h, float16_mulsub_f, float16)
1031 DO_MULADD(gvec_vfms_s, float32_mulsub_f, float32)
1032 
1033 /* For the indexed ops, SVE applies the index per 128-bit vector segment.
1034  * For AdvSIMD, there is of course only one such vector segment.
1035  */
1036 
1037 #define DO_MUL_IDX(NAME, TYPE, H) \
1038 void HELPER(NAME)(void *vd, void *vn, void *vm, uint32_t desc) \
1039 {                                                                          \
1040     intptr_t i, j, oprsz = simd_oprsz(desc), segment = 16 / sizeof(TYPE);  \
1041     intptr_t idx = simd_data(desc);                                        \
1042     TYPE *d = vd, *n = vn, *m = vm;                                        \
1043     for (i = 0; i < oprsz / sizeof(TYPE); i += segment) {                  \
1044         TYPE mm = m[H(i + idx)];                                           \
1045         for (j = 0; j < segment; j++) {                                    \
1046             d[i + j] = n[i + j] * mm;                                      \
1047         }                                                                  \
1048     }                                                                      \
1049     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1050 }
1051 
1052 DO_MUL_IDX(gvec_mul_idx_h, uint16_t, H2)
1053 DO_MUL_IDX(gvec_mul_idx_s, uint32_t, H4)
1054 DO_MUL_IDX(gvec_mul_idx_d, uint64_t, )
1055 
1056 #undef DO_MUL_IDX
1057 
1058 #define DO_MLA_IDX(NAME, TYPE, OP, H) \
1059 void HELPER(NAME)(void *vd, void *vn, void *vm, void *va, uint32_t desc)   \
1060 {                                                                          \
1061     intptr_t i, j, oprsz = simd_oprsz(desc), segment = 16 / sizeof(TYPE);  \
1062     intptr_t idx = simd_data(desc);                                        \
1063     TYPE *d = vd, *n = vn, *m = vm, *a = va;                               \
1064     for (i = 0; i < oprsz / sizeof(TYPE); i += segment) {                  \
1065         TYPE mm = m[H(i + idx)];                                           \
1066         for (j = 0; j < segment; j++) {                                    \
1067             d[i + j] = a[i + j] OP n[i + j] * mm;                          \
1068         }                                                                  \
1069     }                                                                      \
1070     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1071 }
1072 
1073 DO_MLA_IDX(gvec_mla_idx_h, uint16_t, +, H2)
1074 DO_MLA_IDX(gvec_mla_idx_s, uint32_t, +, H4)
1075 DO_MLA_IDX(gvec_mla_idx_d, uint64_t, +,   )
1076 
1077 DO_MLA_IDX(gvec_mls_idx_h, uint16_t, -, H2)
1078 DO_MLA_IDX(gvec_mls_idx_s, uint32_t, -, H4)
1079 DO_MLA_IDX(gvec_mls_idx_d, uint64_t, -,   )
1080 
1081 #undef DO_MLA_IDX
1082 
1083 #define DO_FMUL_IDX(NAME, TYPE, H) \
1084 void HELPER(NAME)(void *vd, void *vn, void *vm, void *stat, uint32_t desc) \
1085 {                                                                          \
1086     intptr_t i, j, oprsz = simd_oprsz(desc), segment = 16 / sizeof(TYPE);  \
1087     intptr_t idx = simd_data(desc);                                        \
1088     TYPE *d = vd, *n = vn, *m = vm;                                        \
1089     for (i = 0; i < oprsz / sizeof(TYPE); i += segment) {                  \
1090         TYPE mm = m[H(i + idx)];                                           \
1091         for (j = 0; j < segment; j++) {                                    \
1092             d[i + j] = TYPE##_mul(n[i + j], mm, stat);                     \
1093         }                                                                  \
1094     }                                                                      \
1095     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1096 }
1097 
1098 DO_FMUL_IDX(gvec_fmul_idx_h, float16, H2)
1099 DO_FMUL_IDX(gvec_fmul_idx_s, float32, H4)
1100 DO_FMUL_IDX(gvec_fmul_idx_d, float64, )
1101 
1102 #undef DO_FMUL_IDX
1103 
1104 #define DO_FMLA_IDX(NAME, TYPE, H)                                         \
1105 void HELPER(NAME)(void *vd, void *vn, void *vm, void *va,                  \
1106                   void *stat, uint32_t desc)                               \
1107 {                                                                          \
1108     intptr_t i, j, oprsz = simd_oprsz(desc), segment = 16 / sizeof(TYPE);  \
1109     TYPE op1_neg = extract32(desc, SIMD_DATA_SHIFT, 1);                    \
1110     intptr_t idx = desc >> (SIMD_DATA_SHIFT + 1);                          \
1111     TYPE *d = vd, *n = vn, *m = vm, *a = va;                               \
1112     op1_neg <<= (8 * sizeof(TYPE) - 1);                                    \
1113     for (i = 0; i < oprsz / sizeof(TYPE); i += segment) {                  \
1114         TYPE mm = m[H(i + idx)];                                           \
1115         for (j = 0; j < segment; j++) {                                    \
1116             d[i + j] = TYPE##_muladd(n[i + j] ^ op1_neg,                   \
1117                                      mm, a[i + j], 0, stat);               \
1118         }                                                                  \
1119     }                                                                      \
1120     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1121 }
1122 
1123 DO_FMLA_IDX(gvec_fmla_idx_h, float16, H2)
1124 DO_FMLA_IDX(gvec_fmla_idx_s, float32, H4)
1125 DO_FMLA_IDX(gvec_fmla_idx_d, float64, )
1126 
1127 #undef DO_FMLA_IDX
1128 
1129 #define DO_SAT(NAME, WTYPE, TYPEN, TYPEM, OP, MIN, MAX) \
1130 void HELPER(NAME)(void *vd, void *vq, void *vn, void *vm, uint32_t desc)   \
1131 {                                                                          \
1132     intptr_t i, oprsz = simd_oprsz(desc);                                  \
1133     TYPEN *d = vd, *n = vn; TYPEM *m = vm;                                 \
1134     bool q = false;                                                        \
1135     for (i = 0; i < oprsz / sizeof(TYPEN); i++) {                          \
1136         WTYPE dd = (WTYPE)n[i] OP m[i];                                    \
1137         if (dd < MIN) {                                                    \
1138             dd = MIN;                                                      \
1139             q = true;                                                      \
1140         } else if (dd > MAX) {                                             \
1141             dd = MAX;                                                      \
1142             q = true;                                                      \
1143         }                                                                  \
1144         d[i] = dd;                                                         \
1145     }                                                                      \
1146     if (q) {                                                               \
1147         uint32_t *qc = vq;                                                 \
1148         qc[0] = 1;                                                         \
1149     }                                                                      \
1150     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1151 }
1152 
1153 DO_SAT(gvec_uqadd_b, int, uint8_t, uint8_t, +, 0, UINT8_MAX)
1154 DO_SAT(gvec_uqadd_h, int, uint16_t, uint16_t, +, 0, UINT16_MAX)
1155 DO_SAT(gvec_uqadd_s, int64_t, uint32_t, uint32_t, +, 0, UINT32_MAX)
1156 
1157 DO_SAT(gvec_sqadd_b, int, int8_t, int8_t, +, INT8_MIN, INT8_MAX)
1158 DO_SAT(gvec_sqadd_h, int, int16_t, int16_t, +, INT16_MIN, INT16_MAX)
1159 DO_SAT(gvec_sqadd_s, int64_t, int32_t, int32_t, +, INT32_MIN, INT32_MAX)
1160 
1161 DO_SAT(gvec_uqsub_b, int, uint8_t, uint8_t, -, 0, UINT8_MAX)
1162 DO_SAT(gvec_uqsub_h, int, uint16_t, uint16_t, -, 0, UINT16_MAX)
1163 DO_SAT(gvec_uqsub_s, int64_t, uint32_t, uint32_t, -, 0, UINT32_MAX)
1164 
1165 DO_SAT(gvec_sqsub_b, int, int8_t, int8_t, -, INT8_MIN, INT8_MAX)
1166 DO_SAT(gvec_sqsub_h, int, int16_t, int16_t, -, INT16_MIN, INT16_MAX)
1167 DO_SAT(gvec_sqsub_s, int64_t, int32_t, int32_t, -, INT32_MIN, INT32_MAX)
1168 
1169 #undef DO_SAT
1170 
1171 void HELPER(gvec_uqadd_d)(void *vd, void *vq, void *vn,
1172                           void *vm, uint32_t desc)
1173 {
1174     intptr_t i, oprsz = simd_oprsz(desc);
1175     uint64_t *d = vd, *n = vn, *m = vm;
1176     bool q = false;
1177 
1178     for (i = 0; i < oprsz / 8; i++) {
1179         uint64_t nn = n[i], mm = m[i], dd = nn + mm;
1180         if (dd < nn) {
1181             dd = UINT64_MAX;
1182             q = true;
1183         }
1184         d[i] = dd;
1185     }
1186     if (q) {
1187         uint32_t *qc = vq;
1188         qc[0] = 1;
1189     }
1190     clear_tail(d, oprsz, simd_maxsz(desc));
1191 }
1192 
1193 void HELPER(gvec_uqsub_d)(void *vd, void *vq, void *vn,
1194                           void *vm, uint32_t desc)
1195 {
1196     intptr_t i, oprsz = simd_oprsz(desc);
1197     uint64_t *d = vd, *n = vn, *m = vm;
1198     bool q = false;
1199 
1200     for (i = 0; i < oprsz / 8; i++) {
1201         uint64_t nn = n[i], mm = m[i], dd = nn - mm;
1202         if (nn < mm) {
1203             dd = 0;
1204             q = true;
1205         }
1206         d[i] = dd;
1207     }
1208     if (q) {
1209         uint32_t *qc = vq;
1210         qc[0] = 1;
1211     }
1212     clear_tail(d, oprsz, simd_maxsz(desc));
1213 }
1214 
1215 void HELPER(gvec_sqadd_d)(void *vd, void *vq, void *vn,
1216                           void *vm, uint32_t desc)
1217 {
1218     intptr_t i, oprsz = simd_oprsz(desc);
1219     int64_t *d = vd, *n = vn, *m = vm;
1220     bool q = false;
1221 
1222     for (i = 0; i < oprsz / 8; i++) {
1223         int64_t nn = n[i], mm = m[i], dd = nn + mm;
1224         if (((dd ^ nn) & ~(nn ^ mm)) & INT64_MIN) {
1225             dd = (nn >> 63) ^ ~INT64_MIN;
1226             q = true;
1227         }
1228         d[i] = dd;
1229     }
1230     if (q) {
1231         uint32_t *qc = vq;
1232         qc[0] = 1;
1233     }
1234     clear_tail(d, oprsz, simd_maxsz(desc));
1235 }
1236 
1237 void HELPER(gvec_sqsub_d)(void *vd, void *vq, void *vn,
1238                           void *vm, uint32_t desc)
1239 {
1240     intptr_t i, oprsz = simd_oprsz(desc);
1241     int64_t *d = vd, *n = vn, *m = vm;
1242     bool q = false;
1243 
1244     for (i = 0; i < oprsz / 8; i++) {
1245         int64_t nn = n[i], mm = m[i], dd = nn - mm;
1246         if (((dd ^ nn) & (nn ^ mm)) & INT64_MIN) {
1247             dd = (nn >> 63) ^ ~INT64_MIN;
1248             q = true;
1249         }
1250         d[i] = dd;
1251     }
1252     if (q) {
1253         uint32_t *qc = vq;
1254         qc[0] = 1;
1255     }
1256     clear_tail(d, oprsz, simd_maxsz(desc));
1257 }
1258 
1259 
1260 #define DO_SRA(NAME, TYPE)                              \
1261 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1262 {                                                       \
1263     intptr_t i, oprsz = simd_oprsz(desc);               \
1264     int shift = simd_data(desc);                        \
1265     TYPE *d = vd, *n = vn;                              \
1266     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1267         d[i] += n[i] >> shift;                          \
1268     }                                                   \
1269     clear_tail(d, oprsz, simd_maxsz(desc));             \
1270 }
1271 
1272 DO_SRA(gvec_ssra_b, int8_t)
1273 DO_SRA(gvec_ssra_h, int16_t)
1274 DO_SRA(gvec_ssra_s, int32_t)
1275 DO_SRA(gvec_ssra_d, int64_t)
1276 
1277 DO_SRA(gvec_usra_b, uint8_t)
1278 DO_SRA(gvec_usra_h, uint16_t)
1279 DO_SRA(gvec_usra_s, uint32_t)
1280 DO_SRA(gvec_usra_d, uint64_t)
1281 
1282 #undef DO_SRA
1283 
1284 #define DO_RSHR(NAME, TYPE)                             \
1285 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1286 {                                                       \
1287     intptr_t i, oprsz = simd_oprsz(desc);               \
1288     int shift = simd_data(desc);                        \
1289     TYPE *d = vd, *n = vn;                              \
1290     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1291         TYPE tmp = n[i] >> (shift - 1);                 \
1292         d[i] = (tmp >> 1) + (tmp & 1);                  \
1293     }                                                   \
1294     clear_tail(d, oprsz, simd_maxsz(desc));             \
1295 }
1296 
1297 DO_RSHR(gvec_srshr_b, int8_t)
1298 DO_RSHR(gvec_srshr_h, int16_t)
1299 DO_RSHR(gvec_srshr_s, int32_t)
1300 DO_RSHR(gvec_srshr_d, int64_t)
1301 
1302 DO_RSHR(gvec_urshr_b, uint8_t)
1303 DO_RSHR(gvec_urshr_h, uint16_t)
1304 DO_RSHR(gvec_urshr_s, uint32_t)
1305 DO_RSHR(gvec_urshr_d, uint64_t)
1306 
1307 #undef DO_RSHR
1308 
1309 #define DO_RSRA(NAME, TYPE)                             \
1310 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1311 {                                                       \
1312     intptr_t i, oprsz = simd_oprsz(desc);               \
1313     int shift = simd_data(desc);                        \
1314     TYPE *d = vd, *n = vn;                              \
1315     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1316         TYPE tmp = n[i] >> (shift - 1);                 \
1317         d[i] += (tmp >> 1) + (tmp & 1);                 \
1318     }                                                   \
1319     clear_tail(d, oprsz, simd_maxsz(desc));             \
1320 }
1321 
1322 DO_RSRA(gvec_srsra_b, int8_t)
1323 DO_RSRA(gvec_srsra_h, int16_t)
1324 DO_RSRA(gvec_srsra_s, int32_t)
1325 DO_RSRA(gvec_srsra_d, int64_t)
1326 
1327 DO_RSRA(gvec_ursra_b, uint8_t)
1328 DO_RSRA(gvec_ursra_h, uint16_t)
1329 DO_RSRA(gvec_ursra_s, uint32_t)
1330 DO_RSRA(gvec_ursra_d, uint64_t)
1331 
1332 #undef DO_RSRA
1333 
1334 #define DO_SRI(NAME, TYPE)                              \
1335 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1336 {                                                       \
1337     intptr_t i, oprsz = simd_oprsz(desc);               \
1338     int shift = simd_data(desc);                        \
1339     TYPE *d = vd, *n = vn;                              \
1340     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1341         d[i] = deposit64(d[i], 0, sizeof(TYPE) * 8 - shift, n[i] >> shift); \
1342     }                                                   \
1343     clear_tail(d, oprsz, simd_maxsz(desc));             \
1344 }
1345 
1346 DO_SRI(gvec_sri_b, uint8_t)
1347 DO_SRI(gvec_sri_h, uint16_t)
1348 DO_SRI(gvec_sri_s, uint32_t)
1349 DO_SRI(gvec_sri_d, uint64_t)
1350 
1351 #undef DO_SRI
1352 
1353 #define DO_SLI(NAME, TYPE)                              \
1354 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1355 {                                                       \
1356     intptr_t i, oprsz = simd_oprsz(desc);               \
1357     int shift = simd_data(desc);                        \
1358     TYPE *d = vd, *n = vn;                              \
1359     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1360         d[i] = deposit64(d[i], shift, sizeof(TYPE) * 8 - shift, n[i]); \
1361     }                                                   \
1362     clear_tail(d, oprsz, simd_maxsz(desc));             \
1363 }
1364 
1365 DO_SLI(gvec_sli_b, uint8_t)
1366 DO_SLI(gvec_sli_h, uint16_t)
1367 DO_SLI(gvec_sli_s, uint32_t)
1368 DO_SLI(gvec_sli_d, uint64_t)
1369 
1370 #undef DO_SLI
1371 
1372 /*
1373  * Convert float16 to float32, raising no exceptions and
1374  * preserving exceptional values, including SNaN.
1375  * This is effectively an unpack+repack operation.
1376  */
1377 static float32 float16_to_float32_by_bits(uint32_t f16, bool fz16)
1378 {
1379     const int f16_bias = 15;
1380     const int f32_bias = 127;
1381     uint32_t sign = extract32(f16, 15, 1);
1382     uint32_t exp = extract32(f16, 10, 5);
1383     uint32_t frac = extract32(f16, 0, 10);
1384 
1385     if (exp == 0x1f) {
1386         /* Inf or NaN */
1387         exp = 0xff;
1388     } else if (exp == 0) {
1389         /* Zero or denormal.  */
1390         if (frac != 0) {
1391             if (fz16) {
1392                 frac = 0;
1393             } else {
1394                 /*
1395                  * Denormal; these are all normal float32.
1396                  * Shift the fraction so that the msb is at bit 11,
1397                  * then remove bit 11 as the implicit bit of the
1398                  * normalized float32.  Note that we still go through
1399                  * the shift for normal numbers below, to put the
1400                  * float32 fraction at the right place.
1401                  */
1402                 int shift = clz32(frac) - 21;
1403                 frac = (frac << shift) & 0x3ff;
1404                 exp = f32_bias - f16_bias - shift + 1;
1405             }
1406         }
1407     } else {
1408         /* Normal number; adjust the bias.  */
1409         exp += f32_bias - f16_bias;
1410     }
1411     sign <<= 31;
1412     exp <<= 23;
1413     frac <<= 23 - 10;
1414 
1415     return sign | exp | frac;
1416 }
1417 
1418 static uint64_t load4_f16(uint64_t *ptr, int is_q, int is_2)
1419 {
1420     /*
1421      * Branchless load of u32[0], u64[0], u32[1], or u64[1].
1422      * Load the 2nd qword iff is_q & is_2.
1423      * Shift to the 2nd dword iff !is_q & is_2.
1424      * For !is_q & !is_2, the upper bits of the result are garbage.
1425      */
1426     return ptr[is_q & is_2] >> ((is_2 & ~is_q) << 5);
1427 }
1428 
1429 /*
1430  * Note that FMLAL requires oprsz == 8 or oprsz == 16,
1431  * as there is not yet SVE versions that might use blocking.
1432  */
1433 
1434 static void do_fmlal(float32 *d, void *vn, void *vm, float_status *fpst,
1435                      uint32_t desc, bool fz16)
1436 {
1437     intptr_t i, oprsz = simd_oprsz(desc);
1438     int is_s = extract32(desc, SIMD_DATA_SHIFT, 1);
1439     int is_2 = extract32(desc, SIMD_DATA_SHIFT + 1, 1);
1440     int is_q = oprsz == 16;
1441     uint64_t n_4, m_4;
1442 
1443     /* Pre-load all of the f16 data, avoiding overlap issues.  */
1444     n_4 = load4_f16(vn, is_q, is_2);
1445     m_4 = load4_f16(vm, is_q, is_2);
1446 
1447     /* Negate all inputs for FMLSL at once.  */
1448     if (is_s) {
1449         n_4 ^= 0x8000800080008000ull;
1450     }
1451 
1452     for (i = 0; i < oprsz / 4; i++) {
1453         float32 n_1 = float16_to_float32_by_bits(n_4 >> (i * 16), fz16);
1454         float32 m_1 = float16_to_float32_by_bits(m_4 >> (i * 16), fz16);
1455         d[H4(i)] = float32_muladd(n_1, m_1, d[H4(i)], 0, fpst);
1456     }
1457     clear_tail(d, oprsz, simd_maxsz(desc));
1458 }
1459 
1460 void HELPER(gvec_fmlal_a32)(void *vd, void *vn, void *vm,
1461                             void *venv, uint32_t desc)
1462 {
1463     CPUARMState *env = venv;
1464     do_fmlal(vd, vn, vm, &env->vfp.standard_fp_status, desc,
1465              get_flush_inputs_to_zero(&env->vfp.fp_status_f16));
1466 }
1467 
1468 void HELPER(gvec_fmlal_a64)(void *vd, void *vn, void *vm,
1469                             void *venv, uint32_t desc)
1470 {
1471     CPUARMState *env = venv;
1472     do_fmlal(vd, vn, vm, &env->vfp.fp_status, desc,
1473              get_flush_inputs_to_zero(&env->vfp.fp_status_f16));
1474 }
1475 
1476 static void do_fmlal_idx(float32 *d, void *vn, void *vm, float_status *fpst,
1477                          uint32_t desc, bool fz16)
1478 {
1479     intptr_t i, oprsz = simd_oprsz(desc);
1480     int is_s = extract32(desc, SIMD_DATA_SHIFT, 1);
1481     int is_2 = extract32(desc, SIMD_DATA_SHIFT + 1, 1);
1482     int index = extract32(desc, SIMD_DATA_SHIFT + 2, 3);
1483     int is_q = oprsz == 16;
1484     uint64_t n_4;
1485     float32 m_1;
1486 
1487     /* Pre-load all of the f16 data, avoiding overlap issues.  */
1488     n_4 = load4_f16(vn, is_q, is_2);
1489 
1490     /* Negate all inputs for FMLSL at once.  */
1491     if (is_s) {
1492         n_4 ^= 0x8000800080008000ull;
1493     }
1494 
1495     m_1 = float16_to_float32_by_bits(((float16 *)vm)[H2(index)], fz16);
1496 
1497     for (i = 0; i < oprsz / 4; i++) {
1498         float32 n_1 = float16_to_float32_by_bits(n_4 >> (i * 16), fz16);
1499         d[H4(i)] = float32_muladd(n_1, m_1, d[H4(i)], 0, fpst);
1500     }
1501     clear_tail(d, oprsz, simd_maxsz(desc));
1502 }
1503 
1504 void HELPER(gvec_fmlal_idx_a32)(void *vd, void *vn, void *vm,
1505                                 void *venv, uint32_t desc)
1506 {
1507     CPUARMState *env = venv;
1508     do_fmlal_idx(vd, vn, vm, &env->vfp.standard_fp_status, desc,
1509                  get_flush_inputs_to_zero(&env->vfp.fp_status_f16));
1510 }
1511 
1512 void HELPER(gvec_fmlal_idx_a64)(void *vd, void *vn, void *vm,
1513                                 void *venv, uint32_t desc)
1514 {
1515     CPUARMState *env = venv;
1516     do_fmlal_idx(vd, vn, vm, &env->vfp.fp_status, desc,
1517                  get_flush_inputs_to_zero(&env->vfp.fp_status_f16));
1518 }
1519 
1520 void HELPER(gvec_sshl_b)(void *vd, void *vn, void *vm, uint32_t desc)
1521 {
1522     intptr_t i, opr_sz = simd_oprsz(desc);
1523     int8_t *d = vd, *n = vn, *m = vm;
1524 
1525     for (i = 0; i < opr_sz; ++i) {
1526         int8_t mm = m[i];
1527         int8_t nn = n[i];
1528         int8_t res = 0;
1529         if (mm >= 0) {
1530             if (mm < 8) {
1531                 res = nn << mm;
1532             }
1533         } else {
1534             res = nn >> (mm > -8 ? -mm : 7);
1535         }
1536         d[i] = res;
1537     }
1538     clear_tail(d, opr_sz, simd_maxsz(desc));
1539 }
1540 
1541 void HELPER(gvec_sshl_h)(void *vd, void *vn, void *vm, uint32_t desc)
1542 {
1543     intptr_t i, opr_sz = simd_oprsz(desc);
1544     int16_t *d = vd, *n = vn, *m = vm;
1545 
1546     for (i = 0; i < opr_sz / 2; ++i) {
1547         int8_t mm = m[i];   /* only 8 bits of shift are significant */
1548         int16_t nn = n[i];
1549         int16_t res = 0;
1550         if (mm >= 0) {
1551             if (mm < 16) {
1552                 res = nn << mm;
1553             }
1554         } else {
1555             res = nn >> (mm > -16 ? -mm : 15);
1556         }
1557         d[i] = res;
1558     }
1559     clear_tail(d, opr_sz, simd_maxsz(desc));
1560 }
1561 
1562 void HELPER(gvec_ushl_b)(void *vd, void *vn, void *vm, uint32_t desc)
1563 {
1564     intptr_t i, opr_sz = simd_oprsz(desc);
1565     uint8_t *d = vd, *n = vn, *m = vm;
1566 
1567     for (i = 0; i < opr_sz; ++i) {
1568         int8_t mm = m[i];
1569         uint8_t nn = n[i];
1570         uint8_t res = 0;
1571         if (mm >= 0) {
1572             if (mm < 8) {
1573                 res = nn << mm;
1574             }
1575         } else {
1576             if (mm > -8) {
1577                 res = nn >> -mm;
1578             }
1579         }
1580         d[i] = res;
1581     }
1582     clear_tail(d, opr_sz, simd_maxsz(desc));
1583 }
1584 
1585 void HELPER(gvec_ushl_h)(void *vd, void *vn, void *vm, uint32_t desc)
1586 {
1587     intptr_t i, opr_sz = simd_oprsz(desc);
1588     uint16_t *d = vd, *n = vn, *m = vm;
1589 
1590     for (i = 0; i < opr_sz / 2; ++i) {
1591         int8_t mm = m[i];   /* only 8 bits of shift are significant */
1592         uint16_t nn = n[i];
1593         uint16_t res = 0;
1594         if (mm >= 0) {
1595             if (mm < 16) {
1596                 res = nn << mm;
1597             }
1598         } else {
1599             if (mm > -16) {
1600                 res = nn >> -mm;
1601             }
1602         }
1603         d[i] = res;
1604     }
1605     clear_tail(d, opr_sz, simd_maxsz(desc));
1606 }
1607 
1608 /*
1609  * 8x8->8 polynomial multiply.
1610  *
1611  * Polynomial multiplication is like integer multiplication except the
1612  * partial products are XORed, not added.
1613  *
1614  * TODO: expose this as a generic vector operation, as it is a common
1615  * crypto building block.
1616  */
1617 void HELPER(gvec_pmul_b)(void *vd, void *vn, void *vm, uint32_t desc)
1618 {
1619     intptr_t i, j, opr_sz = simd_oprsz(desc);
1620     uint64_t *d = vd, *n = vn, *m = vm;
1621 
1622     for (i = 0; i < opr_sz / 8; ++i) {
1623         uint64_t nn = n[i];
1624         uint64_t mm = m[i];
1625         uint64_t rr = 0;
1626 
1627         for (j = 0; j < 8; ++j) {
1628             uint64_t mask = (nn & 0x0101010101010101ull) * 0xff;
1629             rr ^= mm & mask;
1630             mm = (mm << 1) & 0xfefefefefefefefeull;
1631             nn >>= 1;
1632         }
1633         d[i] = rr;
1634     }
1635     clear_tail(d, opr_sz, simd_maxsz(desc));
1636 }
1637 
1638 /*
1639  * 64x64->128 polynomial multiply.
1640  * Because of the lanes are not accessed in strict columns,
1641  * this probably cannot be turned into a generic helper.
1642  */
1643 void HELPER(gvec_pmull_q)(void *vd, void *vn, void *vm, uint32_t desc)
1644 {
1645     intptr_t i, j, opr_sz = simd_oprsz(desc);
1646     intptr_t hi = simd_data(desc);
1647     uint64_t *d = vd, *n = vn, *m = vm;
1648 
1649     for (i = 0; i < opr_sz / 8; i += 2) {
1650         uint64_t nn = n[i + hi];
1651         uint64_t mm = m[i + hi];
1652         uint64_t rhi = 0;
1653         uint64_t rlo = 0;
1654 
1655         /* Bit 0 can only influence the low 64-bit result.  */
1656         if (nn & 1) {
1657             rlo = mm;
1658         }
1659 
1660         for (j = 1; j < 64; ++j) {
1661             uint64_t mask = -((nn >> j) & 1);
1662             rlo ^= (mm << j) & mask;
1663             rhi ^= (mm >> (64 - j)) & mask;
1664         }
1665         d[i] = rlo;
1666         d[i + 1] = rhi;
1667     }
1668     clear_tail(d, opr_sz, simd_maxsz(desc));
1669 }
1670 
1671 /*
1672  * 8x8->16 polynomial multiply.
1673  *
1674  * The byte inputs are expanded to (or extracted from) half-words.
1675  * Note that neon and sve2 get the inputs from different positions.
1676  * This allows 4 bytes to be processed in parallel with uint64_t.
1677  */
1678 
1679 static uint64_t expand_byte_to_half(uint64_t x)
1680 {
1681     return  (x & 0x000000ff)
1682          | ((x & 0x0000ff00) << 8)
1683          | ((x & 0x00ff0000) << 16)
1684          | ((x & 0xff000000) << 24);
1685 }
1686 
1687 static uint64_t pmull_h(uint64_t op1, uint64_t op2)
1688 {
1689     uint64_t result = 0;
1690     int i;
1691 
1692     for (i = 0; i < 8; ++i) {
1693         uint64_t mask = (op1 & 0x0001000100010001ull) * 0xffff;
1694         result ^= op2 & mask;
1695         op1 >>= 1;
1696         op2 <<= 1;
1697     }
1698     return result;
1699 }
1700 
1701 void HELPER(neon_pmull_h)(void *vd, void *vn, void *vm, uint32_t desc)
1702 {
1703     int hi = simd_data(desc);
1704     uint64_t *d = vd, *n = vn, *m = vm;
1705     uint64_t nn = n[hi], mm = m[hi];
1706 
1707     d[0] = pmull_h(expand_byte_to_half(nn), expand_byte_to_half(mm));
1708     nn >>= 32;
1709     mm >>= 32;
1710     d[1] = pmull_h(expand_byte_to_half(nn), expand_byte_to_half(mm));
1711 
1712     clear_tail(d, 16, simd_maxsz(desc));
1713 }
1714 
1715 #ifdef TARGET_AARCH64
1716 void HELPER(sve2_pmull_h)(void *vd, void *vn, void *vm, uint32_t desc)
1717 {
1718     int shift = simd_data(desc) * 8;
1719     intptr_t i, opr_sz = simd_oprsz(desc);
1720     uint64_t *d = vd, *n = vn, *m = vm;
1721 
1722     for (i = 0; i < opr_sz / 8; ++i) {
1723         uint64_t nn = (n[i] >> shift) & 0x00ff00ff00ff00ffull;
1724         uint64_t mm = (m[i] >> shift) & 0x00ff00ff00ff00ffull;
1725 
1726         d[i] = pmull_h(nn, mm);
1727     }
1728 }
1729 #endif
1730 
1731 #define DO_CMP0(NAME, TYPE, OP)                         \
1732 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1733 {                                                       \
1734     intptr_t i, opr_sz = simd_oprsz(desc);              \
1735     for (i = 0; i < opr_sz; i += sizeof(TYPE)) {        \
1736         TYPE nn = *(TYPE *)(vn + i);                    \
1737         *(TYPE *)(vd + i) = -(nn OP 0);                 \
1738     }                                                   \
1739     clear_tail(vd, opr_sz, simd_maxsz(desc));           \
1740 }
1741 
1742 DO_CMP0(gvec_ceq0_b, int8_t, ==)
1743 DO_CMP0(gvec_clt0_b, int8_t, <)
1744 DO_CMP0(gvec_cle0_b, int8_t, <=)
1745 DO_CMP0(gvec_cgt0_b, int8_t, >)
1746 DO_CMP0(gvec_cge0_b, int8_t, >=)
1747 
1748 DO_CMP0(gvec_ceq0_h, int16_t, ==)
1749 DO_CMP0(gvec_clt0_h, int16_t, <)
1750 DO_CMP0(gvec_cle0_h, int16_t, <=)
1751 DO_CMP0(gvec_cgt0_h, int16_t, >)
1752 DO_CMP0(gvec_cge0_h, int16_t, >=)
1753 
1754 #undef DO_CMP0
1755 
1756 #define DO_ABD(NAME, TYPE)                                      \
1757 void HELPER(NAME)(void *vd, void *vn, void *vm, uint32_t desc)  \
1758 {                                                               \
1759     intptr_t i, opr_sz = simd_oprsz(desc);                      \
1760     TYPE *d = vd, *n = vn, *m = vm;                             \
1761                                                                 \
1762     for (i = 0; i < opr_sz / sizeof(TYPE); ++i) {               \
1763         d[i] = n[i] < m[i] ? m[i] - n[i] : n[i] - m[i];         \
1764     }                                                           \
1765     clear_tail(d, opr_sz, simd_maxsz(desc));                    \
1766 }
1767 
1768 DO_ABD(gvec_sabd_b, int8_t)
1769 DO_ABD(gvec_sabd_h, int16_t)
1770 DO_ABD(gvec_sabd_s, int32_t)
1771 DO_ABD(gvec_sabd_d, int64_t)
1772 
1773 DO_ABD(gvec_uabd_b, uint8_t)
1774 DO_ABD(gvec_uabd_h, uint16_t)
1775 DO_ABD(gvec_uabd_s, uint32_t)
1776 DO_ABD(gvec_uabd_d, uint64_t)
1777 
1778 #undef DO_ABD
1779 
1780 #define DO_ABA(NAME, TYPE)                                      \
1781 void HELPER(NAME)(void *vd, void *vn, void *vm, uint32_t desc)  \
1782 {                                                               \
1783     intptr_t i, opr_sz = simd_oprsz(desc);                      \
1784     TYPE *d = vd, *n = vn, *m = vm;                             \
1785                                                                 \
1786     for (i = 0; i < opr_sz / sizeof(TYPE); ++i) {               \
1787         d[i] += n[i] < m[i] ? m[i] - n[i] : n[i] - m[i];        \
1788     }                                                           \
1789     clear_tail(d, opr_sz, simd_maxsz(desc));                    \
1790 }
1791 
1792 DO_ABA(gvec_saba_b, int8_t)
1793 DO_ABA(gvec_saba_h, int16_t)
1794 DO_ABA(gvec_saba_s, int32_t)
1795 DO_ABA(gvec_saba_d, int64_t)
1796 
1797 DO_ABA(gvec_uaba_b, uint8_t)
1798 DO_ABA(gvec_uaba_h, uint16_t)
1799 DO_ABA(gvec_uaba_s, uint32_t)
1800 DO_ABA(gvec_uaba_d, uint64_t)
1801 
1802 #undef DO_ABA
1803 
1804 #define DO_NEON_PAIRWISE(NAME, OP)                                      \
1805     void HELPER(NAME##s)(void *vd, void *vn, void *vm,                  \
1806                          void *stat, uint32_t oprsz)                    \
1807     {                                                                   \
1808         float_status *fpst = stat;                                      \
1809         float32 *d = vd;                                                \
1810         float32 *n = vn;                                                \
1811         float32 *m = vm;                                                \
1812         float32 r0, r1;                                                 \
1813                                                                         \
1814         /* Read all inputs before writing outputs in case vm == vd */   \
1815         r0 = float32_##OP(n[H4(0)], n[H4(1)], fpst);                    \
1816         r1 = float32_##OP(m[H4(0)], m[H4(1)], fpst);                    \
1817                                                                         \
1818         d[H4(0)] = r0;                                                  \
1819         d[H4(1)] = r1;                                                  \
1820     }                                                                   \
1821                                                                         \
1822     void HELPER(NAME##h)(void *vd, void *vn, void *vm,                  \
1823                          void *stat, uint32_t oprsz)                    \
1824     {                                                                   \
1825         float_status *fpst = stat;                                      \
1826         float16 *d = vd;                                                \
1827         float16 *n = vn;                                                \
1828         float16 *m = vm;                                                \
1829         float16 r0, r1, r2, r3;                                         \
1830                                                                         \
1831         /* Read all inputs before writing outputs in case vm == vd */   \
1832         r0 = float16_##OP(n[H2(0)], n[H2(1)], fpst);                    \
1833         r1 = float16_##OP(n[H2(2)], n[H2(3)], fpst);                    \
1834         r2 = float16_##OP(m[H2(0)], m[H2(1)], fpst);                    \
1835         r3 = float16_##OP(m[H2(2)], m[H2(3)], fpst);                    \
1836                                                                         \
1837         d[H4(0)] = r0;                                                  \
1838         d[H4(1)] = r1;                                                  \
1839         d[H4(2)] = r2;                                                  \
1840         d[H4(3)] = r3;                                                  \
1841     }
1842 
1843 DO_NEON_PAIRWISE(neon_padd, add)
1844 DO_NEON_PAIRWISE(neon_pmax, max)
1845 DO_NEON_PAIRWISE(neon_pmin, min)
1846 
1847 #undef DO_NEON_PAIRWISE
1848