xref: /qemu/target/arm/tcg/vec_helper.c (revision 1dc587ee9bfe804406eb3e0bacf47a80644d8abc)
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 #define DO_2OP(NAME, FUNC, TYPE) \
715 void HELPER(NAME)(void *vd, void *vn, void *stat, uint32_t desc)  \
716 {                                                                 \
717     intptr_t i, oprsz = simd_oprsz(desc);                         \
718     TYPE *d = vd, *n = vn;                                        \
719     for (i = 0; i < oprsz / sizeof(TYPE); i++) {                  \
720         d[i] = FUNC(n[i], stat);                                  \
721     }                                                             \
722     clear_tail(d, oprsz, simd_maxsz(desc));                       \
723 }
724 
725 DO_2OP(gvec_frecpe_h, helper_recpe_f16, float16)
726 DO_2OP(gvec_frecpe_s, helper_recpe_f32, float32)
727 DO_2OP(gvec_frecpe_d, helper_recpe_f64, float64)
728 
729 DO_2OP(gvec_frsqrte_h, helper_rsqrte_f16, float16)
730 DO_2OP(gvec_frsqrte_s, helper_rsqrte_f32, float32)
731 DO_2OP(gvec_frsqrte_d, helper_rsqrte_f64, float64)
732 
733 #define WRAP_CMP0_FWD(FN, CMPOP, TYPE)                          \
734     static TYPE TYPE##_##FN##0(TYPE op, float_status *stat)     \
735     {                                                           \
736         return TYPE##_##CMPOP(op, TYPE##_zero, stat);           \
737     }
738 
739 #define WRAP_CMP0_REV(FN, CMPOP, TYPE)                          \
740     static TYPE TYPE##_##FN##0(TYPE op, float_status *stat)    \
741     {                                                           \
742         return TYPE##_##CMPOP(TYPE##_zero, op, stat);           \
743     }
744 
745 #define DO_2OP_CMP0(FN, CMPOP, DIRN)                    \
746     WRAP_CMP0_##DIRN(FN, CMPOP, float16)                \
747     WRAP_CMP0_##DIRN(FN, CMPOP, float32)                \
748     DO_2OP(gvec_f##FN##0_h, float16_##FN##0, float16)   \
749     DO_2OP(gvec_f##FN##0_s, float32_##FN##0, float32)
750 
751 DO_2OP_CMP0(cgt, cgt, FWD)
752 DO_2OP_CMP0(cge, cge, FWD)
753 DO_2OP_CMP0(ceq, ceq, FWD)
754 DO_2OP_CMP0(clt, cgt, REV)
755 DO_2OP_CMP0(cle, cge, REV)
756 
757 #undef DO_2OP
758 #undef DO_2OP_CMP0
759 
760 /* Floating-point trigonometric starting value.
761  * See the ARM ARM pseudocode function FPTrigSMul.
762  */
763 static float16 float16_ftsmul(float16 op1, uint16_t op2, float_status *stat)
764 {
765     float16 result = float16_mul(op1, op1, stat);
766     if (!float16_is_any_nan(result)) {
767         result = float16_set_sign(result, op2 & 1);
768     }
769     return result;
770 }
771 
772 static float32 float32_ftsmul(float32 op1, uint32_t op2, float_status *stat)
773 {
774     float32 result = float32_mul(op1, op1, stat);
775     if (!float32_is_any_nan(result)) {
776         result = float32_set_sign(result, op2 & 1);
777     }
778     return result;
779 }
780 
781 static float64 float64_ftsmul(float64 op1, uint64_t op2, float_status *stat)
782 {
783     float64 result = float64_mul(op1, op1, stat);
784     if (!float64_is_any_nan(result)) {
785         result = float64_set_sign(result, op2 & 1);
786     }
787     return result;
788 }
789 
790 static float16 float16_abd(float16 op1, float16 op2, float_status *stat)
791 {
792     return float16_abs(float16_sub(op1, op2, stat));
793 }
794 
795 static float32 float32_abd(float32 op1, float32 op2, float_status *stat)
796 {
797     return float32_abs(float32_sub(op1, op2, stat));
798 }
799 
800 /*
801  * Reciprocal step. These are the AArch32 version which uses a
802  * non-fused multiply-and-subtract.
803  */
804 static float16 float16_recps_nf(float16 op1, float16 op2, float_status *stat)
805 {
806     op1 = float16_squash_input_denormal(op1, stat);
807     op2 = float16_squash_input_denormal(op2, stat);
808 
809     if ((float16_is_infinity(op1) && float16_is_zero(op2)) ||
810         (float16_is_infinity(op2) && float16_is_zero(op1))) {
811         return float16_two;
812     }
813     return float16_sub(float16_two, float16_mul(op1, op2, stat), stat);
814 }
815 
816 static float32 float32_recps_nf(float32 op1, float32 op2, float_status *stat)
817 {
818     op1 = float32_squash_input_denormal(op1, stat);
819     op2 = float32_squash_input_denormal(op2, stat);
820 
821     if ((float32_is_infinity(op1) && float32_is_zero(op2)) ||
822         (float32_is_infinity(op2) && float32_is_zero(op1))) {
823         return float32_two;
824     }
825     return float32_sub(float32_two, float32_mul(op1, op2, stat), stat);
826 }
827 
828 /* Reciprocal square-root step. AArch32 non-fused semantics. */
829 static float16 float16_rsqrts_nf(float16 op1, float16 op2, float_status *stat)
830 {
831     op1 = float16_squash_input_denormal(op1, stat);
832     op2 = float16_squash_input_denormal(op2, stat);
833 
834     if ((float16_is_infinity(op1) && float16_is_zero(op2)) ||
835         (float16_is_infinity(op2) && float16_is_zero(op1))) {
836         return float16_one_point_five;
837     }
838     op1 = float16_sub(float16_three, float16_mul(op1, op2, stat), stat);
839     return float16_div(op1, float16_two, stat);
840 }
841 
842 static float32 float32_rsqrts_nf(float32 op1, float32 op2, float_status *stat)
843 {
844     op1 = float32_squash_input_denormal(op1, stat);
845     op2 = float32_squash_input_denormal(op2, stat);
846 
847     if ((float32_is_infinity(op1) && float32_is_zero(op2)) ||
848         (float32_is_infinity(op2) && float32_is_zero(op1))) {
849         return float32_one_point_five;
850     }
851     op1 = float32_sub(float32_three, float32_mul(op1, op2, stat), stat);
852     return float32_div(op1, float32_two, stat);
853 }
854 
855 #define DO_3OP(NAME, FUNC, TYPE) \
856 void HELPER(NAME)(void *vd, void *vn, void *vm, void *stat, uint32_t desc) \
857 {                                                                          \
858     intptr_t i, oprsz = simd_oprsz(desc);                                  \
859     TYPE *d = vd, *n = vn, *m = vm;                                        \
860     for (i = 0; i < oprsz / sizeof(TYPE); i++) {                           \
861         d[i] = FUNC(n[i], m[i], stat);                                     \
862     }                                                                      \
863     clear_tail(d, oprsz, simd_maxsz(desc));                                \
864 }
865 
866 DO_3OP(gvec_fadd_h, float16_add, float16)
867 DO_3OP(gvec_fadd_s, float32_add, float32)
868 DO_3OP(gvec_fadd_d, float64_add, float64)
869 
870 DO_3OP(gvec_fsub_h, float16_sub, float16)
871 DO_3OP(gvec_fsub_s, float32_sub, float32)
872 DO_3OP(gvec_fsub_d, float64_sub, float64)
873 
874 DO_3OP(gvec_fmul_h, float16_mul, float16)
875 DO_3OP(gvec_fmul_s, float32_mul, float32)
876 DO_3OP(gvec_fmul_d, float64_mul, float64)
877 
878 DO_3OP(gvec_ftsmul_h, float16_ftsmul, float16)
879 DO_3OP(gvec_ftsmul_s, float32_ftsmul, float32)
880 DO_3OP(gvec_ftsmul_d, float64_ftsmul, float64)
881 
882 DO_3OP(gvec_fabd_h, float16_abd, float16)
883 DO_3OP(gvec_fabd_s, float32_abd, float32)
884 
885 DO_3OP(gvec_fceq_h, float16_ceq, float16)
886 DO_3OP(gvec_fceq_s, float32_ceq, float32)
887 
888 DO_3OP(gvec_fcge_h, float16_cge, float16)
889 DO_3OP(gvec_fcge_s, float32_cge, float32)
890 
891 DO_3OP(gvec_fcgt_h, float16_cgt, float16)
892 DO_3OP(gvec_fcgt_s, float32_cgt, float32)
893 
894 DO_3OP(gvec_facge_h, float16_acge, float16)
895 DO_3OP(gvec_facge_s, float32_acge, float32)
896 
897 DO_3OP(gvec_facgt_h, float16_acgt, float16)
898 DO_3OP(gvec_facgt_s, float32_acgt, float32)
899 
900 DO_3OP(gvec_fmax_h, float16_max, float16)
901 DO_3OP(gvec_fmax_s, float32_max, float32)
902 
903 DO_3OP(gvec_fmin_h, float16_min, float16)
904 DO_3OP(gvec_fmin_s, float32_min, float32)
905 
906 DO_3OP(gvec_fmaxnum_h, float16_maxnum, float16)
907 DO_3OP(gvec_fmaxnum_s, float32_maxnum, float32)
908 
909 DO_3OP(gvec_fminnum_h, float16_minnum, float16)
910 DO_3OP(gvec_fminnum_s, float32_minnum, float32)
911 
912 DO_3OP(gvec_recps_nf_h, float16_recps_nf, float16)
913 DO_3OP(gvec_recps_nf_s, float32_recps_nf, float32)
914 
915 DO_3OP(gvec_rsqrts_nf_h, float16_rsqrts_nf, float16)
916 DO_3OP(gvec_rsqrts_nf_s, float32_rsqrts_nf, float32)
917 
918 #ifdef TARGET_AARCH64
919 
920 DO_3OP(gvec_recps_h, helper_recpsf_f16, float16)
921 DO_3OP(gvec_recps_s, helper_recpsf_f32, float32)
922 DO_3OP(gvec_recps_d, helper_recpsf_f64, float64)
923 
924 DO_3OP(gvec_rsqrts_h, helper_rsqrtsf_f16, float16)
925 DO_3OP(gvec_rsqrts_s, helper_rsqrtsf_f32, float32)
926 DO_3OP(gvec_rsqrts_d, helper_rsqrtsf_f64, float64)
927 
928 #endif
929 #undef DO_3OP
930 
931 /* Non-fused multiply-add (unlike float16_muladd etc, which are fused) */
932 static float16 float16_muladd_nf(float16 dest, float16 op1, float16 op2,
933                                  float_status *stat)
934 {
935     return float16_add(dest, float16_mul(op1, op2, stat), stat);
936 }
937 
938 static float32 float32_muladd_nf(float32 dest, float32 op1, float32 op2,
939                                  float_status *stat)
940 {
941     return float32_add(dest, float32_mul(op1, op2, stat), stat);
942 }
943 
944 static float16 float16_mulsub_nf(float16 dest, float16 op1, float16 op2,
945                                  float_status *stat)
946 {
947     return float16_sub(dest, float16_mul(op1, op2, stat), stat);
948 }
949 
950 static float32 float32_mulsub_nf(float32 dest, float32 op1, float32 op2,
951                                  float_status *stat)
952 {
953     return float32_sub(dest, float32_mul(op1, op2, stat), stat);
954 }
955 
956 /* Fused versions; these have the semantics Neon VFMA/VFMS want */
957 static float16 float16_muladd_f(float16 dest, float16 op1, float16 op2,
958                                 float_status *stat)
959 {
960     return float16_muladd(op1, op2, dest, 0, stat);
961 }
962 
963 static float32 float32_muladd_f(float32 dest, float32 op1, float32 op2,
964                                  float_status *stat)
965 {
966     return float32_muladd(op1, op2, dest, 0, stat);
967 }
968 
969 static float16 float16_mulsub_f(float16 dest, float16 op1, float16 op2,
970                                  float_status *stat)
971 {
972     return float16_muladd(float16_chs(op1), op2, dest, 0, stat);
973 }
974 
975 static float32 float32_mulsub_f(float32 dest, float32 op1, float32 op2,
976                                  float_status *stat)
977 {
978     return float32_muladd(float32_chs(op1), op2, dest, 0, stat);
979 }
980 
981 #define DO_MULADD(NAME, FUNC, TYPE)                                     \
982 void HELPER(NAME)(void *vd, void *vn, void *vm, void *stat, uint32_t desc) \
983 {                                                                          \
984     intptr_t i, oprsz = simd_oprsz(desc);                                  \
985     TYPE *d = vd, *n = vn, *m = vm;                                        \
986     for (i = 0; i < oprsz / sizeof(TYPE); i++) {                           \
987         d[i] = FUNC(d[i], n[i], m[i], stat);                               \
988     }                                                                      \
989     clear_tail(d, oprsz, simd_maxsz(desc));                                \
990 }
991 
992 DO_MULADD(gvec_fmla_h, float16_muladd_nf, float16)
993 DO_MULADD(gvec_fmla_s, float32_muladd_nf, float32)
994 
995 DO_MULADD(gvec_fmls_h, float16_mulsub_nf, float16)
996 DO_MULADD(gvec_fmls_s, float32_mulsub_nf, float32)
997 
998 DO_MULADD(gvec_vfma_h, float16_muladd_f, float16)
999 DO_MULADD(gvec_vfma_s, float32_muladd_f, float32)
1000 
1001 DO_MULADD(gvec_vfms_h, float16_mulsub_f, float16)
1002 DO_MULADD(gvec_vfms_s, float32_mulsub_f, float32)
1003 
1004 /* For the indexed ops, SVE applies the index per 128-bit vector segment.
1005  * For AdvSIMD, there is of course only one such vector segment.
1006  */
1007 
1008 #define DO_MUL_IDX(NAME, TYPE, H) \
1009 void HELPER(NAME)(void *vd, void *vn, void *vm, uint32_t desc) \
1010 {                                                                          \
1011     intptr_t i, j, oprsz = simd_oprsz(desc), segment = 16 / sizeof(TYPE);  \
1012     intptr_t idx = simd_data(desc);                                        \
1013     TYPE *d = vd, *n = vn, *m = vm;                                        \
1014     for (i = 0; i < oprsz / sizeof(TYPE); i += segment) {                  \
1015         TYPE mm = m[H(i + idx)];                                           \
1016         for (j = 0; j < segment; j++) {                                    \
1017             d[i + j] = n[i + j] * mm;                                      \
1018         }                                                                  \
1019     }                                                                      \
1020     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1021 }
1022 
1023 DO_MUL_IDX(gvec_mul_idx_h, uint16_t, H2)
1024 DO_MUL_IDX(gvec_mul_idx_s, uint32_t, H4)
1025 DO_MUL_IDX(gvec_mul_idx_d, uint64_t, )
1026 
1027 #undef DO_MUL_IDX
1028 
1029 #define DO_MLA_IDX(NAME, TYPE, OP, H) \
1030 void HELPER(NAME)(void *vd, void *vn, void *vm, void *va, uint32_t desc)   \
1031 {                                                                          \
1032     intptr_t i, j, oprsz = simd_oprsz(desc), segment = 16 / sizeof(TYPE);  \
1033     intptr_t idx = simd_data(desc);                                        \
1034     TYPE *d = vd, *n = vn, *m = vm, *a = va;                               \
1035     for (i = 0; i < oprsz / sizeof(TYPE); i += segment) {                  \
1036         TYPE mm = m[H(i + idx)];                                           \
1037         for (j = 0; j < segment; j++) {                                    \
1038             d[i + j] = a[i + j] OP n[i + j] * mm;                          \
1039         }                                                                  \
1040     }                                                                      \
1041     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1042 }
1043 
1044 DO_MLA_IDX(gvec_mla_idx_h, uint16_t, +, H2)
1045 DO_MLA_IDX(gvec_mla_idx_s, uint32_t, +, H4)
1046 DO_MLA_IDX(gvec_mla_idx_d, uint64_t, +,   )
1047 
1048 DO_MLA_IDX(gvec_mls_idx_h, uint16_t, -, H2)
1049 DO_MLA_IDX(gvec_mls_idx_s, uint32_t, -, H4)
1050 DO_MLA_IDX(gvec_mls_idx_d, uint64_t, -,   )
1051 
1052 #undef DO_MLA_IDX
1053 
1054 #define DO_FMUL_IDX(NAME, TYPE, H) \
1055 void HELPER(NAME)(void *vd, void *vn, void *vm, void *stat, uint32_t desc) \
1056 {                                                                          \
1057     intptr_t i, j, oprsz = simd_oprsz(desc), segment = 16 / sizeof(TYPE);  \
1058     intptr_t idx = simd_data(desc);                                        \
1059     TYPE *d = vd, *n = vn, *m = vm;                                        \
1060     for (i = 0; i < oprsz / sizeof(TYPE); i += segment) {                  \
1061         TYPE mm = m[H(i + idx)];                                           \
1062         for (j = 0; j < segment; j++) {                                    \
1063             d[i + j] = TYPE##_mul(n[i + j], mm, stat);                     \
1064         }                                                                  \
1065     }                                                                      \
1066     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1067 }
1068 
1069 DO_FMUL_IDX(gvec_fmul_idx_h, float16, H2)
1070 DO_FMUL_IDX(gvec_fmul_idx_s, float32, H4)
1071 DO_FMUL_IDX(gvec_fmul_idx_d, float64, )
1072 
1073 #undef DO_FMUL_IDX
1074 
1075 #define DO_FMLA_IDX(NAME, TYPE, H)                                         \
1076 void HELPER(NAME)(void *vd, void *vn, void *vm, void *va,                  \
1077                   void *stat, uint32_t desc)                               \
1078 {                                                                          \
1079     intptr_t i, j, oprsz = simd_oprsz(desc), segment = 16 / sizeof(TYPE);  \
1080     TYPE op1_neg = extract32(desc, SIMD_DATA_SHIFT, 1);                    \
1081     intptr_t idx = desc >> (SIMD_DATA_SHIFT + 1);                          \
1082     TYPE *d = vd, *n = vn, *m = vm, *a = va;                               \
1083     op1_neg <<= (8 * sizeof(TYPE) - 1);                                    \
1084     for (i = 0; i < oprsz / sizeof(TYPE); i += segment) {                  \
1085         TYPE mm = m[H(i + idx)];                                           \
1086         for (j = 0; j < segment; j++) {                                    \
1087             d[i + j] = TYPE##_muladd(n[i + j] ^ op1_neg,                   \
1088                                      mm, a[i + j], 0, stat);               \
1089         }                                                                  \
1090     }                                                                      \
1091     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1092 }
1093 
1094 DO_FMLA_IDX(gvec_fmla_idx_h, float16, H2)
1095 DO_FMLA_IDX(gvec_fmla_idx_s, float32, H4)
1096 DO_FMLA_IDX(gvec_fmla_idx_d, float64, )
1097 
1098 #undef DO_FMLA_IDX
1099 
1100 #define DO_SAT(NAME, WTYPE, TYPEN, TYPEM, OP, MIN, MAX) \
1101 void HELPER(NAME)(void *vd, void *vq, void *vn, void *vm, uint32_t desc)   \
1102 {                                                                          \
1103     intptr_t i, oprsz = simd_oprsz(desc);                                  \
1104     TYPEN *d = vd, *n = vn; TYPEM *m = vm;                                 \
1105     bool q = false;                                                        \
1106     for (i = 0; i < oprsz / sizeof(TYPEN); i++) {                          \
1107         WTYPE dd = (WTYPE)n[i] OP m[i];                                    \
1108         if (dd < MIN) {                                                    \
1109             dd = MIN;                                                      \
1110             q = true;                                                      \
1111         } else if (dd > MAX) {                                             \
1112             dd = MAX;                                                      \
1113             q = true;                                                      \
1114         }                                                                  \
1115         d[i] = dd;                                                         \
1116     }                                                                      \
1117     if (q) {                                                               \
1118         uint32_t *qc = vq;                                                 \
1119         qc[0] = 1;                                                         \
1120     }                                                                      \
1121     clear_tail(d, oprsz, simd_maxsz(desc));                                \
1122 }
1123 
1124 DO_SAT(gvec_uqadd_b, int, uint8_t, uint8_t, +, 0, UINT8_MAX)
1125 DO_SAT(gvec_uqadd_h, int, uint16_t, uint16_t, +, 0, UINT16_MAX)
1126 DO_SAT(gvec_uqadd_s, int64_t, uint32_t, uint32_t, +, 0, UINT32_MAX)
1127 
1128 DO_SAT(gvec_sqadd_b, int, int8_t, int8_t, +, INT8_MIN, INT8_MAX)
1129 DO_SAT(gvec_sqadd_h, int, int16_t, int16_t, +, INT16_MIN, INT16_MAX)
1130 DO_SAT(gvec_sqadd_s, int64_t, int32_t, int32_t, +, INT32_MIN, INT32_MAX)
1131 
1132 DO_SAT(gvec_uqsub_b, int, uint8_t, uint8_t, -, 0, UINT8_MAX)
1133 DO_SAT(gvec_uqsub_h, int, uint16_t, uint16_t, -, 0, UINT16_MAX)
1134 DO_SAT(gvec_uqsub_s, int64_t, uint32_t, uint32_t, -, 0, UINT32_MAX)
1135 
1136 DO_SAT(gvec_sqsub_b, int, int8_t, int8_t, -, INT8_MIN, INT8_MAX)
1137 DO_SAT(gvec_sqsub_h, int, int16_t, int16_t, -, INT16_MIN, INT16_MAX)
1138 DO_SAT(gvec_sqsub_s, int64_t, int32_t, int32_t, -, INT32_MIN, INT32_MAX)
1139 
1140 #undef DO_SAT
1141 
1142 void HELPER(gvec_uqadd_d)(void *vd, void *vq, void *vn,
1143                           void *vm, uint32_t desc)
1144 {
1145     intptr_t i, oprsz = simd_oprsz(desc);
1146     uint64_t *d = vd, *n = vn, *m = vm;
1147     bool q = false;
1148 
1149     for (i = 0; i < oprsz / 8; i++) {
1150         uint64_t nn = n[i], mm = m[i], dd = nn + mm;
1151         if (dd < nn) {
1152             dd = UINT64_MAX;
1153             q = true;
1154         }
1155         d[i] = dd;
1156     }
1157     if (q) {
1158         uint32_t *qc = vq;
1159         qc[0] = 1;
1160     }
1161     clear_tail(d, oprsz, simd_maxsz(desc));
1162 }
1163 
1164 void HELPER(gvec_uqsub_d)(void *vd, void *vq, void *vn,
1165                           void *vm, uint32_t desc)
1166 {
1167     intptr_t i, oprsz = simd_oprsz(desc);
1168     uint64_t *d = vd, *n = vn, *m = vm;
1169     bool q = false;
1170 
1171     for (i = 0; i < oprsz / 8; i++) {
1172         uint64_t nn = n[i], mm = m[i], dd = nn - mm;
1173         if (nn < mm) {
1174             dd = 0;
1175             q = true;
1176         }
1177         d[i] = dd;
1178     }
1179     if (q) {
1180         uint32_t *qc = vq;
1181         qc[0] = 1;
1182     }
1183     clear_tail(d, oprsz, simd_maxsz(desc));
1184 }
1185 
1186 void HELPER(gvec_sqadd_d)(void *vd, void *vq, void *vn,
1187                           void *vm, uint32_t desc)
1188 {
1189     intptr_t i, oprsz = simd_oprsz(desc);
1190     int64_t *d = vd, *n = vn, *m = vm;
1191     bool q = false;
1192 
1193     for (i = 0; i < oprsz / 8; i++) {
1194         int64_t nn = n[i], mm = m[i], dd = nn + mm;
1195         if (((dd ^ nn) & ~(nn ^ mm)) & INT64_MIN) {
1196             dd = (nn >> 63) ^ ~INT64_MIN;
1197             q = true;
1198         }
1199         d[i] = dd;
1200     }
1201     if (q) {
1202         uint32_t *qc = vq;
1203         qc[0] = 1;
1204     }
1205     clear_tail(d, oprsz, simd_maxsz(desc));
1206 }
1207 
1208 void HELPER(gvec_sqsub_d)(void *vd, void *vq, void *vn,
1209                           void *vm, uint32_t desc)
1210 {
1211     intptr_t i, oprsz = simd_oprsz(desc);
1212     int64_t *d = vd, *n = vn, *m = vm;
1213     bool q = false;
1214 
1215     for (i = 0; i < oprsz / 8; i++) {
1216         int64_t nn = n[i], mm = m[i], dd = nn - mm;
1217         if (((dd ^ nn) & (nn ^ mm)) & INT64_MIN) {
1218             dd = (nn >> 63) ^ ~INT64_MIN;
1219             q = true;
1220         }
1221         d[i] = dd;
1222     }
1223     if (q) {
1224         uint32_t *qc = vq;
1225         qc[0] = 1;
1226     }
1227     clear_tail(d, oprsz, simd_maxsz(desc));
1228 }
1229 
1230 
1231 #define DO_SRA(NAME, TYPE)                              \
1232 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1233 {                                                       \
1234     intptr_t i, oprsz = simd_oprsz(desc);               \
1235     int shift = simd_data(desc);                        \
1236     TYPE *d = vd, *n = vn;                              \
1237     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1238         d[i] += n[i] >> shift;                          \
1239     }                                                   \
1240     clear_tail(d, oprsz, simd_maxsz(desc));             \
1241 }
1242 
1243 DO_SRA(gvec_ssra_b, int8_t)
1244 DO_SRA(gvec_ssra_h, int16_t)
1245 DO_SRA(gvec_ssra_s, int32_t)
1246 DO_SRA(gvec_ssra_d, int64_t)
1247 
1248 DO_SRA(gvec_usra_b, uint8_t)
1249 DO_SRA(gvec_usra_h, uint16_t)
1250 DO_SRA(gvec_usra_s, uint32_t)
1251 DO_SRA(gvec_usra_d, uint64_t)
1252 
1253 #undef DO_SRA
1254 
1255 #define DO_RSHR(NAME, TYPE)                             \
1256 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1257 {                                                       \
1258     intptr_t i, oprsz = simd_oprsz(desc);               \
1259     int shift = simd_data(desc);                        \
1260     TYPE *d = vd, *n = vn;                              \
1261     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1262         TYPE tmp = n[i] >> (shift - 1);                 \
1263         d[i] = (tmp >> 1) + (tmp & 1);                  \
1264     }                                                   \
1265     clear_tail(d, oprsz, simd_maxsz(desc));             \
1266 }
1267 
1268 DO_RSHR(gvec_srshr_b, int8_t)
1269 DO_RSHR(gvec_srshr_h, int16_t)
1270 DO_RSHR(gvec_srshr_s, int32_t)
1271 DO_RSHR(gvec_srshr_d, int64_t)
1272 
1273 DO_RSHR(gvec_urshr_b, uint8_t)
1274 DO_RSHR(gvec_urshr_h, uint16_t)
1275 DO_RSHR(gvec_urshr_s, uint32_t)
1276 DO_RSHR(gvec_urshr_d, uint64_t)
1277 
1278 #undef DO_RSHR
1279 
1280 #define DO_RSRA(NAME, TYPE)                             \
1281 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1282 {                                                       \
1283     intptr_t i, oprsz = simd_oprsz(desc);               \
1284     int shift = simd_data(desc);                        \
1285     TYPE *d = vd, *n = vn;                              \
1286     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1287         TYPE tmp = n[i] >> (shift - 1);                 \
1288         d[i] += (tmp >> 1) + (tmp & 1);                 \
1289     }                                                   \
1290     clear_tail(d, oprsz, simd_maxsz(desc));             \
1291 }
1292 
1293 DO_RSRA(gvec_srsra_b, int8_t)
1294 DO_RSRA(gvec_srsra_h, int16_t)
1295 DO_RSRA(gvec_srsra_s, int32_t)
1296 DO_RSRA(gvec_srsra_d, int64_t)
1297 
1298 DO_RSRA(gvec_ursra_b, uint8_t)
1299 DO_RSRA(gvec_ursra_h, uint16_t)
1300 DO_RSRA(gvec_ursra_s, uint32_t)
1301 DO_RSRA(gvec_ursra_d, uint64_t)
1302 
1303 #undef DO_RSRA
1304 
1305 #define DO_SRI(NAME, TYPE)                              \
1306 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1307 {                                                       \
1308     intptr_t i, oprsz = simd_oprsz(desc);               \
1309     int shift = simd_data(desc);                        \
1310     TYPE *d = vd, *n = vn;                              \
1311     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1312         d[i] = deposit64(d[i], 0, sizeof(TYPE) * 8 - shift, n[i] >> shift); \
1313     }                                                   \
1314     clear_tail(d, oprsz, simd_maxsz(desc));             \
1315 }
1316 
1317 DO_SRI(gvec_sri_b, uint8_t)
1318 DO_SRI(gvec_sri_h, uint16_t)
1319 DO_SRI(gvec_sri_s, uint32_t)
1320 DO_SRI(gvec_sri_d, uint64_t)
1321 
1322 #undef DO_SRI
1323 
1324 #define DO_SLI(NAME, TYPE)                              \
1325 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1326 {                                                       \
1327     intptr_t i, oprsz = simd_oprsz(desc);               \
1328     int shift = simd_data(desc);                        \
1329     TYPE *d = vd, *n = vn;                              \
1330     for (i = 0; i < oprsz / sizeof(TYPE); i++) {        \
1331         d[i] = deposit64(d[i], shift, sizeof(TYPE) * 8 - shift, n[i]); \
1332     }                                                   \
1333     clear_tail(d, oprsz, simd_maxsz(desc));             \
1334 }
1335 
1336 DO_SLI(gvec_sli_b, uint8_t)
1337 DO_SLI(gvec_sli_h, uint16_t)
1338 DO_SLI(gvec_sli_s, uint32_t)
1339 DO_SLI(gvec_sli_d, uint64_t)
1340 
1341 #undef DO_SLI
1342 
1343 /*
1344  * Convert float16 to float32, raising no exceptions and
1345  * preserving exceptional values, including SNaN.
1346  * This is effectively an unpack+repack operation.
1347  */
1348 static float32 float16_to_float32_by_bits(uint32_t f16, bool fz16)
1349 {
1350     const int f16_bias = 15;
1351     const int f32_bias = 127;
1352     uint32_t sign = extract32(f16, 15, 1);
1353     uint32_t exp = extract32(f16, 10, 5);
1354     uint32_t frac = extract32(f16, 0, 10);
1355 
1356     if (exp == 0x1f) {
1357         /* Inf or NaN */
1358         exp = 0xff;
1359     } else if (exp == 0) {
1360         /* Zero or denormal.  */
1361         if (frac != 0) {
1362             if (fz16) {
1363                 frac = 0;
1364             } else {
1365                 /*
1366                  * Denormal; these are all normal float32.
1367                  * Shift the fraction so that the msb is at bit 11,
1368                  * then remove bit 11 as the implicit bit of the
1369                  * normalized float32.  Note that we still go through
1370                  * the shift for normal numbers below, to put the
1371                  * float32 fraction at the right place.
1372                  */
1373                 int shift = clz32(frac) - 21;
1374                 frac = (frac << shift) & 0x3ff;
1375                 exp = f32_bias - f16_bias - shift + 1;
1376             }
1377         }
1378     } else {
1379         /* Normal number; adjust the bias.  */
1380         exp += f32_bias - f16_bias;
1381     }
1382     sign <<= 31;
1383     exp <<= 23;
1384     frac <<= 23 - 10;
1385 
1386     return sign | exp | frac;
1387 }
1388 
1389 static uint64_t load4_f16(uint64_t *ptr, int is_q, int is_2)
1390 {
1391     /*
1392      * Branchless load of u32[0], u64[0], u32[1], or u64[1].
1393      * Load the 2nd qword iff is_q & is_2.
1394      * Shift to the 2nd dword iff !is_q & is_2.
1395      * For !is_q & !is_2, the upper bits of the result are garbage.
1396      */
1397     return ptr[is_q & is_2] >> ((is_2 & ~is_q) << 5);
1398 }
1399 
1400 /*
1401  * Note that FMLAL requires oprsz == 8 or oprsz == 16,
1402  * as there is not yet SVE versions that might use blocking.
1403  */
1404 
1405 static void do_fmlal(float32 *d, void *vn, void *vm, float_status *fpst,
1406                      uint32_t desc, bool fz16)
1407 {
1408     intptr_t i, oprsz = simd_oprsz(desc);
1409     int is_s = extract32(desc, SIMD_DATA_SHIFT, 1);
1410     int is_2 = extract32(desc, SIMD_DATA_SHIFT + 1, 1);
1411     int is_q = oprsz == 16;
1412     uint64_t n_4, m_4;
1413 
1414     /* Pre-load all of the f16 data, avoiding overlap issues.  */
1415     n_4 = load4_f16(vn, is_q, is_2);
1416     m_4 = load4_f16(vm, is_q, is_2);
1417 
1418     /* Negate all inputs for FMLSL at once.  */
1419     if (is_s) {
1420         n_4 ^= 0x8000800080008000ull;
1421     }
1422 
1423     for (i = 0; i < oprsz / 4; i++) {
1424         float32 n_1 = float16_to_float32_by_bits(n_4 >> (i * 16), fz16);
1425         float32 m_1 = float16_to_float32_by_bits(m_4 >> (i * 16), fz16);
1426         d[H4(i)] = float32_muladd(n_1, m_1, d[H4(i)], 0, fpst);
1427     }
1428     clear_tail(d, oprsz, simd_maxsz(desc));
1429 }
1430 
1431 void HELPER(gvec_fmlal_a32)(void *vd, void *vn, void *vm,
1432                             void *venv, uint32_t desc)
1433 {
1434     CPUARMState *env = venv;
1435     do_fmlal(vd, vn, vm, &env->vfp.standard_fp_status, desc,
1436              get_flush_inputs_to_zero(&env->vfp.fp_status_f16));
1437 }
1438 
1439 void HELPER(gvec_fmlal_a64)(void *vd, void *vn, void *vm,
1440                             void *venv, uint32_t desc)
1441 {
1442     CPUARMState *env = venv;
1443     do_fmlal(vd, vn, vm, &env->vfp.fp_status, desc,
1444              get_flush_inputs_to_zero(&env->vfp.fp_status_f16));
1445 }
1446 
1447 static void do_fmlal_idx(float32 *d, void *vn, void *vm, float_status *fpst,
1448                          uint32_t desc, bool fz16)
1449 {
1450     intptr_t i, oprsz = simd_oprsz(desc);
1451     int is_s = extract32(desc, SIMD_DATA_SHIFT, 1);
1452     int is_2 = extract32(desc, SIMD_DATA_SHIFT + 1, 1);
1453     int index = extract32(desc, SIMD_DATA_SHIFT + 2, 3);
1454     int is_q = oprsz == 16;
1455     uint64_t n_4;
1456     float32 m_1;
1457 
1458     /* Pre-load all of the f16 data, avoiding overlap issues.  */
1459     n_4 = load4_f16(vn, is_q, is_2);
1460 
1461     /* Negate all inputs for FMLSL at once.  */
1462     if (is_s) {
1463         n_4 ^= 0x8000800080008000ull;
1464     }
1465 
1466     m_1 = float16_to_float32_by_bits(((float16 *)vm)[H2(index)], fz16);
1467 
1468     for (i = 0; i < oprsz / 4; i++) {
1469         float32 n_1 = float16_to_float32_by_bits(n_4 >> (i * 16), fz16);
1470         d[H4(i)] = float32_muladd(n_1, m_1, d[H4(i)], 0, fpst);
1471     }
1472     clear_tail(d, oprsz, simd_maxsz(desc));
1473 }
1474 
1475 void HELPER(gvec_fmlal_idx_a32)(void *vd, void *vn, void *vm,
1476                                 void *venv, uint32_t desc)
1477 {
1478     CPUARMState *env = venv;
1479     do_fmlal_idx(vd, vn, vm, &env->vfp.standard_fp_status, desc,
1480                  get_flush_inputs_to_zero(&env->vfp.fp_status_f16));
1481 }
1482 
1483 void HELPER(gvec_fmlal_idx_a64)(void *vd, void *vn, void *vm,
1484                                 void *venv, uint32_t desc)
1485 {
1486     CPUARMState *env = venv;
1487     do_fmlal_idx(vd, vn, vm, &env->vfp.fp_status, desc,
1488                  get_flush_inputs_to_zero(&env->vfp.fp_status_f16));
1489 }
1490 
1491 void HELPER(gvec_sshl_b)(void *vd, void *vn, void *vm, uint32_t desc)
1492 {
1493     intptr_t i, opr_sz = simd_oprsz(desc);
1494     int8_t *d = vd, *n = vn, *m = vm;
1495 
1496     for (i = 0; i < opr_sz; ++i) {
1497         int8_t mm = m[i];
1498         int8_t nn = n[i];
1499         int8_t res = 0;
1500         if (mm >= 0) {
1501             if (mm < 8) {
1502                 res = nn << mm;
1503             }
1504         } else {
1505             res = nn >> (mm > -8 ? -mm : 7);
1506         }
1507         d[i] = res;
1508     }
1509     clear_tail(d, opr_sz, simd_maxsz(desc));
1510 }
1511 
1512 void HELPER(gvec_sshl_h)(void *vd, void *vn, void *vm, uint32_t desc)
1513 {
1514     intptr_t i, opr_sz = simd_oprsz(desc);
1515     int16_t *d = vd, *n = vn, *m = vm;
1516 
1517     for (i = 0; i < opr_sz / 2; ++i) {
1518         int8_t mm = m[i];   /* only 8 bits of shift are significant */
1519         int16_t nn = n[i];
1520         int16_t res = 0;
1521         if (mm >= 0) {
1522             if (mm < 16) {
1523                 res = nn << mm;
1524             }
1525         } else {
1526             res = nn >> (mm > -16 ? -mm : 15);
1527         }
1528         d[i] = res;
1529     }
1530     clear_tail(d, opr_sz, simd_maxsz(desc));
1531 }
1532 
1533 void HELPER(gvec_ushl_b)(void *vd, void *vn, void *vm, uint32_t desc)
1534 {
1535     intptr_t i, opr_sz = simd_oprsz(desc);
1536     uint8_t *d = vd, *n = vn, *m = vm;
1537 
1538     for (i = 0; i < opr_sz; ++i) {
1539         int8_t mm = m[i];
1540         uint8_t nn = n[i];
1541         uint8_t res = 0;
1542         if (mm >= 0) {
1543             if (mm < 8) {
1544                 res = nn << mm;
1545             }
1546         } else {
1547             if (mm > -8) {
1548                 res = nn >> -mm;
1549             }
1550         }
1551         d[i] = res;
1552     }
1553     clear_tail(d, opr_sz, simd_maxsz(desc));
1554 }
1555 
1556 void HELPER(gvec_ushl_h)(void *vd, void *vn, void *vm, uint32_t desc)
1557 {
1558     intptr_t i, opr_sz = simd_oprsz(desc);
1559     uint16_t *d = vd, *n = vn, *m = vm;
1560 
1561     for (i = 0; i < opr_sz / 2; ++i) {
1562         int8_t mm = m[i];   /* only 8 bits of shift are significant */
1563         uint16_t nn = n[i];
1564         uint16_t res = 0;
1565         if (mm >= 0) {
1566             if (mm < 16) {
1567                 res = nn << mm;
1568             }
1569         } else {
1570             if (mm > -16) {
1571                 res = nn >> -mm;
1572             }
1573         }
1574         d[i] = res;
1575     }
1576     clear_tail(d, opr_sz, simd_maxsz(desc));
1577 }
1578 
1579 /*
1580  * 8x8->8 polynomial multiply.
1581  *
1582  * Polynomial multiplication is like integer multiplication except the
1583  * partial products are XORed, not added.
1584  *
1585  * TODO: expose this as a generic vector operation, as it is a common
1586  * crypto building block.
1587  */
1588 void HELPER(gvec_pmul_b)(void *vd, void *vn, void *vm, uint32_t desc)
1589 {
1590     intptr_t i, j, opr_sz = simd_oprsz(desc);
1591     uint64_t *d = vd, *n = vn, *m = vm;
1592 
1593     for (i = 0; i < opr_sz / 8; ++i) {
1594         uint64_t nn = n[i];
1595         uint64_t mm = m[i];
1596         uint64_t rr = 0;
1597 
1598         for (j = 0; j < 8; ++j) {
1599             uint64_t mask = (nn & 0x0101010101010101ull) * 0xff;
1600             rr ^= mm & mask;
1601             mm = (mm << 1) & 0xfefefefefefefefeull;
1602             nn >>= 1;
1603         }
1604         d[i] = rr;
1605     }
1606     clear_tail(d, opr_sz, simd_maxsz(desc));
1607 }
1608 
1609 /*
1610  * 64x64->128 polynomial multiply.
1611  * Because of the lanes are not accessed in strict columns,
1612  * this probably cannot be turned into a generic helper.
1613  */
1614 void HELPER(gvec_pmull_q)(void *vd, void *vn, void *vm, uint32_t desc)
1615 {
1616     intptr_t i, j, opr_sz = simd_oprsz(desc);
1617     intptr_t hi = simd_data(desc);
1618     uint64_t *d = vd, *n = vn, *m = vm;
1619 
1620     for (i = 0; i < opr_sz / 8; i += 2) {
1621         uint64_t nn = n[i + hi];
1622         uint64_t mm = m[i + hi];
1623         uint64_t rhi = 0;
1624         uint64_t rlo = 0;
1625 
1626         /* Bit 0 can only influence the low 64-bit result.  */
1627         if (nn & 1) {
1628             rlo = mm;
1629         }
1630 
1631         for (j = 1; j < 64; ++j) {
1632             uint64_t mask = -((nn >> j) & 1);
1633             rlo ^= (mm << j) & mask;
1634             rhi ^= (mm >> (64 - j)) & mask;
1635         }
1636         d[i] = rlo;
1637         d[i + 1] = rhi;
1638     }
1639     clear_tail(d, opr_sz, simd_maxsz(desc));
1640 }
1641 
1642 /*
1643  * 8x8->16 polynomial multiply.
1644  *
1645  * The byte inputs are expanded to (or extracted from) half-words.
1646  * Note that neon and sve2 get the inputs from different positions.
1647  * This allows 4 bytes to be processed in parallel with uint64_t.
1648  */
1649 
1650 static uint64_t expand_byte_to_half(uint64_t x)
1651 {
1652     return  (x & 0x000000ff)
1653          | ((x & 0x0000ff00) << 8)
1654          | ((x & 0x00ff0000) << 16)
1655          | ((x & 0xff000000) << 24);
1656 }
1657 
1658 static uint64_t pmull_h(uint64_t op1, uint64_t op2)
1659 {
1660     uint64_t result = 0;
1661     int i;
1662 
1663     for (i = 0; i < 8; ++i) {
1664         uint64_t mask = (op1 & 0x0001000100010001ull) * 0xffff;
1665         result ^= op2 & mask;
1666         op1 >>= 1;
1667         op2 <<= 1;
1668     }
1669     return result;
1670 }
1671 
1672 void HELPER(neon_pmull_h)(void *vd, void *vn, void *vm, uint32_t desc)
1673 {
1674     int hi = simd_data(desc);
1675     uint64_t *d = vd, *n = vn, *m = vm;
1676     uint64_t nn = n[hi], mm = m[hi];
1677 
1678     d[0] = pmull_h(expand_byte_to_half(nn), expand_byte_to_half(mm));
1679     nn >>= 32;
1680     mm >>= 32;
1681     d[1] = pmull_h(expand_byte_to_half(nn), expand_byte_to_half(mm));
1682 
1683     clear_tail(d, 16, simd_maxsz(desc));
1684 }
1685 
1686 #ifdef TARGET_AARCH64
1687 void HELPER(sve2_pmull_h)(void *vd, void *vn, void *vm, uint32_t desc)
1688 {
1689     int shift = simd_data(desc) * 8;
1690     intptr_t i, opr_sz = simd_oprsz(desc);
1691     uint64_t *d = vd, *n = vn, *m = vm;
1692 
1693     for (i = 0; i < opr_sz / 8; ++i) {
1694         uint64_t nn = (n[i] >> shift) & 0x00ff00ff00ff00ffull;
1695         uint64_t mm = (m[i] >> shift) & 0x00ff00ff00ff00ffull;
1696 
1697         d[i] = pmull_h(nn, mm);
1698     }
1699 }
1700 #endif
1701 
1702 #define DO_CMP0(NAME, TYPE, OP)                         \
1703 void HELPER(NAME)(void *vd, void *vn, uint32_t desc)    \
1704 {                                                       \
1705     intptr_t i, opr_sz = simd_oprsz(desc);              \
1706     for (i = 0; i < opr_sz; i += sizeof(TYPE)) {        \
1707         TYPE nn = *(TYPE *)(vn + i);                    \
1708         *(TYPE *)(vd + i) = -(nn OP 0);                 \
1709     }                                                   \
1710     clear_tail(vd, opr_sz, simd_maxsz(desc));           \
1711 }
1712 
1713 DO_CMP0(gvec_ceq0_b, int8_t, ==)
1714 DO_CMP0(gvec_clt0_b, int8_t, <)
1715 DO_CMP0(gvec_cle0_b, int8_t, <=)
1716 DO_CMP0(gvec_cgt0_b, int8_t, >)
1717 DO_CMP0(gvec_cge0_b, int8_t, >=)
1718 
1719 DO_CMP0(gvec_ceq0_h, int16_t, ==)
1720 DO_CMP0(gvec_clt0_h, int16_t, <)
1721 DO_CMP0(gvec_cle0_h, int16_t, <=)
1722 DO_CMP0(gvec_cgt0_h, int16_t, >)
1723 DO_CMP0(gvec_cge0_h, int16_t, >=)
1724 
1725 #undef DO_CMP0
1726 
1727 #define DO_ABD(NAME, TYPE)                                      \
1728 void HELPER(NAME)(void *vd, void *vn, void *vm, uint32_t desc)  \
1729 {                                                               \
1730     intptr_t i, opr_sz = simd_oprsz(desc);                      \
1731     TYPE *d = vd, *n = vn, *m = vm;                             \
1732                                                                 \
1733     for (i = 0; i < opr_sz / sizeof(TYPE); ++i) {               \
1734         d[i] = n[i] < m[i] ? m[i] - n[i] : n[i] - m[i];         \
1735     }                                                           \
1736     clear_tail(d, opr_sz, simd_maxsz(desc));                    \
1737 }
1738 
1739 DO_ABD(gvec_sabd_b, int8_t)
1740 DO_ABD(gvec_sabd_h, int16_t)
1741 DO_ABD(gvec_sabd_s, int32_t)
1742 DO_ABD(gvec_sabd_d, int64_t)
1743 
1744 DO_ABD(gvec_uabd_b, uint8_t)
1745 DO_ABD(gvec_uabd_h, uint16_t)
1746 DO_ABD(gvec_uabd_s, uint32_t)
1747 DO_ABD(gvec_uabd_d, uint64_t)
1748 
1749 #undef DO_ABD
1750 
1751 #define DO_ABA(NAME, TYPE)                                      \
1752 void HELPER(NAME)(void *vd, void *vn, void *vm, uint32_t desc)  \
1753 {                                                               \
1754     intptr_t i, opr_sz = simd_oprsz(desc);                      \
1755     TYPE *d = vd, *n = vn, *m = vm;                             \
1756                                                                 \
1757     for (i = 0; i < opr_sz / sizeof(TYPE); ++i) {               \
1758         d[i] += n[i] < m[i] ? m[i] - n[i] : n[i] - m[i];        \
1759     }                                                           \
1760     clear_tail(d, opr_sz, simd_maxsz(desc));                    \
1761 }
1762 
1763 DO_ABA(gvec_saba_b, int8_t)
1764 DO_ABA(gvec_saba_h, int16_t)
1765 DO_ABA(gvec_saba_s, int32_t)
1766 DO_ABA(gvec_saba_d, int64_t)
1767 
1768 DO_ABA(gvec_uaba_b, uint8_t)
1769 DO_ABA(gvec_uaba_h, uint16_t)
1770 DO_ABA(gvec_uaba_s, uint32_t)
1771 DO_ABA(gvec_uaba_d, uint64_t)
1772 
1773 #undef DO_ABA
1774 
1775 #define DO_NEON_PAIRWISE(NAME, OP)                                      \
1776     void HELPER(NAME##s)(void *vd, void *vn, void *vm,                  \
1777                          void *stat, uint32_t oprsz)                    \
1778     {                                                                   \
1779         float_status *fpst = stat;                                      \
1780         float32 *d = vd;                                                \
1781         float32 *n = vn;                                                \
1782         float32 *m = vm;                                                \
1783         float32 r0, r1;                                                 \
1784                                                                         \
1785         /* Read all inputs before writing outputs in case vm == vd */   \
1786         r0 = float32_##OP(n[H4(0)], n[H4(1)], fpst);                    \
1787         r1 = float32_##OP(m[H4(0)], m[H4(1)], fpst);                    \
1788                                                                         \
1789         d[H4(0)] = r0;                                                  \
1790         d[H4(1)] = r1;                                                  \
1791     }                                                                   \
1792                                                                         \
1793     void HELPER(NAME##h)(void *vd, void *vn, void *vm,                  \
1794                          void *stat, uint32_t oprsz)                    \
1795     {                                                                   \
1796         float_status *fpst = stat;                                      \
1797         float16 *d = vd;                                                \
1798         float16 *n = vn;                                                \
1799         float16 *m = vm;                                                \
1800         float16 r0, r1, r2, r3;                                         \
1801                                                                         \
1802         /* Read all inputs before writing outputs in case vm == vd */   \
1803         r0 = float16_##OP(n[H2(0)], n[H2(1)], fpst);                    \
1804         r1 = float16_##OP(n[H2(2)], n[H2(3)], fpst);                    \
1805         r2 = float16_##OP(m[H2(0)], m[H2(1)], fpst);                    \
1806         r3 = float16_##OP(m[H2(2)], m[H2(3)], fpst);                    \
1807                                                                         \
1808         d[H4(0)] = r0;                                                  \
1809         d[H4(1)] = r1;                                                  \
1810         d[H4(2)] = r2;                                                  \
1811         d[H4(3)] = r3;                                                  \
1812     }
1813 
1814 DO_NEON_PAIRWISE(neon_padd, add)
1815 DO_NEON_PAIRWISE(neon_pmax, max)
1816 DO_NEON_PAIRWISE(neon_pmin, min)
1817 
1818 #undef DO_NEON_PAIRWISE
1819