aboutsummaryrefslogtreecommitdiff
path: root/vendor/gmp-6.3.0/mpn/s390_64/z13/addmul_1.c
blob: 022e5edcc3150a29cd0748ef80072ffd14107a51 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
/* Addmul_1 / mul_1 for IBM z13 and later
   Contributed by Marius Hillenbrand

Copyright 2021 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of either:

  * the GNU Lesser General Public License as published by the Free
    Software Foundation; either version 3 of the License, or (at your
    option) any later version.

or

  * the GNU General Public License as published by the Free Software
    Foundation; either version 2 of the License, or (at your option) any
    later version.

or both in parallel, as here.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received copies of the GNU General Public License and the
GNU Lesser General Public License along with the GNU MP Library.  If not,
see https://www.gnu.org/licenses/.  */

#include "gmp-impl.h"
#include "s390_64/z13/common-vec.h"

#undef FUNCNAME

#ifdef DO_INLINE
#  ifdef OPERATION_addmul_1
#    define ADD
#    define FUNCNAME inline_addmul_1
#  elif defined(OPERATION_mul_1)
#    define FUNCNAME inline_mul_1
#  endif

#else
#  ifdef OPERATION_addmul_1
#    define ADD
#    define FUNCNAME mpn_addmul_1
#  elif defined(OPERATION_mul_1)
#    define FUNCNAME mpn_mul_1
#  endif
#endif

#ifdef DO_INLINE
static inline mp_limb_t
FUNCNAME (mp_ptr rp, mp_srcptr s1p, mp_size_t n, mp_limb_t s2limb)
    __attribute__ ((always_inline));

