• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 // Auto-generated file. Do not edit!
2 //   Template: src/f32-raddstoreexpminusmax/psimd-p5.c.in
3 //   Generator: tools/xngen
4 //
5 // Copyright 2019 Google LLC
6 //
7 // This source code is licensed under the BSD-style license found in the
8 // LICENSE file in the root directory of this source tree.
9 
10 #include <assert.h>
11 
12 #include <psimd.h>
13 
14 #include <xnnpack/common.h>
15 #include <xnnpack/raddstoreexpminusmax.h>
16 
17 
xnn_f32_raddstoreexpminusmax_ukernel__psimd_p5_x20_acc2(size_t elements,const float * input,float * output,float * sum,float max)18 void xnn_f32_raddstoreexpminusmax_ukernel__psimd_p5_x20_acc2(
19     size_t elements,
20     const float* input,
21     float* output,
22     float* sum,
23     float max)
24 {
25   assert(elements % sizeof(float) == 0);
26 
27   const psimd_f32 vmagic_bias = psimd_splat_f32(0x1.8000FEp23f);
28   // The smallest x for which expf(x) is normalized.
29   const psimd_f32 vdenorm_cutoff = psimd_splat_f32(-0x1.5D589Ep6f);
30   const psimd_f32 vlog2e = psimd_splat_f32(0x1.715476p+0f);
31   // Last 7 bits are zeroes
32   const psimd_f32 vminus_ln2_hi = psimd_splat_f32(-0x1.62E400p-1f);
33   const psimd_f32 vminus_ln2_lo = psimd_splat_f32(-0x1.7F7D1Cp-20f);
34 
35   const psimd_f32 vc1 = psimd_splat_f32(0x1.FFFFF6p-1f);
36   const psimd_f32 vc2 = psimd_splat_f32(0x1.FFFDC6p-2f);
37   const psimd_f32 vc3 = psimd_splat_f32(0x1.555A80p-3f);
38   const psimd_f32 vc4 = psimd_splat_f32(0x1.573A1Ap-5f);
39   const psimd_f32 vc5 = psimd_splat_f32(0x1.0F9F9Cp-7f);
40 
41   const psimd_f32 vi_max = psimd_splat_f32(max);
42 
43   psimd_f32 vacc0 = psimd_zero_f32();
44   psimd_f32 vacc1 = psimd_zero_f32();
45   for (; elements >= 20 * sizeof(float); elements -= 20 * sizeof(float)) {
46     // Load 20 (5x4) inputs at a time.
47     const psimd_f32 vi0123 = psimd_load_f32(input);
48     const psimd_f32 vi4567 = psimd_load_f32(input + 4);
49     const psimd_f32 vi89AB = psimd_load_f32(input + 8);
50     const psimd_f32 viCDEF = psimd_load_f32(input + 12);
51     const psimd_f32 viGHIJ = psimd_load_f32(input + 16);
52     input += 20;
53 
54     // Subtract maximum input x := i - i_max. This implies x <= 0.
55     const psimd_f32 vx0123 = psimd_sub_f32(vi0123, vi_max);
56     const psimd_f32 vx4567 = psimd_sub_f32(vi4567, vi_max);
57     const psimd_f32 vx89AB = psimd_sub_f32(vi89AB, vi_max);
58     const psimd_f32 vxCDEF = psimd_sub_f32(viCDEF, vi_max);
59     const psimd_f32 vxGHIJ = psimd_sub_f32(viGHIJ, vi_max);
60 
61     // Compute reduced argument elements := round(x / log(2)).
62     psimd_f32 vn0123 = psimd_qfma_f32(vmagic_bias, vx0123, vlog2e);
63     psimd_f32 vn4567 = psimd_qfma_f32(vmagic_bias, vx4567, vlog2e);
64     psimd_f32 vn89AB = psimd_qfma_f32(vmagic_bias, vx89AB, vlog2e);
65     psimd_f32 vnCDEF = psimd_qfma_f32(vmagic_bias, vxCDEF, vlog2e);
66     psimd_f32 vnGHIJ = psimd_qfma_f32(vmagic_bias, vxGHIJ, vlog2e);
67 
68     // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e.
69     // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly.
70     const psimd_f32 vs0123 = (psimd_f32) ((psimd_u32) vn0123 << 23);
71     const psimd_f32 vs4567 = (psimd_f32) ((psimd_u32) vn4567 << 23);
72     const psimd_f32 vs89AB = (psimd_f32) ((psimd_u32) vn89AB << 23);
73     const psimd_f32 vsCDEF = (psimd_f32) ((psimd_u32) vnCDEF << 23);
74     const psimd_f32 vsGHIJ = (psimd_f32) ((psimd_u32) vnGHIJ << 23);
75 
76     // Subtract the large number back to get final elements := round(x / log(2)).
77     vn0123 = psimd_sub_f32(vn0123, vmagic_bias);
78     vn4567 = psimd_sub_f32(vn4567, vmagic_bias);
79     vn89AB = psimd_sub_f32(vn89AB, vmagic_bias);
80     vnCDEF = psimd_sub_f32(vnCDEF, vmagic_bias);
81     vnGHIJ = psimd_sub_f32(vnGHIJ, vmagic_bias);
82 
83     // Compute reduced argument t := x - elements * log(2).
84     // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
85     psimd_f32 vt0123 = psimd_qfma_f32(vx0123, vn0123, vminus_ln2_hi);
86     psimd_f32 vt4567 = psimd_qfma_f32(vx4567, vn4567, vminus_ln2_hi);
87     psimd_f32 vt89AB = psimd_qfma_f32(vx89AB, vn89AB, vminus_ln2_hi);
88     psimd_f32 vtCDEF = psimd_qfma_f32(vxCDEF, vnCDEF, vminus_ln2_hi);
89     psimd_f32 vtGHIJ = psimd_qfma_f32(vxGHIJ, vnGHIJ, vminus_ln2_hi);
90 
91     vt0123 = psimd_qfma_f32(vt0123, vn0123, vminus_ln2_lo);
92     vt4567 = psimd_qfma_f32(vt4567, vn4567, vminus_ln2_lo);
93     vt89AB = psimd_qfma_f32(vt89AB, vn89AB, vminus_ln2_lo);
94     vtCDEF = psimd_qfma_f32(vtCDEF, vnCDEF, vminus_ln2_lo);
95     vtGHIJ = psimd_qfma_f32(vtGHIJ, vnGHIJ, vminus_ln2_lo);
96 
97     // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2].
98     psimd_f32 vp0123 = psimd_qfma_f32(vc4, vc5, vt0123);
99     psimd_f32 vp4567 = psimd_qfma_f32(vc4, vc5, vt4567);
100     psimd_f32 vp89AB = psimd_qfma_f32(vc4, vc5, vt89AB);
101     psimd_f32 vpCDEF = psimd_qfma_f32(vc4, vc5, vtCDEF);
102     psimd_f32 vpGHIJ = psimd_qfma_f32(vc4, vc5, vtGHIJ);
103 
104     vp0123 = psimd_qfma_f32(vc3, vp0123, vt0123);
105     vp4567 = psimd_qfma_f32(vc3, vp4567, vt4567);
106     vp89AB = psimd_qfma_f32(vc3, vp89AB, vt89AB);
107     vpCDEF = psimd_qfma_f32(vc3, vpCDEF, vtCDEF);
108     vpGHIJ = psimd_qfma_f32(vc3, vpGHIJ, vtGHIJ);
109 
110     vp0123 = psimd_qfma_f32(vc2, vp0123, vt0123);
111     vp4567 = psimd_qfma_f32(vc2, vp4567, vt4567);
112     vp89AB = psimd_qfma_f32(vc2, vp89AB, vt89AB);
113     vpCDEF = psimd_qfma_f32(vc2, vpCDEF, vtCDEF);
114     vpGHIJ = psimd_qfma_f32(vc2, vpGHIJ, vtGHIJ);
115 
116     vp0123 = psimd_qfma_f32(vc1, vp0123, vt0123);
117     vp4567 = psimd_qfma_f32(vc1, vp4567, vt4567);
118     vp89AB = psimd_qfma_f32(vc1, vp89AB, vt89AB);
119     vpCDEF = psimd_qfma_f32(vc1, vpCDEF, vtCDEF);
120     vpGHIJ = psimd_qfma_f32(vc1, vpGHIJ, vtGHIJ);
121 
122     // Reconstruct the final f value:
123     //   f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
124     //     = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
125     //     = s + (t * s) * p
126     vt0123 = psimd_mul_f32(vt0123, vs0123);
127     vt4567 = psimd_mul_f32(vt4567, vs4567);
128     vt89AB = psimd_mul_f32(vt89AB, vs89AB);
129     vtCDEF = psimd_mul_f32(vtCDEF, vsCDEF);
130     vtGHIJ = psimd_mul_f32(vtGHIJ, vsGHIJ);
131 
132     psimd_f32 vf0123 = psimd_qfma_f32(vs0123, vt0123, vp0123);
133     psimd_f32 vf4567 = psimd_qfma_f32(vs4567, vt4567, vp4567);
134     psimd_f32 vf89AB = psimd_qfma_f32(vs89AB, vt89AB, vp89AB);
135     psimd_f32 vfCDEF = psimd_qfma_f32(vsCDEF, vtCDEF, vpCDEF);
136     psimd_f32 vfGHIJ = psimd_qfma_f32(vsGHIJ, vtGHIJ, vpGHIJ);
137 
138     // For inputs below zero cutoff, replace output with +0.0f.
139     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
140     vf0123 = psimd_andnotmask_f32(vx0123 < vdenorm_cutoff, vf0123);
141     vf4567 = psimd_andnotmask_f32(vx4567 < vdenorm_cutoff, vf4567);
142     vf89AB = psimd_andnotmask_f32(vx89AB < vdenorm_cutoff, vf89AB);
143     vfCDEF = psimd_andnotmask_f32(vxCDEF < vdenorm_cutoff, vfCDEF);
144     vfGHIJ = psimd_andnotmask_f32(vxGHIJ < vdenorm_cutoff, vfGHIJ);
145 
146     // Store 20 (5x4) outputs at a time.
147     psimd_store_f32(output, vf0123);
148     psimd_store_f32(output + 4, vf4567);
149     psimd_store_f32(output + 8, vf89AB);
150     psimd_store_f32(output + 12, vfCDEF);
151     psimd_store_f32(output + 16, vfGHIJ);
152     output += 20;
153 
154     // Accumulate computed exponents.
155     vacc0 = psimd_add_f32(vacc0, vf0123);
156     vacc0 = psimd_add_f32(vacc0, vf4567);
157     vacc0 = psimd_add_f32(vacc0, vf89AB);
158     vacc0 = psimd_add_f32(vacc0, vfCDEF);
159     vacc0 = psimd_add_f32(vacc0, vfGHIJ);
160   }
161   // Add up all accumulators to vacc0
162   vacc0 = psimd_add_f32(vacc0, vacc1);
163 
164   psimd_f32 vacc = vacc0;
165   for (; elements >= 4 * sizeof(float); elements -= 4 * sizeof(float)) {
166     // Load 4 inputs at a time.
167     const psimd_f32 vi = psimd_load_f32(input);
168     input += 4;
169 
170     // Subtract maximum input x := i - i_max. This implies x <= 0.
171     const psimd_f32 vx = psimd_sub_f32(vi, vi_max);
172 
173     // Compute reduced argument elements := round(x / log(2)).
174     psimd_f32 vn = psimd_qfma_f32(vmagic_bias, vx, vlog2e);
175 
176     // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e.
177     // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly.
178     const psimd_f32 vs = (psimd_f32) ((psimd_u32) vn << 23);
179 
180     // Subtract the large number back to get final elements := round(x / log(2)).
181     vn = psimd_sub_f32(vn, vmagic_bias);
182 
183     // Compute reduced argument t := x - elements * log(2).
184     // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
185     psimd_f32 vt = psimd_qfma_f32(vx, vn, vminus_ln2_hi);
186     vt = psimd_qfma_f32(vt, vn, vminus_ln2_lo);
187 
188     // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2].
189     psimd_f32 vp = psimd_qfma_f32(vc4, vc5, vt);
190     vp = psimd_qfma_f32(vc3, vp, vt);
191     vp = psimd_qfma_f32(vc2, vp, vt);
192     vp = psimd_qfma_f32(vc1, vp, vt);
193 
194     // Reconstruct the final f value:
195     //   f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
196     //     = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
197     //     = s + (t * s) * p
198     vt = psimd_mul_f32(vt, vs);
199     psimd_f32 vf = psimd_qfma_f32(vs, vt, vp);
200 
201     // For inputs below zero cutoff, replace output with +0.0f.
202     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
203     vf = psimd_andnotmask_f32(vx < vdenorm_cutoff, vf);
204 
205     // Store 4 outputs at a time.
206     psimd_store_f32(output, vf);
207     output += 4;
208 
209     // Accumulate computed exponents.
210     vacc = psimd_add_f32(vacc, vf);
211   }
212   if (elements != 0) {
213     assert(elements >= 1 * sizeof(float));
214     assert(elements <= 3 * sizeof(float));
215     // Load 4 inputs at a time.
216     const psimd_f32 vi = psimd_load_f32(input);
217 
218     // Subtract maximum input x := i - i_max. This implies x <= 0.
219     const psimd_f32 vx = psimd_sub_f32(vi, vi_max);
220 
221     // Compute reduced argument elements := round(x / log(2)).
222     psimd_f32 vn = psimd_qfma_f32(vmagic_bias, vx, vlog2e);
223 
224     // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e.
225     // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly.
226     const psimd_f32 vs = (psimd_f32) ((psimd_u32) vn << 23);
227 
228     // Subtract the large number back to get final elements := round(x / log(2)).
229     vn = psimd_sub_f32(vn, vmagic_bias);
230 
231     // Compute reduced argument t := x - elements * log(2).
232     // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
233     psimd_f32 vt = psimd_qfma_f32(vx, vn, vminus_ln2_hi);
234     vt = psimd_qfma_f32(vt, vn, vminus_ln2_lo);
235 
236     // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2].
237     psimd_f32 vp = psimd_qfma_f32(vc4, vc5, vt);
238     vp = psimd_qfma_f32(vc3, vp, vt);
239     vp = psimd_qfma_f32(vc2, vp, vt);
240     vp = psimd_qfma_f32(vc1, vp, vt);
241 
242     // Reconstruct the final f value:
243     //   f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
244     //     = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
245     //     = s + (t * s) * p
246     vt = psimd_mul_f32(vt, vs);
247     psimd_f32 vf = psimd_qfma_f32(vs, vt, vp);
248 
249     // For inputs below zero cutoff, replace output with +0.0f.
250     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
251     vf = psimd_andnotmask_f32(vx < vdenorm_cutoff, vf);
252 
253     if (elements & (2 * sizeof(float))) {
254       // Store 2 outputs at a time.
255       psimd_store2_f32(output, vf);
256       output += 2;
257 
258       // Accumulate 2 computed exponents.
259       vacc = psimd_add_f32(vacc, psimd_concat_lo_f32(vf, psimd_zero_f32()));
260 
261       vf = psimd_concat_hi_f32(vf, vf);
262     }
263     if (elements & (1 * sizeof(float))) {
264       // Store 1 output at a time.
265       psimd_store1_f32(output, vf);
266 
267       // Accumulate 1 computed exponent.
268       const psimd_f32 vzero = psimd_zero_f32();
269       vf = psimd_concat_lo_f32(vf, vzero);
270       vf = psimd_concat_even_f32(vf, vzero);
271       vacc = psimd_add_f32(vacc, vf);
272     }
273   }
274   // Reduce 4 elements in the SIMD register
275   *sum = psimd_reduce_sum_f32(vacc);
276 }
277