• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 // Auto-generated file. Do not edit!
2 //   Template: src/f32-raddstoreexpminusmax/neon-lut64-p2.c.in
3 //   Generator: tools/xngen
4 //
5 // Copyright 2020 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 <arm_neon.h>
13 
14 #include <xnnpack/common.h>
15 #include <xnnpack/raddstoreexpminusmax.h>
16 
17 
18 extern XNN_INTERNAL const float xnn_table_exp2_k_over_64[64];
19 
xnn_f32_raddstoreexpminusmax_ukernel__neon_lut64_p2_x8(size_t elements,const float * input,float * output,float * sum,float max)20 void xnn_f32_raddstoreexpminusmax_ukernel__neon_lut64_p2_x8(
21     size_t elements,
22     const float* input,
23     float* output,
24     float* sum,
25     float max)
26 {
27   assert(elements % sizeof(float) == 0);
28 
29   const float32x4_t vmagic_bias = vmovq_n_f32(0x1.800000p23f);
30   // The smallest x for which expf(x) is normalized.
31   const float32x4_t vdenorm_cutoff = vmovq_n_f32(-0x1.5D589Ep6f);
32   const float32x4_t vlog2e_x64  = vmovq_n_f32(0x1.715476p6f);
33   // Last 13 bits are zeroes
34   const float32x4_t vminus_ln2_o64_hi = vmovq_n_f32(-0x1.630000p-7f);
35   const float32x4_t vminus_ln2_o64_lo = vmovq_n_f32(0x1.BD0106p-19f);
36 
37   const float32x4_t vc2 = vmovq_n_f32(0x1.FFFF0Ap-2f);
38 
39   const int32x4_t vindex_mask = vmovq_n_s32(INT32_C(0x3F));
40 
41   const float32x4_t vi_max = vdupq_n_f32(max);
42 
43   float32x4_t vacc0 = vmovq_n_f32(0.0f);
44   for (; elements >= 8 * sizeof(float); elements -= 8 * sizeof(float)) {
45     // Load 8 (2x4) inputs at a time.
46     const float32x4_t vi0123 = vld1q_f32(input); input += 4;
47     const float32x4_t vi4567 = vld1q_f32(input); input += 4;
48 
49     // Subtract maximum input x := i - i_max. This implies x <= 0.
50     const float32x4_t vx0123 = vsubq_f32(vi0123, vi_max);
51     const float32x4_t vx4567 = vsubq_f32(vi4567, vi_max);
52 
53     // Compute reduced argument n := round(x * 64 / log(2)).
54     // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing
55     // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction.
56     // The trick with adding large number is valid only within certain bounds (|x * 64 / log(2)| <= 2**22, i.e.
57     // |x| <= 0x1.62E43p+15 = 45426.09375), but that is acceptable, because inputs outside of [-87.336540, 0.0]
58     // result in denormalized or underflown expf(x). We fixup the result for such inputs at the very end of the
59     // algorithm.
60     float32x4_t vn0123 = vmlaq_f32(vmagic_bias, vx0123, vlog2e_x64);
61     float32x4_t vn4567 = vmlaq_f32(vmagic_bias, vx4567, vlog2e_x64);
62 
63     // Create a floating-point number s (scale) such that s := 2**(n / 64) for such inputs that expf(x) is normalized,
64     // i.e. -87.33642 <= x <= 0.0. As n has 6 fractional bits, we split s == 2**(n / 64) = 2**e * 2**(n / 64 - e), where
65     // e := int(n / 64). We create s in two steps:
66     // 1. Fetch 2**(n / 64 - e) = 2**(n % 64) from the table using the 6 low bits of n, as integer. Note that the
67     //    fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0.
68     // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized
69     //    number, because for -87.33642 <= x <= 0.0 (inputs for which expf(x) is normalized) we have -126 <= e <= 0,
70     //    and thus the adjusted exponent is not lower than -126.
71     //
72     // Extract e from bits 6:14 of n and shift it into bits 23:31 (position of floating-point exponent).
73     const int32x4_t ve0123 = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn0123), vmovq_n_s32(INT32_C(0x3F))), 17);
74     const int32x4_t ve4567 = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn4567), vmovq_n_s32(INT32_C(0x3F))), 17);
75 
76     // Use bits 0:6 bits of n, as integer, as an index for table lookup of l := 2**(n % 64).
77     const uint64x2_t vidx0123 = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn0123), vindex_mask));
78     const uint64_t vidx01 = vgetq_lane_u64(vidx0123, 0);
79     const uint64_t vidx23 = vgetq_lane_u64(vidx0123, 1);
80     const uint64x2_t vidx4567 = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn4567), vindex_mask));
81     const uint64_t vidx45 = vgetq_lane_u64(vidx4567, 0);
82     const uint64_t vidx67 = vgetq_lane_u64(vidx4567, 1);
83 
84     float32x2_t vl01 = vld1_dup_f32(&xnn_table_exp2_k_over_64[(uint32_t) vidx01]);
85     float32x2_t vl23 = vld1_dup_f32(&xnn_table_exp2_k_over_64[(uint32_t) vidx23]);
86     float32x2_t vl45 = vld1_dup_f32(&xnn_table_exp2_k_over_64[(uint32_t) vidx45]);
87     float32x2_t vl67 = vld1_dup_f32(&xnn_table_exp2_k_over_64[(uint32_t) vidx67]);
88 
89     vl01 = vld1_lane_f32(&xnn_table_exp2_k_over_64[(uint32_t) (vidx01 >> 32)], vl01, 1);
90     vl23 = vld1_lane_f32(&xnn_table_exp2_k_over_64[(uint32_t) (vidx23 >> 32)], vl23, 1);
91     const float32x4_t vl0123 = vcombine_f32(vl01, vl23);
92     vl45 = vld1_lane_f32(&xnn_table_exp2_k_over_64[(uint32_t) (vidx45 >> 32)], vl45, 1);
93     vl67 = vld1_lane_f32(&xnn_table_exp2_k_over_64[(uint32_t) (vidx67 >> 32)], vl67, 1);
94     const float32x4_t vl4567 = vcombine_f32(vl45, vl67);
95 
96     // Adjust exponent of the value l fetched from the table to get the final s value.
97     const float32x4_t vs0123 = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl0123), ve0123));
98     const float32x4_t vs4567 = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl4567), ve4567));
99 
100     // Subtract the large number back to get final n := round(x * 64 / log(2)) as a floating-point number.
101     vn0123 = vsubq_f32(vn0123, vmagic_bias);
102     vn4567 = vsubq_f32(vn4567, vmagic_bias);
103 
104     // Compute reduced argument t := x - n * log(2) / 64.
105     // Use Cody-Waite range reduction method (note the two constants representing log(2) / 64) to improve accuracy.
106     float32x4_t vt0123 = vmlaq_f32(vx0123, vn0123, vminus_ln2_o64_hi);
107     float32x4_t vt4567 = vmlaq_f32(vx4567, vn4567, vminus_ln2_o64_hi);
108 
109     vt0123 = vmlaq_f32(vt0123, vn0123, vminus_ln2_o64_lo);
110     vt4567 = vmlaq_f32(vt4567, vn4567, vminus_ln2_o64_lo);
111 
112     // Compute degree-2 polynomial approxiatmion for exp(t) on [-log(2)/128, log(2)/128].
113     float32x4_t vp0123 = vmulq_f32(vt0123, vc2);
114     float32x4_t vp4567 = vmulq_f32(vt4567, vc2);
115 
116     vp0123 = vmlaq_f32(vt0123, vt0123, vp0123);
117     vp4567 = vmlaq_f32(vt4567, vt4567, vp4567);
118 
119     // Reconstruct the final f value:
120     //   f = s * (1 + t * (1 + t * c2))
121     //     = s * (1 + t + t * (t * c2))
122     //     = s + s * (t + t * (t * c2))
123     //     = s + s * p
124     float32x4_t vf0123 = vmlaq_f32(vs0123, vs0123, vp0123);
125     float32x4_t vf4567 = vmlaq_f32(vs4567, vs4567, vp4567);
126 
127     // For inputs below denormal cutoff, replace output with +0.0f.
128     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
129     vf0123 = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf0123), vcltq_f32(vx0123, vdenorm_cutoff)));
130     vf4567 = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf4567), vcltq_f32(vx4567, vdenorm_cutoff)));
131 
132     // Store 8 (2x4) outputs at a time.
133     vst1q_f32(output, vf0123); output += 4;
134     vst1q_f32(output, vf4567); output += 4;
135 
136     // Accumulate computed exponents.
137     vacc0 = vaddq_f32(vacc0, vf0123);
138     vacc0 = vaddq_f32(vacc0, vf4567);
139   }
140 
141   float32x4_t vacc = vacc0;
142   for (; elements >= 4 * sizeof(float); elements -= 4 * sizeof(float)) {
143     // Load 4 inputs at a time.
144     const float32x4_t vi = vld1q_f32(input); input += 4;
145 
146     // Subtract maximum input x := i - i_max. This implies x <= 0.
147     const float32x4_t vx = vsubq_f32(vi, vi_max);
148 
149     // Compute reduced argument n := round(x * 64 / log(2)).
150     // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing
151     // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction.
152     // The trick with adding large number is valid only within certain bounds (|x * 64 / log(2)| <= 2**22, i.e.
153     // |x| <= 0x1.62E43p+15 = 45426.09375), but that is acceptable, because inputs outside of [-87.336540, 0.0]
154     // result in denormalized or underflown expf(x). We fixup the result for such inputs at the very end of the
155     // algorithm.
156     float32x4_t vn = vmlaq_f32(vmagic_bias, vx, vlog2e_x64);
157 
158     // Create a floating-point number s (scale) such that s := 2**(n / 64) for such inputs that expf(x) is normalized,
159     // i.e. -87.33642 <= x <= 0.0. As n has 6 fractional bits, we split s == 2**(n / 64) = 2**e * 2**(n / 64 - e), where
160     // e := int(n / 64). We create s in two steps:
161     // 1. Fetch 2**(n / 64 - e) = 2**(n % 64) from the table using the 6 low bits of n, as integer. Note that the
162     //    fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0.
163     // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized
164     //    number, because for -87.33642 <= x <= 0.0 (inputs for which expf(x) is normalized) we have -126 <= e <= 0,
165     //    and thus the adjusted exponent is not lower than -126.
166     //
167     // Extract e from bits 6:14 of n and shift it into bits 23:31 (position of floating-point exponent).
168     const int32x4_t ve = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn), vmovq_n_s32(INT32_C(0x3F))), 17);
169 
170     // Use bits 0:6 bits of n, as integer, as an index for table lookup of l := 2**(n % 64).
171     const uint64x2_t vidx = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn), vindex_mask));
172     const uint64_t vidx_lo = vgetq_lane_u64(vidx, 0);
173     const uint64_t vidx_hi = vgetq_lane_u64(vidx, 1);
174     float32x2_t vl_lo = vld1_dup_f32(&xnn_table_exp2_k_over_64[(uint32_t) vidx_lo]);
175     float32x2_t vl_hi = vld1_dup_f32(&xnn_table_exp2_k_over_64[(uint32_t) vidx_hi]);
176     vl_lo = vld1_lane_f32(&xnn_table_exp2_k_over_64[(uint32_t) (vidx_lo >> 32)], vl_lo, 1);
177     vl_hi = vld1_lane_f32(&xnn_table_exp2_k_over_64[(uint32_t) (vidx_hi >> 32)], vl_hi, 1);
178     const float32x4_t vl = vcombine_f32(vl_lo, vl_hi);
179     // Adjust exponent of the value l fetched from the table to get the final s value.
180     const float32x4_t vs = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl), ve));
181 
182     // Subtract the large number back to get final n := round(x * 64 / log(2)) as a floating-point number.
183     vn = vsubq_f32(vn, vmagic_bias);
184 
185     // Compute reduced argument t := x - n * log(2) / 64.
186     // Use Cody-Waite range reduction method (note the two constants representing log(2) / 64) to improve accuracy.
187     float32x4_t vt = vmlaq_f32(vx, vn, vminus_ln2_o64_hi);
188     vt = vmlaq_f32(vt, vn, vminus_ln2_o64_lo);
189 
190     // Compute degree-2 polynomial approxiatmion for exp(t) on [-log(2)/128, log(2)/128].
191     float32x4_t vp = vmulq_f32(vt, vc2);
192     vp = vmlaq_f32(vt, vt, vp);
193 
194     // Reconstruct the final f value:
195     //   f = s * (1 + t * (1 + t * c2))
196     //     = s * (1 + t + t * (t * c2))
197     //     = s + s * (t + t * (t * c2))
198     //     = s + s * p
199     float32x4_t vf = vmlaq_f32(vs, vs, vp);
200 
201     // For inputs below denormal cutoff, replace output with +0.0f.
202     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
203     vf = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf), vcltq_f32(vx, vdenorm_cutoff)));
204 
205     // Store 4 outputs at a time.
206     vst1q_f32(output, vf); output += 4;
207 
208     // Accumulate computed exponents.
209     vacc = vaddq_f32(vacc, vf);
210   }
211 #if XNN_ARCH_ARM64
212   float vacc_lo = vaddvq_f32(vacc);
213 #else
214   float32x2_t vacc_lo = vadd_f32(vget_high_f32(vacc), vget_low_f32(vacc));
215 #endif
216   if (elements != 0) {
217     assert(elements >= 1 * sizeof(float));
218     assert(elements <= 3 * sizeof(float));
219     // Load 4 inputs at a time.
220     const float32x4_t vi = vld1q_f32(input); input += 4;
221 
222     // Subtract maximum input x := i - i_max. This implies x <= 0.
223     const float32x4_t vx = vsubq_f32(vi, vi_max);
224 
225     // Compute reduced argument n := round(x * 64 / log(2)).
226     // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing
227     // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction.
228     // The trick with adding large number is valid only within certain bounds (|x * 64 / log(2)| <= 2**22, i.e.
229     // |x| <= 0x1.62E43p+15 = 45426.09375), but that is acceptable, because inputs outside of [-87.336540, 0.0]
230     // result in denormalized or underflown expf(x). We fixup the result for such inputs at the very end of the
231     // algorithm.
232     float32x4_t vn = vmlaq_f32(vmagic_bias, vx, vlog2e_x64);
233 
234     // Create a floating-point number s (scale) such that s := 2**(n / 64) for such inputs that expf(x) is normalized,
235     // i.e. -87.33642 <= x <= 0.0. As n has 6 fractional bits, we split s == 2**(n / 64) = 2**e * 2**(n / 64 - e), where
236     // e := int(n / 64). We create s in two steps:
237     // 1. Fetch 2**(n / 64 - e) = 2**(n % 64) from the table using the 6 low bits of n, as integer. Note that the
238     //    fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0.
239     // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized
240     //    number, because for -87.33642 <= x <= 0.0 (inputs for which expf(x) is normalized) we have -126 <= e <= 0,
241     //    and thus the adjusted exponent is not lower than -126.
242     //
243     // Extract e from bits 6:14 of n and shift it into bits 23:31 (position of floating-point exponent).
244     const int32x4_t ve = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn), vmovq_n_s32(INT32_C(0x3F))), 17);
245 
246     // Use bits 0:6 bits of n, as integer, as an index for table lookup of l := 2**(n % 64).
247     const uint64x2_t vidx = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn), vindex_mask));
248     const uint64_t vidx_lo = vgetq_lane_u64(vidx, 0);
249     const uint64_t vidx_hi = vgetq_lane_u64(vidx, 1);
250     float32x2_t vl_lo = vld1_dup_f32(&xnn_table_exp2_k_over_64[(uint32_t) vidx_lo]);
251     float32x2_t vl_hi = vld1_dup_f32(&xnn_table_exp2_k_over_64[(uint32_t) vidx_hi]);
252     vl_lo = vld1_lane_f32(&xnn_table_exp2_k_over_64[(uint32_t) (vidx_lo >> 32)], vl_lo, 1);
253     vl_hi = vld1_lane_f32(&xnn_table_exp2_k_over_64[(uint32_t) (vidx_hi >> 32)], vl_hi, 1);
254     const float32x4_t vl = vcombine_f32(vl_lo, vl_hi);
255     // Adjust exponent of the value l fetched from the table to get the final s value.
256     const float32x4_t vs = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl), ve));
257 
258     // Subtract the large number back to get final n := round(x * 64 / log(2)) as a floating-point number.
259     vn = vsubq_f32(vn, vmagic_bias);
260 
261     // Compute reduced argument t := x - n * log(2) / 64.
262     // Use Cody-Waite range reduction method (note the two constants representing log(2) / 64) to improve accuracy.
263     float32x4_t vt = vmlaq_f32(vx, vn, vminus_ln2_o64_hi);
264     vt = vmlaq_f32(vt, vn, vminus_ln2_o64_lo);
265 
266     // Compute degree-2 polynomial approxiatmion for exp(t) on [-log(2)/128, log(2)/128].
267     float32x4_t vp = vmulq_f32(vt, vc2);
268     vp = vmlaq_f32(vt, vt, vp);
269 
270     // Reconstruct the final f value:
271     //   f = s * (1 + t * (1 + t * c2))
272     //     = s * (1 + t + t * (t * c2))
273     //     = s + s * (t + t * (t * c2))
274     //     = s + s * p
275     float32x4_t vf = vmlaq_f32(vs, vs, vp);
276 
277     // For inputs below denormal cutoff, replace output with +0.0f.
278     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
279     vf = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf), vcltq_f32(vx, vdenorm_cutoff)));
280 
281     float32x2_t vf_lo = vget_low_f32(vf);
282     if (elements & (2 * sizeof(float))) {
283       // Store 2 outputs at a time.
284       vst1_f32(output, vf_lo); output += 2;
285 
286       // Accumulate 2 computed exponents.
287       #if XNN_ARCH_ARM64
288         vacc_lo += vaddv_f32(vf_lo);
289       #else
290         vacc_lo = vadd_f32(vacc_lo, vf_lo);
291       #endif
292 
293       vf_lo = vget_high_f32(vf);
294     }
295     if (elements & (1 * sizeof(float))) {
296       // Store 1 output at a time.
297       vst1_lane_f32(output, vf_lo, 0);
298 
299       // Accumulate 1 computed exponent.
300       #if XNN_ARCH_ARM64
301         vacc_lo += vget_lane_f32(vf_lo, 0);
302       #else
303         vacc_lo = vadd_f32(vacc_lo, vreinterpret_f32_u64(vshl_n_u64(vreinterpret_u64_f32(vf_lo), 32)));
304       #endif
305     }
306   }
307   // Reduce 4 elements in the SIMD register
308 #if XNN_ARCH_ARM64
309   *sum = vacc_lo;
310 #else
311   vst1_lane_f32(sum, vpadd_f32(vacc_lo, vacc_lo), 0);
312 #endif
313 }
314