Skip to content

Commit fbcf972

Browse files
falseywinchnetfalseywinchnet
authored
Merge pull request #64 from falseywinchnet/cdq5hs-main/add-new-atan_vectoring_network_experiment
Add vec5+cubic vectoring phase path and benchmark comparisons
2 parents 6c5a90a + ca78ab9 commit fbcf972

2 files changed

Lines changed: 261 additions & 9 deletions

File tree

‎examples/atan2_benchmark.cpp‎

Lines changed: 31 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,12 @@ int main(int argc, char** argv) {
7777
return normalize_phase(std::atan2(value.y, value.x));
7878
}, sink);
7979

80-
const double bfft_double_seconds = time_loop(samples, [](const sample& value) {
80+
const double bfft_tree_double_seconds = time_loop(samples, [](const sample& value) {
81+
return bruun::bruun_phase_atan2_core(value.y, value.x, value.mag,
82+
bruun::bruun_phase_first_octant_tree32_degree7);
83+
}, sink);
84+
85+
const double bfft_vec_double_seconds = time_loop(samples, [](const sample& value) {
8186
return bruun::bruun_phase_atan2_mag(value.y, value.x, value.mag);
8287
}, sink);
8388

@@ -87,18 +92,40 @@ int main(int argc, char** argv) {
8792
return normalize_phase_f32(std::atan2(y, x));
8893
}, sink);
8994

90-
const double bfft_float_seconds = time_loop(samples, [](const sample& value) {
95+
const double bfft_tree_float_seconds = time_loop(samples, [](const sample& value) {
96+
const double phase = bruun::bruun_phase_atan2_core(value.y, value.x, value.mag,
97+
bruun::bruun_phase_first_octant_tree32_degree7);
98+
return static_cast<float>(phase);
99+
}, sink);
100+
101+
const double bfft_vec_float_seconds = time_loop(samples, [](const sample& value) {
91102
const float y = static_cast<float>(value.y);
92103
const float x = static_cast<float>(value.x);
93104
const float mag = static_cast<float>(value.mag);
94105
return bruun::bruun_phase_atan2_mag_f32(y, x, mag);
95106
}, sink);
96107

108+
const double std_sincos_seconds = time_loop(samples, [](const sample& value) {
109+
return std::sin(value.x) + std::cos(value.x);
110+
}, sink);
111+
112+
const double bfft_sincos_seconds = time_loop(samples, [](const sample& value) {
113+
double s_value;
114+
double c_value;
115+
double phase = std::fmod(std::fabs(value.x), tau);
116+
bruun::bruun_table256_poly3_sincos(phase, &s_value, &c_value);
117+
return s_value + c_value;
118+
}, sink);
119+
97120
const double scale = static_cast<double>(count) / 1000000.0;
98121
std::printf("samples=%zu sink=%.17g\n", count, sink);
99122
std::printf("std::atan2 double: %.6f s, %.3f Mphase/s\n", std_double_seconds, scale / std_double_seconds);
100-
std::printf("bfft atan2 double: %.6f s, %.3f Mphase/s\n", bfft_double_seconds, scale / bfft_double_seconds);
123+
std::printf("bfft tree32 double: %.6f s, %.3f Mphase/s\n", bfft_tree_double_seconds, scale / bfft_tree_double_seconds);
124+
std::printf("bfft vec5 double : %.6f s, %.3f Mphase/s\n", bfft_vec_double_seconds, scale / bfft_vec_double_seconds);
101125
std::printf("std::atan2 float : %.6f s, %.3f Mphase/s\n", std_float_seconds, scale / std_float_seconds);
102-
std::printf("bfft atan2 float : %.6f s, %.3f Mphase/s\n", bfft_float_seconds, scale / bfft_float_seconds);
126+
std::printf("bfft tree32 float : %.6f s, %.3f Mphase/s\n", bfft_tree_float_seconds, scale / bfft_tree_float_seconds);
127+
std::printf("bfft vec5 float : %.6f s, %.3f Mphase/s\n", bfft_vec_float_seconds, scale / bfft_vec_float_seconds);
128+
std::printf("std sin+cos : %.6f s, %.3f Mpair/s\n", std_sincos_seconds, scale / std_sincos_seconds);
129+
std::printf("bfft table sincos : %.6f s, %.3f Mpair/s\n", bfft_sincos_seconds, scale / bfft_sincos_seconds);
103130
return 0;
104131
}

‎src/detail/bruun_kernel.hpp‎

Lines changed: 230 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,7 @@ typedef __m128d bruun_v2;
5959
# define V2_ADD(a, b) _mm_add_pd((a), (b))
6060
# define V2_SUB(a, b) _mm_sub_pd((a), (b))
6161
# define V2_MUL(a, b) _mm_mul_pd((a), (b))
62+
# define V2_DIV(a, b) _mm_div_pd((a), (b))
6263
# define V2_MADD(a, b, c) V2_ADD((a), V2_MUL((b), (c)))
6364
# define V2_MSUB(a, b, c) V2_SUB((a), V2_MUL((b), (c)))
6465
# define V2_SET1(x) _mm_set1_pd(x)
@@ -68,13 +69,16 @@ typedef __m128d bruun_v2;
6869
# define V2_DUP0(a) _mm_unpacklo_pd((a), (a))
6970
# define V2_DUP1(a) _mm_unpackhi_pd((a), (a))
7071
# define V2_NEGHI(a) _mm_xor_pd((a), _mm_set_pd(-0.0, 0.0))
72+
# define V2_CMPGT(a, b) _mm_cmpgt_pd((a), (b))
73+
# define V2_SELECT(m, t, f) _mm_or_pd(_mm_and_pd((m), (t)), _mm_andnot_pd((m), (f)))
7174
#elif defined(BRUUN_NEON_128)
7275
typedef float64x2_t bruun_v2;
7376
# define V2_LD(p) vld1q_f64(p)
7477
# define V2_ST(p, a) vst1q_f64((p), (a))
7578
# define V2_ADD(a, b) vaddq_f64((a), (b))
7679
# define V2_SUB(a, b) vsubq_f64((a), (b))
7780
# define V2_MUL(a, b) vmulq_f64((a), (b))
81+
# define V2_DIV(a, b) vdivq_f64((a), (b))
7882
# define V2_MADD(a, b, c) vfmaq_f64((a), (b), (c))
7983
# define V2_MSUB(a, b, c) vfmsq_f64((a), (b), (c))
8084
# define V2_SET1(x) vdupq_n_f64(x)
@@ -88,6 +92,14 @@ static inline float64x2_t bruun_neghi(float64x2_t a) {
8892
return vreinterpretq_f64_u64(veorq_u64(vreinterpretq_u64_f64(a), m));
8993
}
9094
# define V2_NEGHI(a) bruun_neghi(a)
95+
static inline float64x2_t bruun_v2_cmpgt(float64x2_t a, float64x2_t b) {
96+
return vreinterpretq_f64_u64(vcgtq_f64(a, b));
97+
}
98+
static inline float64x2_t bruun_v2_select(float64x2_t m, float64x2_t t, float64x2_t f) {
99+
return vbslq_f64(vreinterpretq_u64_f64(m), t, f);
100+
}
101+
# define V2_CMPGT(a, b) bruun_v2_cmpgt((a), (b))
102+
# define V2_SELECT(m, t, f) bruun_v2_select((m), (t), (f))
91103
#endif
92104

93105
// 4-lane float vector primitive for internal float32 helper loops.
@@ -151,6 +163,12 @@ namespace bruun {
151163
#define BRUUN_ASSUME_ALIGNED(ptr, bytes) (ptr)
152164
#endif
153165

166+
#if defined(__GNUC__) || defined(__clang__)
167+
#define BRUUN_ALWAYS_INLINE inline __attribute__((always_inline))
168+
#else
169+
#define BRUUN_ALWAYS_INLINE inline
170+
#endif
171+
154172
static constexpr std::size_t bruun_cache_alignment = 64;
155173

156174
template <typename T>
@@ -463,6 +481,11 @@ static inline double bruun_atan_leaf_odd(double h) {
463481
return h * p;
464482
}
465483

484+
static BRUUN_ALWAYS_INLINE double bruun_atan_leaf_cubic(double h) {
485+
const double s = h * h;
486+
return h * (1.0 - s * (1.0 / 3.0));
487+
}
488+
466489
static inline int bruun_phase_child_index(int idx, int base, double major, double minor) {
467490
const bruun_phase_node node = bruun_phase_nodes[base + idx];
468491
const double side = bruun_fused_msub(minor, node.c, major, node.s);
@@ -473,7 +496,7 @@ static inline int bruun_phase_child_index(int idx, int base, double major, doubl
473496
return idx;
474497
}
475498

476-
static inline double bruun_phase_first_octant(double major, double minor, double mag) {
499+
static inline double bruun_phase_first_octant_tree32_degree7(double major, double minor, double mag) {
477500
int idx = 0;
478501
idx = bruun_phase_child_index(idx, 0, major, minor);
479502
idx = bruun_phase_child_index(idx, 1, major, minor);
@@ -489,7 +512,74 @@ static inline double bruun_phase_first_octant(double major, double minor, double
489512
return leaf.center + delta;
490513
}
491514

492-
static inline double bruun_phase_atan2_mag(double y, double x, double mag) {
515+
static BRUUN_ALWAYS_INLINE double bruun_phase_first_octant_vec5_cubic(double major, double minor, double mag) {
516+
double theta = M_PI / 8.0;
517+
double c = 0.92387953251128674;
518+
double s = 0.38268343236508978;
519+
520+
const double side0 = bruun_fused_msub(minor, c, major, s);
521+
const double negative0 = static_cast<double>(side0 <= 0.0);
522+
const double sign0 = 1.0 - 2.0 * negative0;
523+
theta = bruun_fused_madd(sign0, 0.19634954084936207, theta);
524+
const double signed_sstep0 = sign0 * 0.19509032201612825;
525+
const double nc0 = bruun_fused_madd(-signed_sstep0, s, c * 0.98078528040323043);
526+
const double ns0 = bruun_fused_madd(signed_sstep0, c, s * 0.98078528040323043);
527+
c = nc0;
528+
s = ns0;
529+
530+
const double side1 = bruun_fused_msub(minor, c, major, s);
531+
const double negative1 = static_cast<double>(side1 <= 0.0);
532+
const double sign1 = 1.0 - 2.0 * negative1;
533+
theta = bruun_fused_madd(sign1, 0.09817477042468103, theta);
534+
const double signed_sstep1 = sign1 * 0.098017140329560604;
535+
const double nc1 = bruun_fused_madd(-signed_sstep1, s, c * 0.99518472667219693);
536+
const double ns1 = bruun_fused_madd(signed_sstep1, c, s * 0.99518472667219693);
537+
c = nc1;
538+
s = ns1;
539+
540+
const double side2 = bruun_fused_msub(minor, c, major, s);
541+
const double negative2 = static_cast<double>(side2 <= 0.0);
542+
const double sign2 = 1.0 - 2.0 * negative2;
543+
theta = bruun_fused_madd(sign2, 0.04908738521234052, theta);
544+
const double signed_sstep2 = sign2 * 0.049067674327418015;
545+
const double nc2 = bruun_fused_madd(-signed_sstep2, s, c * 0.99879545620517241);
546+
const double ns2 = bruun_fused_madd(signed_sstep2, c, s * 0.99879545620517241);
547+
c = nc2;
548+
s = ns2;
549+
550+
const double side3 = bruun_fused_msub(minor, c, major, s);
551+
const double negative3 = static_cast<double>(side3 <= 0.0);
552+
const double sign3 = 1.0 - 2.0 * negative3;
553+
theta = bruun_fused_madd(sign3, 0.02454369260617026, theta);
554+
const double signed_sstep3 = sign3 * 0.024541228522912288;
555+
const double nc3 = bruun_fused_madd(-signed_sstep3, s, c * 0.99969881869620425);
556+
const double ns3 = bruun_fused_madd(signed_sstep3, c, s * 0.99969881869620425);
557+
c = nc3;
558+
s = ns3;
559+
560+
const double side4 = bruun_fused_msub(minor, c, major, s);
561+
const double negative4 = static_cast<double>(side4 <= 0.0);
562+
const double sign4 = 1.0 - 2.0 * negative4;
563+
theta = bruun_fused_madd(sign4, 0.01227184630308513, theta);
564+
const double signed_sstep4 = sign4 * 0.012271538285719925;
565+
const double nc4 = bruun_fused_madd(-signed_sstep4, s, c * 0.9999247018391445);
566+
const double ns4 = bruun_fused_madd(signed_sstep4, c, s * 0.9999247018391445);
567+
c = nc4;
568+
s = ns4;
569+
570+
const double dot = bruun_fused_madd(major, c, minor * s);
571+
const double cross = bruun_fused_msub(minor, c, major, s);
572+
const double h = cross / (mag + dot);
573+
const double delta = 2.0 * bruun_atan_leaf_cubic(h);
574+
return theta + delta;
575+
}
576+
577+
static BRUUN_ALWAYS_INLINE double bruun_phase_first_octant(double major, double minor, double mag) {
578+
return bruun_phase_first_octant_vec5_cubic(major, minor, mag);
579+
}
580+
581+
template <class FirstOctant>
582+
static inline double bruun_phase_atan2_core(double y, double x, double mag, FirstOctant first_octant) {
493583
const double ax = std::fabs(x);
494584
const double ay = std::fabs(y);
495585
if (ay == 0.0) {
@@ -510,9 +600,9 @@ static inline double bruun_phase_atan2_mag(double y, double x, double mag) {
510600

511601
double alpha = 0.0;
512602
if (ax >= ay) {
513-
alpha = bruun_phase_first_octant(ax, ay, mag);
603+
alpha = first_octant(ax, ay, mag);
514604
} else {
515-
alpha = bruun_pio2 - bruun_phase_first_octant(ay, ax, mag);
605+
alpha = bruun_pio2 - first_octant(ay, ax, mag);
516606
}
517607

518608
if (x < 0.0) {
@@ -527,6 +617,121 @@ static inline double bruun_phase_atan2_mag(double y, double x, double mag) {
527617
return alpha;
528618
}
529619

620+
621+
#if BRUUN_LEVEL >= 1
622+
static BRUUN_ALWAYS_INLINE bruun_v2 bruun_phase_first_octant_vec5_cubic_v2(bruun_v2 major,
623+
bruun_v2 minor,
624+
bruun_v2 mag) {
625+
bruun_v2 theta = V2_SET1(M_PI / 8.0);
626+
bruun_v2 c = V2_SET1(0.92387953251128674);
627+
bruun_v2 s = V2_SET1(0.38268343236508978);
628+
const bruun_v2 one = V2_SET1(1.0);
629+
const bruun_v2 neg_one = V2_SET1(-1.0);
630+
631+
#define BRUUN_VEC5_V2_STEP(STEP, CSTEP, SSTEP) \
632+
do { \
633+
const bruun_v2 side = V2_MSUB(V2_MUL(minor, c), major, s); \
634+
const bruun_v2 sign = V2_SELECT(V2_CMPGT(side, V2_SET1(0.0)), one, neg_one); \
635+
theta = V2_ADD(theta, V2_MUL(sign, V2_SET1(STEP))); \
636+
const bruun_v2 signed_sstep = V2_MUL(sign, V2_SET1(SSTEP)); \
637+
const bruun_v2 nc = V2_SUB(V2_MUL(c, V2_SET1(CSTEP)), V2_MUL(s, signed_sstep)); \
638+
const bruun_v2 ns = V2_ADD(V2_MUL(c, signed_sstep), V2_MUL(s, V2_SET1(CSTEP))); \
639+
c = nc; \
640+
s = ns; \
641+
} while (false)
642+
643+
BRUUN_VEC5_V2_STEP(0.19634954084936207, 0.98078528040323043, 0.19509032201612825);
644+
BRUUN_VEC5_V2_STEP(0.09817477042468103, 0.99518472667219693, 0.098017140329560604);
645+
BRUUN_VEC5_V2_STEP(0.04908738521234052, 0.99879545620517241, 0.049067674327418015);
646+
BRUUN_VEC5_V2_STEP(0.02454369260617026, 0.99969881869620425, 0.024541228522912288);
647+
BRUUN_VEC5_V2_STEP(0.01227184630308513, 0.9999247018391445, 0.012271538285719925);
648+
#undef BRUUN_VEC5_V2_STEP
649+
650+
const bruun_v2 dot = V2_ADD(V2_MUL(major, c), V2_MUL(minor, s));
651+
const bruun_v2 cross = V2_SUB(V2_MUL(minor, c), V2_MUL(major, s));
652+
const bruun_v2 h = V2_DIV(cross, V2_ADD(mag, dot));
653+
const bruun_v2 h2 = V2_MUL(h, h);
654+
const bruun_v2 cubic = V2_MUL(h, V2_SUB(one, V2_MUL(h2, V2_SET1(1.0 / 3.0))));
655+
return V2_ADD(theta, V2_MUL(V2_SET1(2.0), cubic));
656+
}
657+
#endif
658+
659+
660+
#if BRUUN_LEVEL >= 1
661+
static BRUUN_ALWAYS_INLINE void bruun_phase_atan2_mag_pair(double y0, double x0, double mag0,
662+
double y1, double x1, double mag1,
663+
double* out0, double* out1) {
664+
const double ax0 = std::fabs(x0);
665+
const double ay0 = std::fabs(y0);
666+
const double ax1 = std::fabs(x1);
667+
const double ay1 = std::fabs(y1);
668+
if (ax0 == 0.0 || ay0 == 0.0 || ax1 == 0.0 || ay1 == 0.0) {
669+
*out0 = bruun_phase_atan2_core(y0, x0, mag0, bruun_phase_first_octant);
670+
*out1 = bruun_phase_atan2_core(y1, x1, mag1, bruun_phase_first_octant);
671+
return;
672+
}
673+
674+
double major0 = ax0;
675+
double minor0 = ay0;
676+
bool swap0 = false;
677+
if (ay0 > ax0) {
678+
major0 = ay0;
679+
minor0 = ax0;
680+
swap0 = true;
681+
}
682+
683+
double major1 = ax1;
684+
double minor1 = ay1;
685+
bool swap1 = false;
686+
if (ay1 > ax1) {
687+
major1 = ay1;
688+
minor1 = ax1;
689+
swap1 = true;
690+
}
691+
692+
const bruun_v2 major = V2_SETLH(major0, major1);
693+
const bruun_v2 minor = V2_SETLH(minor0, minor1);
694+
const bruun_v2 mag = V2_SETLH(mag0, mag1);
695+
double phase_pair[2];
696+
V2_ST(phase_pair, bruun_phase_first_octant_vec5_cubic_v2(major, minor, mag));
697+
698+
double alpha0 = phase_pair[0];
699+
if (swap0) {
700+
alpha0 = bruun_pio2 - alpha0;
701+
}
702+
if (x0 < 0.0) {
703+
alpha0 = M_PI - alpha0;
704+
}
705+
if (y0 < 0.0) {
706+
alpha0 = -alpha0;
707+
}
708+
if (alpha0 < 0.0) {
709+
alpha0 += bruun_tau;
710+
}
711+
712+
double alpha1 = phase_pair[1];
713+
if (swap1) {
714+
alpha1 = bruun_pio2 - alpha1;
715+
}
716+
if (x1 < 0.0) {
717+
alpha1 = M_PI - alpha1;
718+
}
719+
if (y1 < 0.0) {
720+
alpha1 = -alpha1;
721+
}
722+
if (alpha1 < 0.0) {
723+
alpha1 += bruun_tau;
724+
}
725+
726+
*out0 = alpha0;
727+
*out1 = alpha1;
728+
}
729+
#endif
730+
731+
static BRUUN_ALWAYS_INLINE double bruun_phase_atan2_mag(double y, double x, double mag) {
732+
return bruun_phase_atan2_core(y, x, mag, bruun_phase_first_octant);
733+
}
734+
530735
static inline float bruun_phase_atan2_mag_f32(float y, float x, float mag) {
531736
const double phase = bruun_phase_atan2_mag(static_cast<double>(y),
532737
static_cast<double>(x),
@@ -2531,7 +2736,27 @@ class RFFT {
25312736
}
25322737

25332738
const int* RESTRICT kin = KINV.data();
2534-
for (int k = 1; k < N / 2; ++k) {
2739+
int k = 1;
2740+
#if BRUUN_LEVEL >= 1
2741+
for (; k + 1 < N / 2; k += 2) {
2742+
const int m0 = kin[k];
2743+
const int m1 = kin[k + 1];
2744+
const double re0 = work[2*m0];
2745+
const double im0 = -work[2*m0 + 1];
2746+
const double re1 = work[2*m1];
2747+
const double im1 = -work[2*m1 + 1];
2748+
const double mag0 = std::sqrt(re0 * re0 + im0 * im0);
2749+
const double mag1 = std::sqrt(re1 * re1 + im1 * im1);
2750+
double phase0;
2751+
double phase1;
2752+
bruun_phase_atan2_mag_pair(im0, re0, mag0, im1, re1, mag1, &phase0, &phase1);
2753+
output[k].re = mag0;
2754+
output[k].im = phase0;
2755+
output[k + 1].re = mag1;
2756+
output[k + 1].im = phase1;
2757+
}
2758+
#endif
2759+
for (; k < N / 2; ++k) {
25352760
const int m = kin[k];
25362761
const double re = work[2*m];
25372762
const double im = -work[2*m + 1];

0 commit comments

Comments
 (0)