• 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_nr1recps1fma_x8(size_t n,const float * x,float * y,const void * params)20 void xnn_f32_sigmoid_ukernel__neonfma_rr1_lut2048_p1_nr1recps1fma_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     // Use Newton-Raphson method (2 iterations) to compute reciprocal of denominator.
126     // Note: 1 < d <= 2, because z >= 0.0 and 0 < exp(-z) <= 1.0.
127     // Thus the reciprocal of the denominator never overflows.
128     float32x4_t vr0123 = vrecpeq_f32(vd0123);
129     float32x4_t vr4567 = vrecpeq_f32(vd4567);
130 
131     vr0123 = vmulq_f32(vr0123, vrecpsq_f32(vr0123, vd0123));
132     vr4567 = vmulq_f32(vr4567, vrecpsq_f32(vr4567, vd4567));
133 
134     vr0123 = vfmaq_f32(vr0123, vr0123, vfmsq_f32(vone, vr0123, vd0123));
135     vr4567 = vfmaq_f32(vr4567, vr4567, vfmsq_f32(vone, vr4567, vd4567));
136 
137     // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z))
138     float32x4_t vf0123 = vmulq_f32(vy0123, vr0123);
139     float32x4_t vf4567 = vmulq_f32(vy4567, vr4567);
140 
141     // For inputs below denormal cutoff, replace output with +0.0f.
142     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
143     vf0123 = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf0123), vcagtq_f32(vx0123, vdenorm_cutoff)));
144     vf4567 = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf4567), vcagtq_f32(vx4567, vdenorm_cutoff)));
145 
146     // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z)
147     const uint32x4_t vm0123 = vcltq_f32(vx0123, vmovq_n_f32(0.0f));
148     const uint32x4_t vm4567 = vcltq_f32(vx4567, vmovq_n_f32(0.0f));
149 
150     vf0123 = vbslq_f32(vm0123, vf0123, vsubq_f32(vone, vf0123));
151     vf4567 = vbslq_f32(vm4567, vf4567, vsubq_f32(vone, vf4567));
152 
153     vst1q_f32(y, vf0123); y += 4;
154     vst1q_f32(y, vf4567); y += 4;
155   }
156   for (; n >= 4 * sizeof(float); n -= 4 * sizeof(float)) {
157     const float32x4_t vx = vld1q_f32(x); x += 4;
158 
159     // General structure of the algorithm:
160     //           / exp(x) / (1 + exp(x)) if x <= 0
161     //   f[x] :=
162     //           \ 1 - f[-x] if x >= 0
163     //
164     // First we compute f[-z] := exp(-z) / (1 + exp(-z)) where z = abs(x),
165     // then replace result with 1 - f[-z] if x >= 0.
166     const float32x4_t vz = vabsq_f32(vx);
167 
168     // Compute reduced argument n := round(-z * 2048 / log(2)).
169     // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing
170     // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction.
171     // The trick with adding large number is valid only within certain bounds (|z * 2048 / log(2)| <= 2**22, i.e.
172     // |z| <= 0x1.62E43p+10 = 1419.5654296875), but that is acceptable, because inputs x outside of
173     // [-87.336544, 17.328678] (i.e. z outsize [0, 87.336544]) underflow or saturate sigmoidf(x). We fixup the result
174     // for such inputs at the very end of the algorithm.
175     float32x4_t vn = vfmaq_f32(vmagic_bias, vz, vminus_log2e_x2048);
176 
177     // Create a floating-point number s (scale) such that s := 2**(n / 2048) for such inputs that sigmoidf(-z) is
178     // normalized, i.e. 0 <= z <= 87.33642. As n has 11 fractional bits, we split s == 2**(n / 2048) =
179     // = 2**e * 2**(n / 2048 - e), where e := int(n / 2048). We create s in two steps:
180     // 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
181     //    fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0.
182     // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized
183     //    number, because for 0 <= z <= 87.33642 (inputs for which sigmoidf(-z) is normalized) we have -126 <= e <= 0,
184     //    and thus the adjusted exponent is not lower than -126.
185     //
186     // Extract e from bits 11:19 of n and shift it into bits 23:31 (position of floating-point exponent).
187     const int32x4_t ve = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn), vmovq_n_s32(INT32_C(0x7FF))), 12);
188 
189     // Use bits 0:11 bits of n, as integer, as an index for table lookup of l := 2**(n % 2048).
190     const uint64x2_t vidx = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn), vindex_mask));
191     const uint64_t vidx_lo = vgetq_lane_u64(vidx, 0);
192     const uint64_t vidx_hi = vgetq_lane_u64(vidx, 1);
193     float32x2_t vl_lo = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx_lo]);
194     float32x2_t vl_hi = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx_hi]);
195     vl_lo = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx_lo >> 32)], vl_lo, 1);
196     vl_hi = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx_hi >> 32)], vl_hi, 1);
197     const float32x4_t vl = vcombine_f32(vl_lo, vl_hi);
198     // Adjust exponent of the value l fetched from the exp2_k_over_2048_table to get the final s value.
199     const float32x4_t vs = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl), ve));
200 
201     // Subtract the large number back to get the final n := round(-z * 2048 / log(2)) as a floating-point number.
202     vn = vsubq_f32(vn, vmagic_bias);
203 
204     // Compute reduced argument t := (z + n * log(2) / 2048). Note that -t = -z - n * log(2) / 2048.
205     float32x4_t vt = vfmaq_f32(vz, vn, vln2_o2048);
206 
207     // Compute degree-1 polynomial approximation for exp(-t) on [-log(2)/2048, log(2)/2048]:
208     //   P1(t) = 1 + t * c1
209     const float32x4_t vp = vmulq_f32(vt, vc1);
210 
211     // Reconstruct the exp(-z) value:
212     //   y = s * (1 + t * c1)
213     //     = s + s * (t * c1))
214     //     = s + s * p
215     const float32x4_t vy = vfmaq_f32(vs, vs, vp);
216 
217     // Denominator of the sigmoid fraction: 1.0 + exp(-z)
218     const float32x4_t vd = vaddq_f32(vy, vone);
219 
220     // Use Newton-Raphson method (2 iterations) to compute reciprocal of denominator.
221     // Note: 1 < d <= 2, because z >= 0.0 and 0 < exp(-z) <= 1.0.
222     // Thus the reciprocal of the denominator never overflows.
223     float32x4_t vr = vrecpeq_f32(vd);
224 
225     vr = vmulq_f32(vr, vrecpsq_f32(vr, vd));
226 
227     vr = vfmaq_f32(vr, vr, vfmsq_f32(vone, vr, vd));
228 
229     // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z))
230     float32x4_t vf = vmulq_f32(vy, vr);
231 
232     // For inputs below denormal cutoff, replace output with +0.0f.
233     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
234     vf = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf), vcagtq_f32(vx, vdenorm_cutoff)));
235 
236     // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z)
237     const uint32x4_t vm = vcltq_f32(vx, vmovq_n_f32(0.0f));
238     vf = vbslq_f32(vm, vf, vsubq_f32(vone, vf));
239 
240     vst1q_f32(y, vf); y += 4;
241   }
242   if XNN_UNLIKELY(n != 0) {
243     const float32x4_t vx = vld1q_f32(x);
244 
245     // General structure of the algorithm:
246     //           / exp(x) / (1 + exp(x)) if x <= 0
247     //   f[x] :=
248     //           \ 1 - f[-x] if x >= 0
249     //
250     // First we compute f[-z] := exp(-z) / (1 + exp(-z)) where z = abs(x),
251     // then replace result with 1 - f[-z] if x >= 0.
252     const float32x4_t vz = vabsq_f32(vx);
253 
254     // Compute reduced argument n := round(-z * 2048 / log(2)).
255     // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing
256     // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction.
257     // The trick with adding large number is valid only within certain bounds (|z * 2048 / log(2)| <= 2**22, i.e.
258     // |z| <= 0x1.62E43p+10 = 1419.5654296875), but that is acceptable, because inputs x outside of
259     // [-87.336544, 17.328678] (i.e. z outsize [0, 87.336544]) underflow or saturate sigmoidf(x). We fixup the result
260     // for such inputs at the very end of the algorithm.
261     float32x4_t vn = vfmaq_f32(vmagic_bias, vz, vminus_log2e_x2048);
262 
263     // Create a floating-point number s (scale) such that s := 2**(n / 2048) for such inputs that sigmoidf(-z) is
264     // normalized, i.e. 0 <= z <= 87.33642. As n has 11 fractional bits, we split s == 2**(n / 2048) =
265     // = 2**e * 2**(n / 2048 - e), where e := int(n / 2048). We create s in two steps:
266     // 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
267     //    fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0.
268     // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized
269     //    number, because for 0 <= z <= 87.33642 (inputs for which sigmoidf(-z) is normalized) we have -126 <= e <= 0,
270     //    and thus the adjusted exponent is not lower than -126.
271     //
272     // Extract e from bits 11:19 of n and shift it into bits 23:31 (position of floating-point exponent).
273     const int32x4_t ve = vshlq_n_s32(vbicq_s32(vreinterpretq_s32_f32(vn), vmovq_n_s32(INT32_C(0x7FF))), 12);
274 
275     // Use bits 0:11 bits of n, as integer, as an index for table lookup of l := 2**(n % 2048).
276     const uint64x2_t vidx = vreinterpretq_u64_s32(vandq_s32(vreinterpretq_s32_f32(vn), vindex_mask));
277     const uint64_t vidx_lo = vgetq_lane_u64(vidx, 0);
278     const uint64_t vidx_hi = vgetq_lane_u64(vidx, 1);
279     float32x2_t vl_lo = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx_lo]);
280     float32x2_t vl_hi = vld1_dup_f32(&xnn_table_exp2_k_over_2048[(uint32_t) vidx_hi]);
281     vl_lo = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx_lo >> 32)], vl_lo, 1);
282     vl_hi = vld1_lane_f32(&xnn_table_exp2_k_over_2048[(uint32_t) (vidx_hi >> 32)], vl_hi, 1);
283     const float32x4_t vl = vcombine_f32(vl_lo, vl_hi);
284     // Adjust exponent of the value l fetched from the exp2_k_over_2048_table to get the final s value.
285     const float32x4_t vs = vreinterpretq_f32_s32(vaddq_s32(vreinterpretq_s32_f32(vl), ve));
286 
287     // Subtract the large number back to get the final n := round(-z * 2048 / log(2)) as a floating-point number.
288     vn = vsubq_f32(vn, vmagic_bias);
289 
290     // Compute reduced argument t := (z + n * log(2) / 2048). Note that -t = -z - n * log(2) / 2048.
291     float32x4_t vt = vfmaq_f32(vz, vn, vln2_o2048);
292 
293     // Compute degree-1 polynomial approximation for exp(-t) on [-log(2)/2048, log(2)/2048]:
294     //   P1(t) = 1 + t * c1
295     const float32x4_t vp = vmulq_f32(vt, vc1);
296 
297     // Reconstruct the exp(-z) value:
298     //   y = s * (1 + t * c1)
299     //     = s + s * (t * c1))
300     //     = s + s * p
301     const float32x4_t vy = vfmaq_f32(vs, vs, vp);
302 
303     // Denominator of the sigmoid fraction: 1.0 + exp(-z)
304     const float32x4_t vd = vaddq_f32(vy, vone);
305 
306     // Use Newton-Raphson method (2 iterations) to compute reciprocal of denominator.
307     // Note: 1 < d <= 2, because z >= 0.0 and 0 < exp(-z) <= 1.0.
308     // Thus the reciprocal of the denominator never overflows.
309     float32x4_t vr = vrecpeq_f32(vd);
310 
311     vr = vmulq_f32(vr, vrecpsq_f32(vr, vd));
312 
313     vr = vfmaq_f32(vr, vr, vfmsq_f32(vone, vr, vd));
314 
315     // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z))
316     float32x4_t vf = vmulq_f32(vy, vr);
317 
318     // For inputs below denormal cutoff, replace output with +0.0f.
319     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
320     vf = vreinterpretq_f32_u32(vbicq_u32(vreinterpretq_u32_f32(vf), vcagtq_f32(vx, vdenorm_cutoff)));
321 
322     // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z)
323     const uint32x4_t vm = vcltq_f32(vx, vmovq_n_f32(0.0f));
324     vf = vbslq_f32(vm, vf, vsubq_f32(vone, vf));
325 
326     float32x2_t vf_lo = vget_low_f32(vf);
327     if (n & (2 * sizeof(float))) {
328       vst1_f32(y, vf_lo); y += 2;
329       vf_lo = vget_high_f32(vf);
330     }
331     if (n & (1 * sizeof(float))) {
332       vst1_lane_f32(y, vf_lo, 0);
333     }
334   }
335 }
336