在使用汇编实现 cosf 时(参考math-neon library),遇到:调用cosf_neon,得到的值与sinf_neon的值一致。例如:
假设参数为0.366519,C math 得到的正确的值为 cos(0.366519)=0.933580,sin(0.366519)=0.358368,而sinf_neon与cosf_neon 都是0.358367。如下
贴出cosf_neon的代码。请问大家有无遇到这种类似的问题?谢谢。
float sinf_neon_hfp(float x) { #ifdef __MATH_NEON asm volatile ( "vld1.32 d3, [%0] \n\t" //d3 = {invrange, range} "vdup.f32 d0, d0[0] \n\t" //d0 = {x, x} "vabs.f32 d1, d0 \n\t" //d1 = {ax, ax} "vmul.f32 d2, d1, d3[0] \n\t" //d2 = d1 * d3[0] "vcvt.u32.f32 d2, d2 \n\t" //d2 = (int) d2 "vmov.i32 d5, #1 \n\t" //d5 = 1 "vcvt.f32.u32 d4, d2 \n\t" //d4 = (float) d2 "vshr.u32 d7, d2, #1 \n\t" //d7 = d2 >> 1 "vmls.f32 d1, d4, d3[1] \n\t" //d1 = d1 - d4 * d3[1] "vand.i32 d5, d2, d5 \n\t" //d5 = d2 & d5 "vclt.f32 d18, d0, #0 \n\t" //d18 = (d0 < 0.0) "vcvt.f32.u32 d6, d5 \n\t" //d6 = (float) d5 "vmls.f32 d1, d6, d3[1] \n\t" //d1 = d1 - d6 * d3[1] "veor.i32 d5, d5, d7 \n\t" //d5 = d5 ^ d7 "vmul.f32 d2, d1, d1 \n\t" //d2 = d1*d1 = {x^2, x^2} "vld1.32 {d16, d17}, [%1] \n\t" //q8 = {p7, p3, p5, p1} "veor.i32 d5, d5, d18 \n\t" //d5 = d5 ^ d18 "vshl.i32 d5, d5, #31 \n\t" //d5 = d5 << 31 "veor.i32 d1, d1, d5 \n\t" //d1 = d1 ^ d5 "vmul.f32 d3, d2, d2 \n\t" //d3 = d2*d2 = {x^4, x^4} "vmul.f32 q0, q8, d1[0] \n\t" //q0 = q8 * d1[0] = {p7x, p3x, p5x, p1x} "vmla.f32 d1, d0, d2[0] \n\t" //d1 = d1 + d0*d2 = {p5x + p7x^3, p1x + p3x^3} "vmla.f32 d1, d3, d1[0] \n\t" //d1 = d1 + d3*d0 = {...., p1x + p3x^3 + p5x^5 + p7x^7} "vmov.f32 s0, s3 \n\t" //s0 = s3 : : "r"(__sinf_rng), "r"(__sinf_lut) : "q0", "q1", "q2", "q3", "q8", "q9" ); #endif } float cosf_neon_hfp(float x) { #ifdef __MATH_NEON float xx = x + M_PI_2; return sinf_neon_hfp(xx); #endif } float cosf_neon_sfp(float x) { #ifdef __MATH_NEON asm volatile ("vdup.f32 d0, r0 \n\t"); cosf_neon_hfp(x); asm volatile ("vmov.f32 r0, s0 \n\t"); #else return cosf_c(x); #endif };
这是算法相关的问题。因为三角函数的实现分为两步:1: range reduction 2: 用Taylor级数计算该三角函数在某个区间内的近似值。
我没有看到上面代码中有range reduction的实现。所以我有一个问题:在上面的Sin函数实现中,支持的Input的范围是多少?
是的。有定义一个
static const float __sinf_rng[2] = { 2.0 / M_PI, M_PI / 2.0 } ALIGN(16);
另外,现在传x+pi/2给sinf_neon 来得到cosf_neon,这样得到的值就与cos(C math)相同了。
如果你的sin实现的输入范围是[2/pi, pi/2], 那你的45行:xx = x + M_PI_2 是不是已经超出范围了?
好的,这个部分我再review这段code,谢谢你。不过我还是蛮纠结为何给sinf_neon 代入 x + pi/2,得到的值就正常了。这部分我再确认一下。谢谢。
不客气。三角函数的实现算法其实比较复杂。如果对实现算法熟悉了,再Debug就应该会比较顺利。