xref: /aosp_15_r20/external/XNNPACK/src/math/sqrt-u32-scalar-clz-newton.c (revision 4bdc94577ba0e567308109d787f7fec7b531ce36)
1 // Copyright 2022 Google LLC
2 //
3 // This source code is licensed under the BSD-style license found in the
4 // LICENSE file in the root directory of this source tree.
5 
6 #include <assert.h>
7 #include <stddef.h>
8 
9 #include <xnnpack/math.h>
10 #include <xnnpack/math-stubs.h>
11 
12 
xnn_math_u32_sqrt__scalar_clz_newton(size_t n,const uint32_t * input,uint32_t * output)13 void xnn_math_u32_sqrt__scalar_clz_newton(
14     size_t n,
15     const uint32_t* input,
16     uint32_t* output)
17 {
18   assert(n % sizeof(uint32_t) == 0);
19 
20   for (; n != 0; n -= sizeof(uint32_t)) {
21     const uint32_t vx = *input++;
22 
23     uint32_t vy = vx;
24 
25     // Based on Hacker's Delight, Figure 11-1.
26     if (vx != 0) {
27       const uint32_t vs = 16 - (math_clz_nonzero_u32(vx - 1) >> 1);
28 
29       uint32_t vg0 = UINT32_C(1) << vs;
30       uint32_t vg1 = (vg0 + (vx >> vs)) >> 1;
31       while XNN_LIKELY(vg1 < vg0) {
32         vg0 = vg1;
33         vg1 = (vg0 + vx / vg0) >> 1;
34       }
35 
36       // vg0 is sqrt(vx) rounded down. Do the final rounding up if needed.
37       if XNN_UNPREDICTABLE(vg0 * vg0 < vx - vg0) {
38         vg0 += 1;
39       }
40       vy = vg0;
41     }
42 
43     *output++ = vy;
44   }
45 }
46