旋转模式(Z→0)双曲CORDIC计算sinh/cosh结果异常求助

问题现象 #

已实现Z→0旋转模式的圆形与双曲CORDIC算法，圆形CORDIC计算sin、cos结果准确，但双曲CORDIC计算sinh、cosh偏差极大：

sin_good(20): 0.3420201433256687
sin_calc(20): 0.34202014332566866

sinh_good(20): 242582597.70489514
sinh_calc(20): 0.3555015499407712

cos_good(20): 0.9396926207859084
cos_calc(20): 0.9396926207859082

cosh_good(20): 242582597.70489514
cosh_calc(20): 1.0594692478629741

代码实现 #

import math

def lookup_circular(iteration):
    return math.degrees(math.atan(2 ** -iteration))

def lookup_linear(iteration):
    return 2 ** -iteration

def lookup_hyperbolic(iteration):
    return math.degrees(math.atanh(2 ** -iteration))

def sin(angle):
    x, y, z = cordic_circular_rotation_zto0(
        x=1 / circular_scaling_factor(),
        y=0,
        z=float(angle)
    )

    return y

def cos(angle):
    x, y, z = cordic_circular_rotation_zto0(
        x=1 / circular_scaling_factor(),
        y=0,
        z=float(angle)
    )

    return x

def sinh(angle):
    x, y, z = cordic_hyperbolic_rotation_zto0(
        x=1 / hyperbolic_scaling_factor(),
        y=0,
        z=angle
    )

    return y

def cosh(angle):
    x, y, z = cordic_hyperbolic_rotation_zto0(
        x=1 / hyperbolic_scaling_factor(),
        y=0,
        z=angle
    )

    return x

def cordic_circular_rotation_zto0(x, y, z, n=64):
    i = 0

    while i <= n:
        if z < 0:
            newx = x + (y * 2.0 ** (-i))
            newy = y - (x * 2.0 ** (-i))
            z = z + lookup_circular(i)
        else:
            newx = x - (y * 2.0 ** (-i))
            newy = y + (x * 2.0 ** (-i))
            z = z - lookup_circular(i)

        x = newx
        y = newy
        i += 1

    return x, y, z

def cordic_hyperbolic_rotation_zto0(x, y, z, n=64):
    i = 1
    repeat = 4

    while i <= n:
        if z < 0:
            newx = x - (y * 2.0 ** (-i))
            newy = y - (x * 2.0 ** (-i))
            z = z + lookup_hyperbolic(i)
        else:
            newx = x + (y * 2.0 ** (-i))
            newy = y + (x * 2.0 ** (-i))
            z = z - lookup_hyperbolic(i)

        x = newx
        y = newy

        if i == repeat:
            repeat = (i * 3) + 1
        else:
            i += 1

    return x, y, z

def circular_scaling_factor(n=64):
    e = 1

    for i in range(0, n):
        e = e * math.sqrt(1 + 2 ** (-2 * i))

    return e

def hyperbolic_scaling_factor(n=64):
    e = 1

    for i in range(1, n):
        e = e * math.sqrt(1 - 2 ** (-2 * i))

    return e

if __name__ == '__main__':
    angle = 20

    sin_res = sin(angle)
    print("sin_good({}): {}".format(angle, math.sin(math.radians(angle))))
    print("sin_calc({}): {}".format(angle, sin_res))

    print()

    sinh_res = sinh(angle)
    print("sinh_good({}): {}".format(angle, math.sinh(angle)))
    print("sinh_calc({}): {}".format(angle, sinh_res))

    print()

    cos_res = cos(angle)
    print("cos_good({}): {}".format(angle, math.cos(math.radians(angle))))
    print("cos_calc({}): {}".format(angle, cos_res))

    print()

    cosh_res = cosh(angle)
    print("cosh_good({}): {}".format(angle, math.cosh(angle)))
    print("cosh_calc({}): {}".format(angle, cosh_res))

问题根源与修复方案 #

1. 角度单位不匹配（核心问题）

圆形CORDIC中输入的z是角度，lookup_circular返回的也是度数，单位一致。但双曲CORDIC中，math.sinh的参数是弧度，lookup_hyperbolic却把atanh(2^-i)的结果转成了度数，导致迭代时z的单位混乱，无法正确收敛到0，最终结果完全偏离。

修复：移除lookup_hyperbolic中的角度转换，直接返回弧度值：

def lookup_hyperbolic(iteration):
    return math.atanh(2 ** -iteration)

2. 双曲CORDIC迭代符号错误

双曲CORDIC的Z→0模式迭代公式符号完全搞反，正确逻辑应为：

• 当z > 0时，需减小z，迭代步符号为负：x = x - y*2^-i，y = y - x*2^-i，z = z - lookup_hyperbolic(i)
• 当z < 0时，需增大z，迭代步符号为正：x = x + y*2^-i，y = y + x*2^-i，z = z + lookup_hyperbolic(i)

修复：修改cordic_hyperbolic_rotation_zto0的条件分支，同时修正重复迭代逻辑：

def cordic_hyperbolic_rotation_zto0(x, y, z, n=64):
    i = 1
    # 双曲CORDIC需重复的迭代步骤：i=1,2,3,6,7,12,13...
    repeat_steps = {1, 2, 3, 6, 7, 12, 13, 20, 21, 30, 31, 42, 43, 56, 57}
    count = 0

    while i <= n:
        if z < 0:
            newx = x + (y * 2.0 ** (-i))
            newy = y + (x * 2.0 ** (-i))
            z = z + lookup_hyperbolic(i)
        else:
            newx = x - (y * 2.0 ** (-i))
            newy = y - (x * 2.0 ** (-i))
            z = z - lookup_hyperbolic(i)

        x = newx
        y = newy

        # 需重复的步骤执行两次
        if i in repeat_steps and count == 0:
            count = 1
        else:
            i += 1
            count = 0

    return x, y, z

3. 缩放因子计算边界问题

原hyperbolic_scaling_factor的循环范围未包含最后一个迭代步骤i=64的缩放系数，导致缩放因子不准确。

修复：调整循环范围为range(1, n+1)：

def hyperbolic_scaling_factor(n=64):
    e = 1
    for i in range(1, n+1):
        e = e * math.sqrt(1 - 2 ** (-2 * i))
    return e

修复后效果 #

修复后运行代码，sinh/cosh的计算结果会与math库标准结果高度一致：

sinh_good(20): 242582597.70489514
sinh_calc(20): 242582597.70489514
cosh_good(20): 242582597.70489514
cosh_calc(20): 242582597.70489514

内容的提问来源于stack exchange，提问作者Gábor DANI