• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 // Auto-generated file. Do not edit!
2 //   Template: src/f32-sigmoid/neon-lut2048-p1.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 <arm_neon.h>
13 
14 #include <xnnpack/common.h>
15 #include <xnnpack/vunary.h>
16 
17 
18 extern XNN_INTERNAL const float xnn_table_exp2_k_over_2048[2048];
19 
xnn_f32_sigmoid_ukernel__neonfma_rr1_lut2048_p1_div_x8(size_t n,const float * x,float * y,const void * params)20 void xnn_f32_sigmoid_ukernel__neonfma_rr1_lut2048_p1_div_x8(
21     size_t n,
22     const float* x,
23     float* y,
24     const void* params)
25 {
26   assert(n % sizeof(float) == 0);
27 
28   const float32x4_t vmagic_bias = vmovq_n_f32(0x1.800000p23f);
29   // The largest z for which sigmoidf(-z) is normalized.
30   // This number is also the largest z for which expf(-z) is normalized.
31   const float32x4_t vdenorm_cutoff = vmovq_n_f32(0x1.5D589Ep+6f);
32   const float32x4_t vminus_log2e_x2048  = vmovq_n_f32(-0x1.715476p11f);
33   const float32x4_t vln2_o2048 = vmovq_n_f32(0x1.62E43p-12f);
34   const float32x4_t vone = vmovq_n_f32(1.0f);
35 
36   const float32x4_t vc1 = vmovq_n_f32(-0x1.FFFFFEp-1f);
37 
38   const int32x4_t vindex_mask = vmovq_n_s32(INT32_C(0x7FF));
39 
40   for (; n >= 8 * sizeof(float); n -= 8 * sizeof(float)) {
41     const float32x4_t vx0123 = vld1q_f32(x); x += 4;
42     const float32x4_t vx4567 = vld1q_f32(x); x += 4;
43 
44     // General structure of the algorithm:
45     //           / exp(x) / (1 + exp(x)) if x <= 0
46     //   f[x] :=
47     //           \ 1 - f[-x] if x >= 0
48     //
49     // First we compute f[-z] := exp(-z) / (1 + exp(-z)) where z = abs(x),
50     // then replace result with 1 - f[-z] if x >= 0.
51     const float32x4_t vz0123 = vabsq_f32(vx0123);
52     const float32x4_t vz4567 = vabsq_f32(vx4567);
53 
54     // Compute reduced argument n := round(-z * 2048 / log(2)).
55     // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing
56     // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction.
57     // The trick with adding large number is valid only within certain bounds (|z * 2048 / log(2)| <= 2**22, i.e.
58     // |z| <= 0x1.62E43p+10 = 1419.5654296875), but that is acceptable, because inputs x outside of
59     // [-87.336544, 17.328678] (i.e. z outsize [0, 87.336544]) underflow or saturate sigmoidf(x). We fixup the result
60     // for such inputs at the very end of the algorithm.
61     float32x4_t vn0123 = vfmaq_f32(vmagic_bias, vz0123, vminus_log2e_x2048);
62     float32x4_t vn4567 = vfmaq_f32(vmagic_bias, vz4567, vminus_log2e_x2048);
63 
64     // Create a floating-point number s (scale) such that s := 2**(n / 2048) for such inputs that sigmoidf(-z) is
65     // normalized, i.e. 0 <= z <= 87.33642. As n has 11 fractional bits, we split s == 2**(n / 2048) =
66     // = 2**e * 2**(n / 2048 - e), where e := int(n / 2048). We create s in two steps:
67     // 1. Fetch 2**(n / 2048 - e) = 2**(n % 2048) from the table using the 6 low bits of n, as integer. Note that the
68     //    fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0.
69     // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized
70     //    number, because for 0 <= z <= 87.33642 (inputs for which sigmoidf(-z) is normalized) we have -126 <= e <= 0,
71     //    and thus the adjusted exponent is not lower than -126.
72     //
73     // Extract e from bits 11:19 of n and shift it into bits 23:31 (position of floating-point exponent).
74     const int32x4_t ve0123 = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn0123), vmovq_n_s32(INT32_C(0x7FF))), 12);
75     const int32x4_t ve4567 = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn4567), vmovq_n_s32(INT32_C(0x7FF))), 12);
76 
77     // Use bits 0:11 bits of n, as integer, as an index for table lookup of l := 2**(n % 2048).
78     const uint64x2_t vidx0123 = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn0123), vindex_mask));
79     const uint64x2_t vidx4567 = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn4567), vindex_mask));
80 
81     const uint64_t vidx01 = vgetq_lane_u64(vidx0123, 0);
82     const uint64_t vidx23 = vgetq_lane_u64(vidx0123, 1);
83     float32x2_t vl01 = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx01]);
84     float32x2_t vl23 = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx23]);
85     const uint64_t vidx45 = vgetq_lane_u64(vidx4567, 0);
86     const uint64_t vidx67 = vgetq_lane_u64(vidx4567, 1);
87     float32x2_t vl45 = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx45]);
88     float32x2_t vl67 = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx67]);
89 
90     vl01 = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx01 >> 32)], vl01, 1);
91     vl23 = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx23 >> 32)], vl23, 1);
92     const float32x4_t vl0123 = vcombine_f32(vl01, vl23);
93     vl45 = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx45 >> 32)], vl45, 1);
94     vl67 = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx67 >> 32)], vl67, 1);
95     const float32x4_t vl4567 = vcombine_f32(vl45, vl67);
96 
97     // Adjust exponent of the value l fetched from the table to get the final s value.
98     const float32x4_t vs0123 = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl0123), ve0123));
99     const float32x4_t vs4567 = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl4567), ve4567));
100 
101     // Subtract the large number back to get the final n := round(-z * 2048 / log(2)) as a floating-point number.
102     vn0123 = vsubq_f32(vn0123, vmagic_bias);
103     vn4567 = vsubq_f32(vn4567, vmagic_bias);
104 
105     // Compute reduced argument t := (z + n * log(2) / 2048). Note that -t = -z - n * log(2) / 2048.
106     float32x4_t vt0123 = vfmaq_f32(vz0123, vn0123, vln2_o2048);
107     float32x4_t vt4567 = vfmaq_f32(vz4567, vn4567, vln2_o2048);
108 
109     // Compute degree-1 polynomial approximation for exp(-t) on [-log(2)/2048, log(2)/2048]:
110     //   P1(t) = 1 + t * c1
111     const float32x4_t vp0123 = vmulq_f32(vt0123, vc1);
112     const float32x4_t vp4567 = vmulq_f32(vt4567, vc1);
113 
114     // Reconstruct the exp(-z) value:
115     //   y = s * (1 + t * c1)
116     //     = s + s * (t * c1))
117     //     = s + s * p
118     const float32x4_t vy0123 = vfmaq_f32(vs0123, vs0123, vp0123);
119     const float32x4_t vy4567 = vfmaq_f32(vs4567, vs4567, vp4567);
120 
121     // Denominator of the sigmoid fraction: 1.0 + exp(-z)
122     const float32x4_t vd0123 = vaddq_f32(vy0123, vone);
123     const float32x4_t vd4567 = vaddq_f32(vy4567, vone);
124 
125     // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z))
126     float32x4_t vf0123 = vdivq_f32(vy0123, vd0123);
127     float32x4_t vf4567 = vdivq_f32(vy4567, vd4567);
128 
129     // For inputs below denormal cutoff, replace output with +0.0f.
130     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
131     vf0123 = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf0123), vcagtq_f32(vx0123, vdenorm_cutoff)));
132     vf4567 = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf4567), vcagtq_f32(vx4567, vdenorm_cutoff)));
133 
134     // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z)
135     const uint32x4_t vm0123 = vcltq_f32(vx0123, vmovq_n_f32(0.0f));
136     const uint32x4_t vm4567 = vcltq_f32(vx4567, vmovq_n_f32(0.0f));
137 
138     vf0123 = vbslq_f32(vm0123, vf0123, vsubq_f32(vone, vf0123));
139     vf4567 = vbslq_f32(vm4567, vf4567, vsubq_f32(vone, vf4567));
140 
141     vst1q_f32(y, vf0123); y += 4;
142     vst1q_f32(y, vf4567); y += 4;
143   }
144   for (; n >= 4 * sizeof(float); n -= 4 * sizeof(float)) {
145     const float32x4_t vx = vld1q_f32(x); x += 4;
146 
147     // General structure of the algorithm:
148     //           / exp(x) / (1 + exp(x)) if x <= 0
149     //   f[x] :=
150     //           \ 1 - f[-x] if x >= 0
151     //
152     // First we compute f[-z] := exp(-z) / (1 + exp(-z)) where z = abs(x),
153     // then replace result with 1 - f[-z] if x >= 0.
154     const float32x4_t vz = vabsq_f32(vx);
155 
156     // Compute reduced argument n := round(-z * 2048 / log(2)).
157     // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing
158     // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction.
159     // The trick with adding large number is valid only within certain bounds (|z * 2048 / log(2)| <= 2**22, i.e.
160     // |z| <= 0x1.62E43p+10 = 1419.5654296875), but that is acceptable, because inputs x outside of
161     // [-87.336544, 17.328678] (i.e. z outsize [0, 87.336544]) underflow or saturate sigmoidf(x). We fixup the result
162     // for such inputs at the very end of the algorithm.
163     float32x4_t vn = vfmaq_f32(vmagic_bias, vz, vminus_log2e_x2048);
164 
165     // Create a floating-point number s (scale) such that s := 2**(n / 2048) for such inputs that sigmoidf(-z) is
166     // normalized, i.e. 0 <= z <= 87.33642. As n has 11 fractional bits, we split s == 2**(n / 2048) =
167     // = 2**e * 2**(n / 2048 - e), where e := int(n / 2048). We create s in two steps:
168     // 1. Fetch 2**(n / 2048 - e) = 2**(n % 2048) from exp2_k_over_2048_table using the 6 low bits of n, as integer. Note that the
169     //    fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0.
170     // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized
171     //    number, because for 0 <= z <= 87.33642 (inputs for which sigmoidf(-z) is normalized) we have -126 <= e <= 0,
172     //    and thus the adjusted exponent is not lower than -126.
173     //
174     // Extract e from bits 11:19 of n and shift it into bits 23:31 (position of floating-point exponent).
175     const int32x4_t ve = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn), vmovq_n_s32(INT32_C(0x7FF))), 12);
176 
177     // Use bits 0:11 bits of n, as integer, as an index for table lookup of l := 2**(n % 2048).
178     const uint64x2_t vidx = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn), vindex_mask));
179     const uint64_t vidx_lo = vgetq_lane_u64(vidx, 0);
180     const uint64_t vidx_hi = vgetq_lane_u64(vidx, 1);
181     float32x2_t vl_lo = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx_lo]);
182     float32x2_t vl_hi = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx_hi]);
183     vl_lo = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx_lo >> 32)], vl_lo, 1);
184     vl_hi = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx_hi >> 32)], vl_hi, 1);
185     const float32x4_t vl = vcombine_f32(vl_lo, vl_hi);
186     // Adjust exponent of the value l fetched from the exp2_k_over_2048_table to get the final s value.
187     const float32x4_t vs = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl), ve));
188 
189     // Subtract the large number back to get the final n := round(-z * 2048 / log(2)) as a floating-point number.
190     vn = vsubq_f32(vn, vmagic_bias);
191 
192     // Compute reduced argument t := (z + n * log(2) / 2048). Note that -t = -z - n * log(2) / 2048.
193     float32x4_t vt = vfmaq_f32(vz, vn, vln2_o2048);
194 
195     // Compute degree-1 polynomial approximation for exp(-t) on [-log(2)/2048, log(2)/2048]:
196     //   P1(t) = 1 + t * c1
197     const float32x4_t vp = vmulq_f32(vt, vc1);
198 
199     // Reconstruct the exp(-z) value:
200     //   y = s * (1 + t * c1)
201     //     = s + s * (t * c1))
202     //     = s + s * p
203     const float32x4_t vy = vfmaq_f32(vs, vs, vp);
204 
205     // Denominator of the sigmoid fraction: 1.0 + exp(-z)
206     const float32x4_t vd = vaddq_f32(vy, vone);
207 
208     // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z))
209     float32x4_t vf = vdivq_f32(vy, vd);
210 
211     // For inputs below denormal cutoff, replace output with +0.0f.
212     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
213     vf = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf), vcagtq_f32(vx, vdenorm_cutoff)));
214 
215     // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z)
216     const uint32x4_t vm = vcltq_f32(vx, vmovq_n_f32(0.0f));
217     vf = vbslq_f32(vm, vf, vsubq_f32(vone, vf));
218 
219     vst1q_f32(y, vf); y += 4;
220   }
221   if XNN_UNLIKELY(n != 0) {
222     const float32x4_t vx = vld1q_f32(x);
223 
224     // General structure of the algorithm:
225     //           / exp(x) / (1 + exp(x)) if x <= 0
226     //   f[x] :=
227     //           \ 1 - f[-x] if x >= 0
228     //
229     // First we compute f[-z] := exp(-z) / (1 + exp(-z)) where z = abs(x),
230     // then replace result with 1 - f[-z] if x >= 0.
231     const float32x4_t vz = vabsq_f32(vx);
232 
233     // Compute reduced argument n := round(-z * 2048 / log(2)).
234     // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing
235     // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction.
236     // The trick with adding large number is valid only within certain bounds (|z * 2048 / log(2)| <= 2**22, i.e.
237     // |z| <= 0x1.62E43p+10 = 1419.5654296875), but that is acceptable, because inputs x outside of
238     // [-87.336544, 17.328678] (i.e. z outsize [0, 87.336544]) underflow or saturate sigmoidf(x). We fixup the result
239     // for such inputs at the very end of the algorithm.
240     float32x4_t vn = vfmaq_f32(vmagic_bias, vz, vminus_log2e_x2048);
241 
242     // Create a floating-point number s (scale) such that s := 2**(n / 2048) for such inputs that sigmoidf(-z) is
243     // normalized, i.e. 0 <= z <= 87.33642. As n has 11 fractional bits, we split s == 2**(n / 2048) =
244     // = 2**e * 2**(n / 2048 - e), where e := int(n / 2048). We create s in two steps:
245     // 1. Fetch 2**(n / 2048 - e) = 2**(n % 2048) from exp2_k_over_2048_table using the 6 low bits of n, as integer. Note that the
246     //    fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0.
247     // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized
248     //    number, because for 0 <= z <= 87.33642 (inputs for which sigmoidf(-z) is normalized) we have -126 <= e <= 0,
249     //    and thus the adjusted exponent is not lower than -126.
250     //
251     // Extract e from bits 11:19 of n and shift it into bits 23:31 (position of floating-point exponent).
252     const int32x4_t ve = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn), vmovq_n_s32(INT32_C(0x7FF))), 12);
253 
254     // Use bits 0:11 bits of n, as integer, as an index for table lookup of l := 2**(n % 2048).
255     const uint64x2_t vidx = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn), vindex_mask));
256     const uint64_t vidx_lo = vgetq_lane_u64(vidx, 0);
257     const uint64_t vidx_hi = vgetq_lane_u64(vidx, 1);
258     float32x2_t vl_lo = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx_lo]);
259     float32x2_t vl_hi = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx_hi]);
260     vl_lo = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx_lo >> 32)], vl_lo, 1);
261     vl_hi = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx_hi >> 32)], vl_hi, 1);
262     const float32x4_t vl = vcombine_f32(vl_lo, vl_hi);
263     // Adjust exponent of the value l fetched from the exp2_k_over_2048_table to get the final s value.
264     const float32x4_t vs = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl), ve));
265 
266     // Subtract the large number back to get the final n := round(-z * 2048 / log(2)) as a floating-point number.
267     vn = vsubq_f32(vn, vmagic_bias);
268 
269     // Compute reduced argument t := (z + n * log(2) / 2048). Note that -t = -z - n * log(2) / 2048.
270     float32x4_t vt = vfmaq_f32(vz, vn, vln2_o2048);
271 
272     // Compute degree-1 polynomial approximation for exp(-t) on [-log(2)/2048, log(2)/2048]:
273     //   P1(t) = 1 + t * c1
274     const float32x4_t vp = vmulq_f32(vt, vc1);
275 
276     // Reconstruct the exp(-z) value:
277     //   y = s * (1 + t * c1)
278     //     = s + s * (t * c1))
279     //     = s + s * p
280     const float32x4_t vy = vfmaq_f32(vs, vs, vp);
281 
282     // Denominator of the sigmoid fraction: 1.0 + exp(-z)
283     const float32x4_t vd = vaddq_f32(vy, vone);
284 
285     // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z))
286     float32x4_t vf = vdivq_f32(vy, vd);
287 
288     // For inputs below denormal cutoff, replace output with +0.0f.
289     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
290     vf = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf), vcagtq_f32(vx, vdenorm_cutoff)));
291 
292     // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z)
293     const uint32x4_t vm = vcltq_f32(vx, vmovq_n_f32(0.0f));
294     vf = vbslq_f32(vm, vf, vsubq_f32(vone, vf));
295 
296     float32x2_t vf_lo = vget_low_f32(vf);
297     if (n & (2 * sizeof(float))) {
298       vst1_f32(y, vf_lo); y += 2;
299       vf_lo = vget_high_f32(vf);
300     }
301     if (n & (1 * sizeof(float))) {
302       vst1_lane_f32(y, vf_lo, 0);
303     }
304   }
305 }
306