static inline
#endif
mp_limb_t
FUNCNAME (mp_ptr rp, mp_srcptr s1p, mp_size_t n, mp_limb_t s2limb)
{
  ASSERT (n >= 1);
  ASSERT (MPN_SAME_OR_INCR_P(rp, s1p, n));

  /* Combine 64x64 multiplication into GPR pairs (MLGR) with 128-bit adds in
     VRs (using each VR as a single 128-bit accumulator).
     The inner loop is unrolled to four limbs, with two blocks of four
     multiplications each. Since the MLGR operation operates on even/odd GPR
     pairs, pin the products appropriately. */

  /* products as GPR pairs */
  register mp_limb_t p0_high asm("r0");
  register mp_limb_t p0_low asm("r1");

  register mp_limb_t p1_high asm("r8");
  register mp_limb_t p1_low asm("r9");

  register mp_limb_t p2_high asm("r6");
  register mp_limb_t p2_low asm("r7");

  register mp_limb_t p3_high asm("r10");
  register mp_limb_t p3_low asm("r11");

  /* carry flag for 128-bit add in VR for first carry chain */
  vec_t carry_vec0 = { .dw = vec_splat_u64 (0) };
  mp_limb_t carry_limb = 0;

#ifdef ADD
  /* 2nd carry flag for 2nd carry chain with addmul */
  vec_t carry_vec1 = { .dw = vec_splat_u64 (0) };
  vec_t sum0;
  vec_t rp0_addend, rp1_addend;
  rp0_addend.dw = vec_splat_u64 (0);
  rp1_addend.dw = vec_splat_u64 (0);
#endif
  vec_t sum1;

  vec_t carry_prod = { .dw = vec_splat_u64 (0) };

  /* The scalar multiplications compete with pointer and index increments for
   * issue ports. Thus, increment the loop index in the middle of the loop so
   * that the operations for the next iteration's multiplications can be
   * loaded in time (looks horrible, yet helps performance) and make sure we
   * use addressing with base reg + index reg + immediate displacement
   * (so that only the single index needs incrementing, instead of multiple
   * pointers). */
#undef LOOP_ADVANCE
#undef IDX_OFFSET

#define LOOP_ADVANCE 4 * sizeof (mp_limb_t)
#define IDX_OFFSET (LOOP_ADVANCE)
  register ssize_t idx = 0 - IDX_OFFSET;

  /*
   * branch-on-count implicitly hint to the branch prediction as taken, while
   * compare-and-branch hints as not taken. currently, using branch-on-count
   * has a performance advantage, but it is not clear that it is generally the
   * better choice (e.g., branch-on-count requires decrementing the separate
   * counter). so, allow switching the loop condition to enable either
   * category of branch instructions:
   * - idx is less than an upper bound, for compare-and-branch
   * - iteration counter greater than zero, for branch-on-count
   */
#define BRCTG
#ifdef BRCTG
  ssize_t iterations = (size_t)n / 4;
#else
  ssize_t const idx_bound = n * sizeof (mp_limb_t) - IDX_OFFSET;
#endif

  /* products will be transferred into VRs before adding up.
   * see main loop below for comments on accumulation scheme. */
  vec_t product0, product1, product2;

  product0.dw = vec_splat_u64 (0);

  switch ((size_t)n % 4)
    {
    case 0:
      break;

    case 1:
      idx = 1 * sizeof (mp_limb_t) - IDX_OFFSET;

      p3_low = s1p[0];
      s390_umul_ppmm (p3_high, p3_low, s2limb);

#ifdef ADD
      rp0_addend.dw[1] = rp[0];
      product0.dw[1] = p3_low;

      sum0.sw = vec_add_u128 (product0.sw, rp0_addend.sw);
      carry_vec1.dw = vec_permi (sum0.dw, sum0.dw, 0);

      rp[0] = sum0.dw[1];
#else
      rp[0] = p3_low;
#endif

      carry_limb = p3_high;
      break;

    case 2:
      p0_low = s1p[0];
      p3_low = s1p[1];
      idx = 2 * sizeof (mp_limb_t) - IDX_OFFSET;

      s390_double_umul_ppmm (p0_high, p0_low, p3_high, p3_low, s2limb);

      carry_prod.dw[0] = p3_low;

      product0.dw = vec_load_2di_as_pair (p0_high, p0_low);

      carry_limb = p3_high;

#ifdef ADD
      rp0_addend = vec_load_elements_reversed (rp, 0);
      sum0.sw = vec_add_u128 (carry_prod.sw, rp0_addend.sw);
      carry_vec0.sw = vec_addc_u128 (carry_prod.sw, rp0_addend.sw);

      sum1.sw = vec_add_u128 (sum0.sw, product0.sw);
      carry_vec1.sw = vec_addc_u128 (sum0.sw, product0.sw);
#else
      sum1.sw = vec_add_u128 (carry_prod.sw, product0.sw);
      carry_vec0.sw = vec_addc_u128 (carry_prod.sw, product0.sw);
#endif

      vec_store_elements_reversed (rp, 0, sum1);

      break;

    case 3:
      idx = 3 * sizeof (mp_limb_t) - IDX_OFFSET;

      p0_low = s1p[0];
      s390_umul_ppmm (p0_high, p0_low, s2limb);

#ifdef ADD
      rp0_addend.dw[1] = rp[0];
      product0.dw[1] = p0_low;

      sum0.sw = vec_add_u128 (product0.sw, rp0_addend.sw);
      carry_vec1.dw = vec_permi (sum0.dw, sum0.dw, 0);

      rp[0] = sum0.dw[1];
#else
      rp[0] = p0_low;
#endif
      carry_limb = p0_high;

      p1_low = s1p[1];
      p3_low = s1p[2];

      s390_double_umul_ppmm (p1_high, p1_low, p3_high, p3_low, s2limb);

      carry_prod.dw = vec_load_2di_as_pair (p3_low, carry_limb);
      product1.dw = vec_load_2di_as_pair (p1_high, p1_low);
      carry_limb = p3_high;

#ifdef ADD
      rp0_addend = vec_load_elements_reversed (rp, 8);
      sum0.sw = vec_add_u128 (carry_prod.sw, rp0_addend.sw);
      carry_vec0.sw = vec_addc_u128 (carry_prod.sw, rp0_addend.sw);

      sum1.sw = vec_adde_u128 (sum0.sw, product1.sw, carry_vec1.sw);
      carry_vec1.sw = vec_addec_u128 (sum0.sw, product1.sw, carry_vec1.sw);
#else
      sum1.sw = vec_adde_u128 (carry_prod.sw, product1.sw, carry_vec0.sw);
      carry_vec0.sw
          = vec_addec_u128 (carry_prod.sw, product1.sw, carry_vec0.sw);
#endif
      vec_store_elements_reversed (rp, 8, sum1);
      break;
    }

#ifdef BRCTG
  for (; iterations > 0; iterations--)
    {
#else
  while (idx < idx_bound)
    {
#endif
      vec_t overlap_addend0;
      vec_t overlap_addend1;

      /* The 64x64->128 MLGR multiplies two factors in GPRs and stores the
       * result in a GPR pair. One of the factors is taken from the GPR pair
       * and overwritten.
       * To reuse factors, it turned out cheaper to load limbs multiple times
       * than copying GPR contents. Enforce that and the use of addressing by
       * base + index gpr + immediate displacement via inline asm.
       */
      ASM_LOADGPR (p0_low, s1p, idx, 0 + IDX_OFFSET);
      ASM_LOADGPR (p1_low, s1p, idx, 8 + IDX_OFFSET);
      ASM_LOADGPR (p2_low, s1p, idx, 16 + IDX_OFFSET);
      ASM_LOADGPR (p3_low, s1p, idx, 24 + IDX_OFFSET);

      /*
       * accumulate products as follows (for addmul):
       *                       | rp[i+3] | rp[i+2] | rp[i+1] | rp[i]   |
       *                                             p0_high | p0_low  |
       *                                   p1_high | p1_low  | carry-limb in
       *                         p2_high | p2_low  |
       * c-limb out <- p3_high | p3_low  |
       *                       | <    128-bit VR   > <   128-bit VR    >
       *
       *                         <   rp1_addend    > <  rp0_addend     >
       *     carry-chain 0  <-   +           <-      +  <- carry_vec0[127]
       *                         <   product1      > <  product0       >
       *     carry-chain 1  <-   +           <-      +  <- carry_vec1[127]
       *                         < overlap_addend1 > < overlap_addend0 >
       *
       * note that a 128-bit add with carry in + out is built from two insns
       * - vec_adde_u128 (vacq) provides sum
       * - vec_addec_u128 (vacccq) provides the new carry bit
       */

      s390_double_umul_ppmm (p0_high, p0_low, p1_high, p1_low, s2limb);

      /*
       * "barrier" to enforce scheduling loads for all limbs and first round
       * of MLGR before anything else.
       */
      asm volatile("");

      product0.dw = vec_load_2di_as_pair (p0_high, p0_low);

#ifdef ADD
      rp0_addend = vec_load_elements_reversed_idx (rp, idx, 0 + IDX_OFFSET);
      rp1_addend = vec_load_elements_reversed_idx (rp, idx, 16 + IDX_OFFSET);
#endif
      /* increment loop index to unblock dependant loads of limbs for the next
       * iteration (see above at #define LOOP_ADVANCE) */
      idx += LOOP_ADVANCE;

      s390_double_umul_ppmm (p2_high, p2_low, p3_high, p3_low, s2limb);

      overlap_addend0.dw = vec_load_2di_as_pair (p1_low, carry_limb);
      asm volatile("");

#ifdef ADD
      sum0.sw = vec_adde_u128 (product0.sw, rp0_addend.sw, carry_vec0.sw);
      sum1.sw = vec_adde_u128 (sum0.sw, overlap_addend0.sw, carry_vec1.sw);

      carry_vec0.sw
          = vec_addec_u128 (product0.sw, rp0_addend.sw, carry_vec0.sw);
      carry_vec1.sw
          = vec_addec_u128 (sum0.sw, overlap_addend0.sw, carry_vec1.sw);
#else
      sum1.sw = vec_adde_u128 (product0.sw, overlap_addend0.sw, carry_vec0.sw);
      carry_vec0.sw
          = vec_addec_u128 (product0.sw, overlap_addend0.sw, carry_vec0.sw);
#endif

      asm volatile("");
      product2.dw = vec_load_2di_as_pair (p2_high, p2_low);
      overlap_addend1.dw = vec_load_2di_as_pair (p3_low, p1_high);

      vec_t sum4;

#ifdef ADD
      vec_t sum3;
      sum3.sw = vec_adde_u128 (product2.sw, rp1_addend.sw, carry_vec0.sw);
      sum4.sw = vec_adde_u128 (sum3.sw, overlap_addend1.sw, carry_vec1.sw);

      carry_vec0.sw
          = vec_addec_u128 (product2.sw, rp1_addend.sw, carry_vec0.sw);
      carry_vec1.sw
          = vec_addec_u128 (sum3.sw, overlap_addend1.sw, carry_vec1.sw);
#else
      sum4.sw = vec_adde_u128 (product2.sw, overlap_addend1.sw, carry_vec0.sw);
      carry_vec0.sw
          = vec_addec_u128 (product2.sw, overlap_addend1.sw, carry_vec0.sw);
#endif
      vec_store_elements_reversed_idx (rp, idx, IDX_OFFSET - LOOP_ADVANCE,
                                       sum1);
      vec_store_elements_reversed_idx (rp, idx, 16 + IDX_OFFSET - LOOP_ADVANCE,
                                       sum4);

      carry_limb = p3_high;
    }

#ifdef ADD
  carry_vec0.dw += carry_vec1.dw;
  carry_limb += carry_vec0.dw[1];
#else
  carry_limb += carry_vec0.dw[1];
#endif

  return carry_limb;
}

#undef OPERATION_addmul_1
#undef OPERATION_mul_1
#undef FUNCNAME
#undef ADD