math - TMS320C55X 中的反正切实现

Question

我正在学习 TMS320C55x 中的反正切实现这是源代码:

;* AR0   assigned to _x
;* AR1   assigned to _r
;* T0   assigned to _nx
        PSH     T3
     || BSET    FRCT            ;fractional mode
    SUB #1, T0          ;nx-1
    MOV T0, BRC0        ;repeat nx times
        MOV     #2596 << #16, AC3       ; AC3.Hi = C5
        MOV #-9464 << #16, AC1  ; AC1.Hi = C3
    MOV #32617 << #16, AC2  ; AC2.Hi = C1
*
* Note: loading T3 on the instruction before a multiply that uses it will
* cause a 1-cycle delay.
*
    MPYMR   T3=*AR0+, AC3, AC0      ; (Prime the Pump)
    
     || RPTBLOCAL   loop1-1
            MACR    AC0, T3, AC1, AC0
            MPYR    T3, AC0
          ||MOV *AR0+, T1               ; (for next iteration)
            MACR    AC0, T3, AC2, AC0
            MPYR    T3, AC0
              ||MOV T1, T3
            MOV HI(AC0), *AR1+      ;save result
          ||MPYR    T1, AC3, AC0            ; (for next iteration)
loop1:

        POP     T3
     || BCLR    FRCT            ;return to standard C
    MOV #0, T0          ;return OK value (no possible error)
     || RET

其中 _x 是输入向量，_r 是输出向量。 nx 是元素的数量。问题是关于分配给 AC3、AC1、AC2 的常量。我想这是多项式逼近的系数，但我不明白如何计算它们

Answer 1

我不看汇编代码，但我能猜到那些魔法系数是从哪里来的。

代码注释建议C1、C3、C5是多项式逼近的系数，arctan是一个奇函数，所以它围绕 0 的泰勒展开确实只有 x 的奇次方。将 C1 = 32617 与泰勒展开式中的 1 进行比较 y = x - 1/3 x^3 + 1/5 x^5 - 1/7 x^ 7 + ...，并且给定计算上下文，这进一步表明计算结果按 2^15 = 32768 缩放。

事实证明，y = (32617 x - 9464 x^3 + 2596 x^5)/32768 实际上是 arctan(x) 在区间 [-1, 1] 内。如下所示(在wolfram alpha中验证)近似值的最大绝对误差小于1/1000，并且在端点x = ±1处可以忽略不计y = ±π/4，这在图形计算中可能是可取的。

至于系数的实际推导方式，仅使用 9 个控制点的粗多项式最佳拟合 gives多项式 y = 32613 x - 9443 x^3 + 2573 x^5 其系数已经接近发布代码中使用的系数。更多的控制点和/或额外的条件来最小化端点处的错误会导致系数略有不同，但如果没有任何关于实际使用的优化标准的文档或线索，很难猜测如何与代码中的那些完全匹配.